--- a
+++ b/partyMod/src/Splits.c
@@ -0,0 +1,345 @@
+
+/**
+    Binary splits 
+    *\file Splits.c
+    *\author $Author$
+    *\date $Date$
+*/
+                
+#include "party.h"
+
+
+/**
+    Search for a cutpoint in a ordered variable x maximizing a two-sample
+    statistic w.r.t. (the influence function of ) the response variable y.
+    *\param x raw numeric measurements 
+    *\param p dimension of the transformation
+    *\param y values of the influence function
+    *\param q dimension of the influence function
+    *\param weights case weights
+    *\param n number of observations
+    *\param orderx the ordering of the transformations, i.e. R> order(x)
+    *\param splitctrl an object of class `SplitControl'
+    *\param linexpcov2sample an (uninitialized) object of class 
+                             `LinStatExpectCovar' with p = 1
+    *\param expcovinf an initialized object of class `ExpectCovarInfluence'
+    *\param cutpoint return value; pointer to a double for the cutpoint in x
+    *\param maxstat return value; pointer to a double for the 
+                    maximal test statistic
+    *\param statistics return value; pointer to a n-dim double 
+                       for the statistics
+*/
+                                
+void C_split(const double *x, int p,
+             const double *y, int q,
+             const double *weights, int n,
+             const int *orderx,
+             SEXP splitctrl, SEXP linexpcov2sample, 
+             SEXP expcovinf, double *cutpoint, double *maxstat, 
+             double *statistics) {
+
+    double *dExp_y, *dCov_y, *dlinstat, *dexpect, *dcovar, 
+           tol, sweights, minprob, minbucket, w, tx, f1, f2, f1w, f2ww, tmp;
+    double minobs, maxobs, xmax;
+    int lastj, i, j, k;
+
+    if (p != 1) error("C_split: p not equal to one");
+    tol = get_tol(splitctrl);
+
+    /* init statistics and determine the maximal value with positive weight 
+       since we can't choose this one as cutpoint
+    */
+    xmax = 0.0;
+    for (i = 0; i < n; i++) {
+        statistics[i] = 0.0;
+        if (weights[i] > 0.0 && x[i] > xmax) xmax = x[i];
+    }
+
+    /* we already have expecation and covariance of the response
+     * values and the sum of the weights */
+    dExp_y = REAL(GET_SLOT(expcovinf, PL2_expectationSym));
+    dCov_y = REAL(GET_SLOT(expcovinf, PL2_covarianceSym));
+    sweights = REAL(GET_SLOT(expcovinf, PL2_sumweightsSym))[0];
+
+    /* if there is something to split */
+    if (sweights > 1) {
+
+        /* we need to ensure that at least minbucket weights 
+           are there to split (either left or right) */
+        minprob = get_minprob(splitctrl);
+        minbucket = get_minbucket(splitctrl);
+        minobs = sweights * minprob + 1.0;
+
+        if (minobs < minbucket) 
+            minobs = minbucket; 
+        maxobs = sweights * (1 - minprob) - 1.0;
+        if (maxobs > sweights - minbucket) 
+            maxobs = sweights - minbucket; 
+
+        f1 = (double) sweights / (sweights - 1);
+        f2 = 1.0 / (sweights - 1);
+        w = 0.0;
+    
+        /* pointers to the R-objects */
+        dlinstat = REAL(GET_SLOT(linexpcov2sample, PL2_linearstatisticSym));
+        for (k = 0; k < q; k++) dlinstat[k] = 0.0;
+        dexpect = REAL(GET_SLOT(linexpcov2sample, PL2_expectationSym));
+        dcovar = REAL(GET_SLOT(linexpcov2sample, PL2_covarianceSym));
+
+        tx = 0.0;
+        lastj = 0;
+
+        /* for all possible cutpoints (defined by the observations x) */
+        for (i = 0; i < (n - 1); i++) {
+    
+            /* the ordering of the ith observation */
+            j = orderx[i] - 1;
+        
+            /* if the corresponding weight is zero */
+            if (weights[j] == 0.0) continue;
+
+            /* just a check: can be removed later */
+            if (w > 0 && x[j] < tx)
+                warning("C_split: inconsistent ordering: %f < %f!\n", 
+                        x[j], tx);
+        
+            /* handle ties: delete the entry of the last visited observation
+               (take care of zero weights!) */
+            if (w > 0 && x[j] == tx)
+                statistics[lastj] = 0.0; 
+
+            /* store the value and position of the j smallest observation */
+            tx = x[j];
+            lastj = j;
+        
+            w += weights[j];
+
+            /* do not consider those splits */
+            if (w > maxobs || x[j] >= xmax) break;
+
+            /* compute the linear statistic and expectation and 
+             * covariance if needed */
+            for (k = 0; k < q; k++)
+                dlinstat[k] += y[n * k + j] * weights[j];
+ 
+            if (w >= minobs) {
+                for (k = 0; k < q; k++)
+                    dexpect[k] = w * dExp_y[k];
+
+                f1w = f1 * w;
+                f2ww = f2 * w * w;
+                for (k = 0; k < q*q; k++)
+                    dcovar[k] = f1w * dCov_y[k] - f2ww * dCov_y[k];
+            } else {
+                continue;
+            }
+        
+            /* the absolute standardized test statistic, to be maximized */
+            /* statistics[j] = C_maxabsTestStatistic(dlinstat, 
+                   dexpect, dcovar, q, tol); */
+
+            /* much faster but uses maxabs always*/
+            statistics[j] = 0.0;
+            for (k = 0; k < q; k++) {
+                if (dcovar[k * q + k] <= tol) continue;
+                tmp = fabs(dlinstat[k] - dexpect[k]) / sqrt(dcovar[k * q + k]);
+                if (statistics[j] < tmp) statistics[j] = tmp;
+            }
+
+        }
+    
+        /* search for the maximum and the best separating cutpoint */
+        /* <FIXME> the result might differ between 32 and 64bit systems 
+                   because of rounding errors in 'statistics' */
+        maxstat[0] = 0.0;        
+        for (i = 0; i < n; i++) {
+            if (statistics[i] > maxstat[0]) {
+                maxstat[0] = statistics[i];
+                cutpoint[0] = x[i];
+            }
+        }
+        /* </FIXME> */
+    }
+}
+
+
+/**
+    R-interface to C_split (does not handle ordered y's)
+    *\param x values of the transformation
+    *\param y values of the influence function
+    *\param weights case weights
+    *\param orderx the ordering of the transformations, i.e. R> order(x)
+    *\param linexpcov2sample an (uninitialized) object of class 
+                             `LinStatExpectCovar' with p = 1
+    *\param expcovinf an initialized object of class `ExpectCovarInfluence'
+    *\param splitctrl an object of class `SplitControl'
+*/
+
+SEXP R_split(SEXP x, SEXP y, SEXP weights, SEXP orderx, SEXP linexpcov2sample, 
+             SEXP expcovinf, SEXP splitctrl) {
+             
+    SEXP ans, cutpoint, maxstat, statistics;
+    
+    PROTECT(ans = allocVector(VECSXP, 3));
+    SET_VECTOR_ELT(ans, 0, cutpoint = allocVector(REALSXP, 1));
+    SET_VECTOR_ELT(ans, 1, maxstat = allocVector(REALSXP, 1));
+    SET_VECTOR_ELT(ans, 2, statistics = allocVector(REALSXP, nrow(x)));
+    
+    C_split(REAL(x), ncol(x), REAL(y), ncol(y), REAL(weights), nrow(x),
+            INTEGER(orderx), splitctrl, linexpcov2sample, expcovinf,
+            REAL(cutpoint), REAL(maxstat), REAL(statistics));
+    UNPROTECT(1);
+    return(ans);
+}
+
+
+/**
+    Search for a cutpoint in a unordered factor x maximizing a two-sample
+    statistic w.r.t. (the influence function of ) the response variable y.
+    *\param codingx the coding of x, i.e. as.numeric(x)
+    *\param p dimension of the transformation
+    *\param y values of the influence function
+    *\param q dimension of the influence function
+    *\param weights case weights
+    *\param n number of observations
+    *\param codingx the coding of x, i.e. as.numeric(x)
+    *\param standstat the vector of the standardized statistics for x, y, 
+                      weights 
+    *\param splitctrl an object of class `SplitControl'
+    *\param linexpcov2sample an (uninitialized) object of class 
+                             `LinStatExpectCovar' with p = 1
+    *\param expcovinf an initialized object of class `ExpectCovarInfluence'
+    *\param cutpoint return value; pointer to a double for the cutpoint in x
+    *\param levelset return value; pointer to a p-dim 0/1 integer
+    *\param maxstat return value; pointer to a double for the 
+                    maximal test statistic
+    *\param statistics return value; pointer to a n-dim double for 
+                       the statistics
+*/
+
+void C_splitcategorical(const int *codingx, int p,
+                        const double *y, int q,
+                        const double *weights, int n,
+                        double *standstat,
+                        SEXP splitctrl, SEXP linexpcov2sample, 
+                        SEXP expcovinf, double *cutpoint, int *levelset, 
+                        double *maxstat, double *statistics) {
+
+    double *tmpx, *tmptmpx, tmp = 0.0;
+    int *irank, *ordertmpx, i, j, k, l, jp, chk;
+
+    /* allocate memory */
+    tmpx = Calloc(n, double);
+    ordertmpx = Calloc(n, int);
+    irank = Calloc(p, int);
+    tmptmpx = Calloc(n, double);
+
+    /* for all response variables (aka: dummy variables) */
+    for (j = 0; j < q; j++) {
+    
+        jp = j * p;
+
+        /* determine the ranking of the kth level among 
+           the standardized statistic: This induced an ordering of the 
+           observations */
+        for (k = 0; k < p; k++) {
+            irank[k] = 1;
+            for (l = 0; l < p; l++)
+                if (standstat[jp + l] < standstat[jp + k]) irank[k]++;
+        }
+        
+        /* a temporary response variable: the rank of the level */
+        for (i = 0; i < n; i++) {
+            /* <FIXME> do we have to adjust weights for missing values here??? */
+            if (weights[i] == 0.0) {
+                tmpx[i] = 0.0;
+            } else {
+                tmpx[i] = (double) irank[codingx[i] - 1];
+            }
+            /* </FIXME> */
+            tmptmpx[i] = tmpx[i];
+            ordertmpx[i] = i + 1;
+        }
+        
+        /* order(dtmpx) */
+        rsort_with_index(tmptmpx, ordertmpx, n);
+
+        /* search for a cutpoint (now we do have an ordering) */
+        C_split(tmpx, 1, y, q, weights, n, ordertmpx,
+                splitctrl, linexpcov2sample,
+                expcovinf, cutpoint, maxstat, statistics);
+
+        /* if we have seen an improvement: save this segmentation 
+           note: there may be splits with equal goodness */
+        chk = 0;
+        if (maxstat[0] > tmp) {
+            for (k = 0; k < p; k++) {
+                if (irank[k] > cutpoint[0]) {
+                    levelset[k] = 1;
+                    chk += 1;
+                } else {
+                    levelset[k] = 0;
+                }
+            }
+            tmp = maxstat[0];
+        }
+        /* <FIXME> make sure that at least one level goes left,
+                   C_split may end up with cutpoint > max(irank), why?
+           </FIXME>
+        */
+        /* hm, why did I added 
+        if (chk == 0) tmp = 0.0; 
+        ??? */
+    }
+    maxstat[0] = tmp;
+
+    /* free memory */
+    Free(tmpx); Free(ordertmpx); Free(irank); Free(tmptmpx);
+}
+
+
+/**
+    R-interface to C_splitcategorical (does not handle ordered y's)
+    *\param x the values of the x-transformation
+    *\param codingx the coding of x, i.e. as.numeric(x)
+    *\param y values of the influence function
+    *\param weights case weights
+    *\param linexpcov2sample an (uninitialized) object of class 
+                             `LinStatExpectCovar' with p = 1
+    *\param linexpcov an initialized object of class `LinStatExpectCovar'
+    *\param expcovinf an initialized object of class `ExpectCovarInfluence'
+    *\param splitctrl an object of class `SplitControl'
+*/
+
+SEXP R_splitcategorical(SEXP x, SEXP codingx, SEXP y, SEXP weights, 
+                        SEXP linexpcov2sample, SEXP linexpcov, 
+                        SEXP expcovinf, SEXP splitctrl) {
+             
+    SEXP ans, cutpoint, maxstat, statistics, levelset;
+    double *standstat;
+
+    C_LinStatExpCov(REAL(x), ncol(x), REAL(y), ncol(y), REAL(weights), nrow(x),
+                    1, GET_SLOT(linexpcov, PL2_expcovinfSym), linexpcov);
+
+    standstat = Calloc(get_dimension(linexpcov), double);
+    C_standardize(REAL(GET_SLOT(linexpcov, PL2_linearstatisticSym)),
+                  REAL(GET_SLOT(linexpcov, PL2_expectationSym)),
+                  REAL(GET_SLOT(linexpcov, PL2_covarianceSym)),
+                  get_dimension(linexpcov), get_tol(splitctrl), standstat);
+
+    PROTECT(ans = allocVector(VECSXP, 4));
+    SET_VECTOR_ELT(ans, 0, cutpoint = allocVector(REALSXP, 1));
+    SET_VECTOR_ELT(ans, 1, maxstat = allocVector(REALSXP, 1));
+    SET_VECTOR_ELT(ans, 2, statistics = allocVector(REALSXP, nrow(x)));
+    SET_VECTOR_ELT(ans, 3, levelset = allocVector(INTSXP, ncol(x)));
+    
+    C_splitcategorical(INTEGER(codingx), ncol(x), REAL(y), ncol(y), REAL(weights), 
+                       nrow(x), standstat, 
+                       splitctrl, linexpcov2sample, expcovinf, 
+                       REAL(cutpoint), INTEGER(levelset), REAL(maxstat), 
+                       REAL(statistics));
+
+    UNPROTECT(1);
+    Free(standstat);
+    return(ans);
+}