--- a
+++ b/partyMod/src/Utils.c
@@ -0,0 +1,867 @@
+
+/**
+    Some commonly needed utility functions.
+    *\file Utils.c
+    *\author $Author$
+    *\date $Date$
+*/
+                
+#include "party.h"
+                
+                
+/**
+    Computes the Kronecker product of two matrices\n
+    *\param A matrix
+    *\param m nrow(A)
+    *\param n ncol(A)
+    *\param B matrix
+    *\param r nrow(B)
+    *\param s ncol(B)
+    *\param ans return value; a pointer to a REALSXP-vector of length (mr x ns)
+*/
+
+void C_kronecker (const double *A, const int m, const int n,
+                  const double *B, const int r, const int s,
+                  double *ans) {
+
+    int i, j, k, l, mr, js, ir;
+    double y;
+
+    mr = m * r;
+    for (i = 0; i < m; i++) {
+        ir = i * r;
+        for (j = 0; j < n; j++) {
+            js = j * s;
+            y = A[j*m + i];
+            for (k = 0; k < r; k++) {
+                for (l = 0; l < s; l++) {
+                    ans[(js + l) * mr + ir + k] = y * B[l * r + k];
+                }
+            }
+        }
+    }
+}  
+
+
+/**
+    R-interface to C_kronecker
+    *\param A matrix
+    *\param B matrix
+*/
+                
+SEXP R_kronecker (SEXP A, SEXP B) {
+
+    /*  The Kronecker product, a real (mr x ns) matrix */
+    SEXP ans; 
+    int *adim, *bdim;
+
+    if (!isReal(A) || !isReal(B)) 
+        error("R_kronecker: A and B are not of type REALSXP");
+
+    if (isMatrix(A)) {
+        adim = INTEGER(getAttrib(A, R_DimSymbol));
+    } else {
+        /* assume row vectors */
+        adim = Calloc(2, int);
+        adim[0] = 1;
+        adim[1] = LENGTH(A);
+    }
+    
+    if (isMatrix(B)) {
+        bdim = INTEGER(getAttrib(B, R_DimSymbol));
+    } else {
+        /* assume row vectors */
+        bdim = Calloc(2, int);
+        bdim[0] = 1;
+        bdim[1] = LENGTH(B);
+    }
+
+    PROTECT(ans = allocMatrix(REALSXP, 
+                              adim[0] * bdim[0], 
+                              adim[1] * bdim[1]));
+    C_kronecker(REAL(A), adim[0], adim[1], 
+                REAL(B), bdim[0], bdim[1], REAL(ans));
+    if (!isMatrix(A)) Free(adim); 
+    if (!isMatrix(B)) Free(bdim);
+    UNPROTECT(1);
+    return(ans);
+}
+
+
+/**
+    C- and R-interface to La_svd 
+    *\param jobu
+    *\param jobv
+    *\param x
+    *\param s
+    *\param u
+    *\param v
+    *\param method
+*/
+
+void CR_La_svd(int dim, SEXP jobu, SEXP jobv, SEXP x, SEXP s, SEXP u, SEXP v,
+               SEXP method)
+{
+    int *xdims, n, p, lwork, info = 0;
+    int *iwork;
+    double *work, *xvals, tmp;
+    /* const char * meth; not used*/
+
+    if (!(isString(jobu) && isString(jobv)))
+	error(("'jobu' and 'jobv' must be character strings"));
+    if (!isString(method))
+	error(("'method' must be a character string"));
+    /* meth = CHAR(STRING_ELT(method, 0)); not used */
+    xdims = INTEGER(coerceVector(getAttrib(x, R_DimSymbol), INTSXP));
+    n = xdims[0]; p = xdims[1];
+    xvals = Calloc(n * p, double);
+    /* work on a copy of x */
+    Memcpy(xvals, REAL(x), (size_t) (n * p));
+
+    {
+	int ldu = INTEGER(getAttrib(u, R_DimSymbol))[0],
+	    ldvt = INTEGER(getAttrib(v, R_DimSymbol))[0];
+        /* this only works for SQUARE matrices, potentially
+           with reduced dimension, so use dim x dim for input and
+           output */
+        ldu = dim;
+        ldvt = dim;
+	iwork= (int *) Calloc(8*(n<p ? n : p), int);
+
+	/* ask for optimal size of work array */
+	lwork = -1;
+	F77_CALL(dgesdd)(CHAR(STRING_ELT(jobu, 0)),
+/* was			 &n, &p, xvals, &n, REAL(s), for the non-square case */
+			 &dim, &dim, xvals, &dim, REAL(s),
+			 REAL(u), &ldu,
+			 REAL(v), &ldvt,
+			 &tmp, &lwork, iwork, &info);
+	if (info != 0)
+	    error(("error code %d from Lapack routine '%s'"), info, "dgesdd");
+	lwork = (int) tmp;
+	work = Calloc(lwork, double);
+	F77_CALL(dgesdd)(CHAR(STRING_ELT(jobu, 0)),
+			 &dim, &dim, xvals, &dim, REAL(s),
+/* was			 &n, &p, xvals, &n, REAL(s), for the non-square case */
+			 REAL(u), &ldu,
+			 REAL(v), &ldvt,
+			 work, &lwork, iwork, &info);
+	if (info != 0)
+	    error(("error code %d from Lapack routine '%s'"), info, "dgesdd");
+    }
+    Free(work); Free(xvals); Free(iwork);
+}
+
+/**
+    C-interface to CR_La_svd
+    *\param x matrix
+    *\param svdmem an object of class `svd_mem'
+*/
+
+void C_svd (SEXP x, SEXP svdmem) {
+
+    int p, i;
+    double *du, *dv;
+
+    if (!isMatrix(x) || !isReal(x))
+        error("x is not a real matrix");
+
+    du = REAL(GET_SLOT(svdmem, PL2_uSym));
+    dv = REAL(GET_SLOT(svdmem, PL2_vSym));
+    p = INTEGER(GET_SLOT(svdmem, PL2_pSym))[0];
+    if (nrow(x) < p) error("svd p x error");
+    for (i = 0; i < p*p; i++) {
+        du[i] = 0.0;
+        dv[i] = 0.0;
+    }
+    CR_La_svd(p, GET_SLOT(svdmem, PL2_jobuSym), 
+        GET_SLOT(svdmem, PL2_jobvSym), x, GET_SLOT(svdmem, PL2_sSym), 
+        GET_SLOT(svdmem, PL2_uSym), GET_SLOT(svdmem, PL2_vSym), 
+        GET_SLOT(svdmem, PL2_methodSym));
+    /* return(R_NilValue); */
+}
+
+
+/**
+    R-interface to CR_La_svd
+    *\param x matrix
+    *\param svdmem an object of class `svd_mem'
+*/
+
+SEXP R_svd (SEXP x, SEXP svdmem) {
+
+    C_svd(x, svdmem);
+    return(R_NilValue);
+}
+
+
+/**
+    Reorder the linear statistic, expectation, and covariance
+    in a way that elements with zero variance come last. These
+    will be ignored in later computation of the MP inverse
+    *\param x an object of class `LinStatExpectCovarMPinv'
+*/
+
+void C_linexpcovReduce (SEXP x) {
+
+    double *dlinstat, *dexp, *dcov;
+    double *dlinstat2, *dexp2, *dcov2;
+    int pq, pqn, *zerovar, i, j, itmp, jtmp, sumzv = 0, *dim;
+
+    /* the statistic, expectation, covariance in original dimension */    
+    pq = INTEGER(GET_SLOT(x, PL2_dimensionSym))[0];
+    dlinstat = REAL(GET_SLOT(x, PL2_linearstatisticSym));
+    dexp = REAL(GET_SLOT(x, PL2_expectationSym));
+    dcov = REAL(GET_SLOT(x, PL2_covarianceSym));
+
+    /* indicator of zero variance */
+    zerovar = Calloc(pq, int);
+
+    /* identify and count zero variances (we can use 0.0 because variances
+       corresponding to empty levels are exactly zero */
+    for (i = 0; i < pq; i++) {
+        if (dcov[i + i * pq] > 0.0) {
+            zerovar[i] = 0;
+        } else {
+            zerovar[i] = 1;
+            sumzv++;
+        }
+    }
+
+    /* if there is any such element and not all variances are zero*/
+    if ((sumzv > 0) & (sumzv < pq)) {
+
+        /* do we really need a copy ? */
+        pqn = pq - sumzv;
+        dlinstat2 = Calloc(pq, double);
+        dexp2 = Calloc(pq, double);
+        dcov2 = Calloc(pq * pq, double);
+        
+        /* init */
+        for (i = 0; i < pq; i++) {
+            dlinstat2[i] = 0.0; 
+            dexp2[i] = 0.0; 
+            for (j = 0; j < pq; j++)
+                dcov2[i + j*pq] = 0.0;
+        }
+        
+        /* overwrite zero variance elements with subsequent elements */
+        itmp = 0;
+        for (i = 0; i < pq; i++) {
+            if (zerovar[i] == 0) {
+                dlinstat2[itmp] = dlinstat[i];
+                dexp2[itmp] = dexp[i];
+                jtmp = 0;
+                for (j = 0; j < pq; j++) {
+                    if (zerovar[j] == 0) {
+                        dcov2[itmp + jtmp * pqn] = dcov[i + j * pq];
+                        jtmp++;
+                    }
+                }
+                itmp++;
+            }
+        }
+                                        
+        for (i = 0; i < pq; i++) {
+            dlinstat[i] = dlinstat2[i]; 
+            dexp[i] = dexp2[i]; 
+            for (j = 0; j < pq; j++)
+                dcov[i + j*pq] = dcov2[i + j * pq];
+        }
+
+        /* ATTENTION: This is dangerous but the only way to tell
+           svd and friends that we want to use the first pqn elements only! */
+        dim = INTEGER(GET_SLOT(x, PL2_dimensionSym));
+        dim[0] = pqn;
+        /* we reset the original dimension in C_TestStatistic */
+        
+        Free(dlinstat2);
+        Free(dexp2);
+        Free(dcov2);
+    }
+    Free(zerovar);
+}
+
+
+/**
+    R-interface to C_linexpcovReduce
+    *\param x an object of class `LinStatExpectCovarMPinv'
+*/
+
+SEXP R_linexpcovReduce (SEXP x) {
+
+    C_linexpcovReduce(x);
+    return(R_NilValue);
+}
+        
+
+/**
+    Moore-Penrose inverse of a matrix
+    *\param x matrix
+    *\param tol a tolerance bound
+    *\param svdmem an object of class `svd_mem'
+    *\param ans return value; an object of class `ExpectCovarMPinv'
+*/
+
+void C_MPinv (SEXP x, double tol, SEXP svdmem, SEXP ans) {
+
+    SEXP d, u, vt;
+    int i, j, p, pn, k, *positive;
+    double *dd, *du, *dvt, *dMPinv;
+    double *drank;
+    
+    drank = REAL(GET_SLOT(ans, PL2_rankSym));
+    dMPinv = REAL(GET_SLOT(ans, PL2_MPinvSym));
+
+    C_svd(x, svdmem);
+    d = GET_SLOT(svdmem, PL2_sSym);
+    dd = REAL(d);
+    u = GET_SLOT(svdmem, PL2_uSym);
+    du = REAL(u);
+    vt = GET_SLOT(svdmem, PL2_vSym);
+    dvt = REAL(vt);
+    /* this may be the reduced dimension! Use the first p elements only!!!*/
+    p = INTEGER(GET_SLOT(svdmem, PL2_pSym))[0];
+
+    if (tol * dd[0] > tol) tol = tol * dd[0];
+
+    positive = Calloc(p, int); 
+    
+    drank[0] = 0.0;
+    for (i = 0; i < p; i++) {
+        if (dd[i] > tol) {
+            positive[i] = 1;
+            drank[0] += 1.0;
+        } 
+    }
+    
+    for (j = 0; j < p; j++) {
+        if (positive[j]) {
+            for (i = 0; i < p; i++)
+                du[j * p + i] *= (1 / dd[j]);
+        }
+    }
+    
+    for (i = 0; i < p; i++) {
+        for (j = 0; j < p; j++) {
+            dMPinv[j * p + i] = 0.0;
+            for (k = 0; k < p; k++) {
+                if (positive[k])
+                    dMPinv[j * p + i] += dvt[i * p + k] * du[p * k + j];
+            }
+        }
+    }
+
+    Free(positive);
+}
+
+
+/**
+    R-interface to C_MPinv 
+    *\param x matrix
+    *\param tol a tolerance bound
+    *\param svdmem an object of class `svd_mem'
+*/
+
+SEXP R_MPinv (SEXP x, SEXP tol, SEXP svdmem) {
+
+    SEXP ans;
+    int p;
+
+    if (!isMatrix(x) || !isReal(x))
+        error("R_MPinv: x is not a real matrix");
+
+    if (nrow(x) != ncol(x)) 
+        error("R_MPinv: x is not a square matrix");
+
+    if (!isReal(tol) || LENGTH(tol) != 1)
+        error("R_MPinv: tol is not a scalar real");
+
+    p = nrow(x);
+    /* potentially, the effective dimension was reduced
+    if (p != INTEGER(GET_SLOT(svdmem, PL2_pSym))[0])
+        error("R_MPinv: dimensions don't match");
+    */
+
+    PROTECT(ans = NEW_OBJECT(MAKE_CLASS("LinStatExpectCovarMPinv")));
+    SET_SLOT(ans, PL2_MPinvSym, PROTECT(allocMatrix(REALSXP, p, p)));
+    SET_SLOT(ans, PL2_rankSym, PROTECT(allocVector(REALSXP, 1)));
+    
+    C_MPinv(x, REAL(tol)[0], svdmem, ans);
+    
+    UNPROTECT(3);
+    return(ans);
+}
+
+
+/**
+    the maximum of a double vector
+    *\param x vector
+    *\param n its length
+*/
+
+
+double C_max(const double *x, const int n) {
+   double tmp = 0.0;
+   int i;
+   
+   for (i = 0; i < n; i++) {
+       if (x[i] > tmp) tmp = x[i];
+   }
+   return(tmp);
+}
+
+
+/**
+    R-interface to C_max
+    *\param x numeric vector
+*/
+
+SEXP R_max(SEXP x) {
+
+    SEXP ans;
+    int n;
+    
+    if (!isReal(x)) 
+        error("R_max: x is not of type REALSXP");
+    n = LENGTH(x);
+    PROTECT(ans = allocVector(REALSXP, 1));
+    REAL(ans)[0] = C_max(REAL(x), n);
+    UNPROTECT(1);
+    return(ans);
+}
+
+
+/**
+    absolute value 
+    *\param x numeric vector
+    *\param n length(x)
+*/
+
+void C_abs(double *x, int n) {
+
+    int i;
+    for (i = 0; i < n; i++) x[i] = fabs(x[i]);
+}
+
+
+/**
+    R-interface to C_abs
+    *\param x numeric vector
+*/
+
+SEXP R_abs(SEXP x) {
+
+    SEXP ans;
+    int n;
+    
+    if (!isReal(x)) 
+        error("R_max: x is not of type REALSXP");
+    n = LENGTH(x);
+    PROTECT(ans = duplicate(x));
+    C_abs(REAL(ans), n);
+    UNPROTECT(1);
+    return(ans);
+}
+
+
+/**
+    matrix product x %*% y
+    *\param x a matrix
+    *\param nrx number of rows of x
+    *\param ncx number of cols of x
+    *\param y a matrix
+    *\param nry number of rows of y
+    *\param ncy number of cols of y
+    *\param z a matrix of dimension nrx x ncy
+*/
+
+void C_matprod(double *x, int nrx, int ncx,
+               double *y, int nry, int ncy, double *z)
+{
+    const char *transa = "N", *transb = "N";
+    double one = 1.0, zero = 0.0;
+    int i;
+
+    if (nrx > 0 && ncx > 0 && nry > 0 && ncy > 0) {
+        F77_CALL(dgemm)(transa, transb, &nrx, &ncy, &ncx, &one,
+	                x, &nrx, y, &nry, &zero, z, &nrx);
+    } else /* zero-extent operations should return zeroes */
+	for(i = 0; i < nrx*ncy; i++) z[i] = 0;
+}
+
+
+/**
+    R-interface to C_matprod
+    *\param x a matrix
+    *\param y a matrix
+*/
+
+SEXP R_matprod(SEXP x, SEXP y) {
+
+    SEXP ans;
+    
+    int nrx, ncx, nry, ncy;
+    
+    nrx = nrow(x);
+    ncx = ncol(x);
+    nry = nrow(y);
+    ncy = ncol(y);
+
+    if (ncx != nry)
+        error("R_matprod: dimensions don't match");
+    PROTECT(ans = allocMatrix(REALSXP, nrx, ncy));
+    C_matprod(REAL(x), nrx, ncx, REAL(y), nry, ncy, REAL(ans));
+    UNPROTECT(1);
+    return(ans);
+}
+
+
+/**
+    matrix product x %*% t(y)
+    *\param x a matrix
+    *\param nrx number of rows of x
+    *\param ncx number of cols of x
+    *\param y a matrix
+    *\param nry number of rows of y
+    *\param ncy number of cols of y
+    *\param z a matrix of dimension nrx x ncy
+*/
+
+void C_matprodT(double *x, int nrx, int ncx,
+                double *y, int nry, int ncy, double *z)
+{
+    const char *transa = "N", *transb = "T";
+    double one = 1.0, zero = 0.0;
+    int i;
+
+    if (nrx > 0 && ncx > 0 && nry > 0 && ncy > 0) {
+        F77_CALL(dgemm)(transa, transb, &nrx, &nry, &ncy, &one,
+	                x, &nrx, y, &nry, &zero, z, &nrx);
+    } else /* zero-extent operations should return zeroes */
+	for(i = 0; i < nrx*nry; i++) z[i] = 0;
+}
+
+
+/**
+    R-interface to C_matprodT
+    *\param x a matrix
+    *\param y a matrix
+*/
+
+SEXP R_matprodT(SEXP x, SEXP y) {
+
+    SEXP ans;
+    int nrx, ncx, nry, ncy;
+    
+    nrx = nrow(x);
+    ncx = ncol(x);
+    nry = nrow(y);
+    ncy = ncol(y);
+
+    if (ncx != ncy)
+        error("R_matprod: dimensions don't match");
+    PROTECT(ans = allocMatrix(REALSXP, nrx, nry));
+    C_matprodT(REAL(x), nrx, ncx, REAL(y), nry, ncy, REAL(ans));
+    UNPROTECT(1);
+    return(ans);
+}
+
+
+/**
+    compute a permutation of a (random subset of) 0:(m-1)
+    *\param x an integer vector of length m
+    *\param m integer
+    *\param k integer
+    *\param ans an integer vector of length k
+*/    
+
+void C_SampleNoReplace(int *x, int m, int k, int *ans) {
+     
+    int i, j, n = m;
+
+    for (i = 0; i < m; i++)
+        x[i] = i;
+    for (i = 0; i < k; i++) {
+        j = floor((double) n * unif_rand()); 
+        ans[i] = x[j];
+        x[j] = x[--n];  
+    }
+}
+
+
+/**
+    R-interface to C_SampleNoReplace: the permutation case
+    *\param m integer
+*/    
+
+SEXP R_permute(SEXP m) {
+    
+    SEXP x, ans;
+    int n;
+    
+    n = INTEGER(m)[0];
+    PROTECT(x = allocVector(INTSXP, n));
+    PROTECT(ans = allocVector(INTSXP, n));
+    C_SampleNoReplace(INTEGER(x), n, n, INTEGER(ans));
+    UNPROTECT(2);
+    return(ans);
+}
+
+
+/**
+    R-interface to C_SampleNoReplace: the subset case
+    *\param m integer
+    *\param k integer
+*/    
+
+SEXP R_rsubset(SEXP m, SEXP k) {
+    
+    SEXP x, ans;
+    int n, j;
+    
+    n = INTEGER(m)[0];
+    j = INTEGER(k)[0];
+    PROTECT(x = allocVector(INTSXP, n));
+    PROTECT(ans = allocVector(INTSXP, j));
+    C_SampleNoReplace(INTEGER(x), n, j, INTEGER(ans));
+    UNPROTECT(2);
+    return(ans);
+}
+
+/* Unequal probability sampling; without-replacement case */
+
+void C_ProbSampleNoReplace(int n, double *p, int *perm,
+                           int nans, int *ans)
+{
+    double rT, mass, totalmass;
+    int i, j, k, n1;
+
+    /* Record element identities */
+    for (i = 0; i < n; i++)
+	perm[i] = i + 1;
+
+    /* Sort probabilities into descending order */
+    /* Order element identities in parallel */
+    revsort(p, perm, n);
+
+    /* Compute the sample */
+    totalmass = 1;
+    for (i = 0, n1 = n-1; i < nans; i++, n1--) {
+	rT = totalmass * unif_rand();
+	mass = 0;
+	for (j = 0; j < n1; j++) {
+	    mass += p[j];
+	    if (rT <= mass)
+		break;
+	}
+	ans[i] = perm[j];
+	totalmass -= p[j];
+	for(k = j; k < n1; k++) {
+	    p[k] = p[k + 1];
+	    perm[k] = perm[k + 1];
+	}
+    }
+}
+
+
+/**
+    determine if i is element of the integer vector set
+    *\param i an integer
+    *\param iset a pointer to an integer vector
+    *\param p length(iset)
+*/
+
+int i_in_set(int i, int *iset, int p) {
+
+    int j, is = 0;
+        
+    if (p == 0) return(0);
+                    
+    for (j = 0; j < p; j++) {
+        if (iset[j] == i) {  
+            is = 1;
+            break; 
+        }
+    }
+    return(is);
+}
+
+int C_i_in_set(int i, SEXP set) {
+    if (LENGTH(set) > 0)
+        return(i_in_set(i, INTEGER(set), LENGTH(set)));
+    else 
+        return(0);
+}
+    
+int nrow(SEXP x) {
+    return(INTEGER(getAttrib(x, R_DimSymbol))[0]);
+}
+
+int ncol(SEXP x) {
+    return(INTEGER(getAttrib(x, R_DimSymbol))[1]);
+}
+
+/* compute index of variable with smallest p-value 
+   (and largest test statistic in case two or more p-values coincide -- 
+    should not happen anymore since we use 1 - (1 - p)^k for Bonferroni adjustment)
+*/
+int C_whichmax(double *pvalue, double *teststat, int ninputs) {
+
+    int ans = -1, j;
+    double tmppval = 0.0, tmptstat = 0.0;
+       
+    /* <FIXME> can we switch to the log scale here? </FIXME> */
+
+    tmppval = 0.0;
+    tmptstat = 0.0;
+    for (j = 0; j < ninputs; j++) {
+        if (pvalue[j] > tmppval) {
+            ans = j;
+            tmppval = pvalue[j];
+            tmptstat = teststat[j];
+        } else {
+            if (pvalue[j] == tmppval && teststat[j] > tmptstat) {  
+                ans = j;
+                tmppval = pvalue[j];
+                tmptstat = teststat[j];
+            }
+        }
+    }
+    return(ans);
+}
+
+SEXP R_whichmax(SEXP x, SEXP y) {
+    SEXP ans;
+    
+    if (LENGTH(x) != LENGTH(y)) error("different length");
+    PROTECT(ans = allocVector(INTSXP, 1));
+    INTEGER(ans)[0] = C_whichmax(REAL(x), REAL(y), LENGTH(x));
+    UNPROTECT(1);
+    return(ans);
+}
+
+SEXP R_listplus(SEXP a, SEXP b, SEXP which) {
+
+    int na, nb, i, j, *iwhich;
+    double *dae, *dbe;
+    SEXP ae, be;
+
+    na = LENGTH(a);
+    nb = LENGTH(b);
+    if (na != nb) error("a and b are of different length");
+    
+    iwhich = LOGICAL(which);
+    
+    for (i = 0; i < na; i++) {
+        if (iwhich[i]) continue;
+        
+        ae = VECTOR_ELT(a, i);
+        be = VECTOR_ELT(b, i);
+
+        if (LENGTH(ae) != LENGTH(be)) 
+            error("elements %d are of different length", i);
+            
+        if (!isReal(ae) || !isReal(be))
+            error("elements %d are not of type double", i);
+            
+        dae = REAL(ae);
+        dbe = REAL(be);
+        for (j = 0; j < LENGTH(ae); j++) 
+            dae[j] += dbe[j];
+    }
+    return(a);
+}
+
+SEXP R_modify_response(SEXP x, SEXP vf) {
+
+    double *src, *tar;
+    int i, n;
+    
+    src = REAL(x);
+    n = LENGTH(x);
+
+    tar = REAL(get_transformation(vf, 1));
+    for (i = 0; i < n; i++)
+        tar[i] = src[i];
+
+    tar = REAL(get_test_trafo(vf));
+    for (i = 0; i < n; i++)
+        tar[i] = src[i];
+
+    tar = REAL(get_predict_trafo(vf));
+    for (i = 0; i < n; i++)
+        tar[i] = src[i];
+
+    tar = REAL(get_variable(vf, 1));
+    for (i = 0; i < n; i++)
+        tar[i] = src[i];
+                                          
+    return(R_NilValue);
+}
+
+double F77_SUB(unifrnd)(void) { return unif_rand(); }
+
+void C_SampleSplitting(int n, double *prob, int *weights, int k) {
+
+    int i;
+    double *tmpprob;
+    int *ans, *perm;
+
+    tmpprob = Calloc(n, double);
+    perm = Calloc(n, int);
+    ans = Calloc(k, int);
+    for (i = 0; i < n; i++) tmpprob[i] = prob[i];
+
+    C_ProbSampleNoReplace(n, tmpprob, perm, k, ans);
+    for (i = 0; i < n; i++) weights[i] = 0;
+    for (i = 0; i < k; i++)
+        weights[ans[i] - 1] = 1;
+    Free(tmpprob); Free(perm); Free(ans);
+}
+
+/**
+    Remove weights vector from each node of a tree (in order to save memory)
+    \*param subtree a tree
+*/ 
+
+void C_remove_weights(SEXP subtree, int removestats) {
+
+    SET_VECTOR_ELT(subtree, S3_WEIGHTS, R_NilValue);
+    
+    if (!S3get_nodeterminal(subtree)) {
+        if (removestats) {
+            SET_VECTOR_ELT(VECTOR_ELT(subtree, S3_CRITERION), 
+                           S3_iCRITERION, R_NilValue);
+            SET_VECTOR_ELT(VECTOR_ELT(subtree, S3_CRITERION), 
+                           S3_STATISTICS, R_NilValue);
+        }
+        C_remove_weights(S3get_leftnode(subtree), removestats);
+        C_remove_weights(S3get_rightnode(subtree), removestats);
+    }
+}
+
+SEXP R_remove_weights(SEXP subtree, SEXP removestats) {
+
+    C_remove_weights(subtree, LOGICAL(removestats)[0]);
+    return(R_NilValue);
+}
+
+double* C_tempweights(int j, SEXP weights, SEXP fitmem, SEXP inputs) {
+
+    int nobs, *iNAs, i, k;
+    double *dw, *dweights;
+    SEXP NAs;
+    
+    dw = REAL(get_weights(fitmem));
+    nobs = LENGTH(weights);
+    dweights = REAL(weights);
+    NAs = get_missings(inputs, j);
+    iNAs = INTEGER(NAs);
+    if (length(NAs) == 0) return(dw);
+    for (i = 0; i < nobs; i++) dw[i] = dweights[i];
+    for (k = 0; k < LENGTH(NAs); k++)
+        dw[iNAs[k] - 1] = 0.0;
+    
+    return(dw);
+}