/**

    Node computations

    *\file Node.c

    *\author $Author$

    *\date $Date$

*/

#include "party.h"

/**

    Compute prediction of a node

    *\param y the response variable (raw numeric values or dummy encoded factor)

    *\param n number of observations

    *\param q number of columns of y

    *\param weights case weights

    *\param sweights sum of case weights

    *\param ans return value; the q-dimensional predictions

*/

void C_prediction(const double *y, int n, int q, const double *weights, 

                  const double sweights, double *ans) {

    int i, j, jn;

    for (j = 0; j < q; j++) {

        ans[j] = 0.0;

        jn = j * n;

        for (i = 0; i < n; i++) 

            ans[j] += weights[i] * y[jn + i];

        ans[j] = ans[j] / sweights;

    }

}

void mask_pvalue(double *pvalue, int *variables_to_ignore, int ninputs) {

  int i;

  if(variables_to_ignore == NULL) {

    return;

  }

  for(i=0;i<ninputs;i++) {

    if(variables_to_ignore[i]) {

        pvalue[i] = R_NegInf;

    }

  }

}

/**

    The main function for all node computations

    *\param node an initialized node (an S3 object!)

    *\param learnsample an object of class `LearningSample'

    *\param weights case weights

    *\param fitmem an object of class `TreeFitMemory'

    *\param controls an object of class `TreeControl'

    *\param TERMINAL logical indicating if this node will

                     be a terminal node

    *\param depth an integer giving the depth of the current node

*/

void C_Node(SEXP node, SEXP learnsample, SEXP weights, 

            SEXP fitmem, SEXP controls, int TERMINAL, int depth, int *variables_to_ignore) {

    int nobs, ninputs, jselect, q, j, k, i;

    double mincriterion, sweights, *dprediction;

    double *teststat, *pvalue, smax, cutpoint = 0.0, maxstat = 0.0;

    double *standstat, *splitstat;

    SEXP responses, inputs, x, expcovinf, linexpcov;

    SEXP varctrl, splitctrl, gtctrl, tgctrl, split, testy, predy;

    double *dxtransf, *thisweights;

    int *itable;

    nobs = get_nobs(learnsample);

    ninputs = get_ninputs(learnsample);

    varctrl = get_varctrl(controls);

    splitctrl = get_splitctrl(controls);

    gtctrl = get_gtctrl(controls);

    tgctrl = get_tgctrl(controls);

    mincriterion = get_mincriterion(gtctrl);

    responses = GET_SLOT(learnsample, PL2_responsesSym);

    inputs = GET_SLOT(learnsample, PL2_inputsSym);

    testy = get_test_trafo(responses);

    predy = get_predict_trafo(responses);

    q = ncol(testy);

    /* <FIXME> we compute C_GlobalTest even for TERMINAL nodes! </FIXME> */

    /* compute the test statistics and the node criteria for each input */        

    C_GlobalTest(learnsample, weights, fitmem, varctrl,

                 gtctrl, get_minsplit(splitctrl), 

                 REAL(S3get_teststat(node)), REAL(S3get_criterion(node)), depth);

    /* sum of weights: C_GlobalTest did nothing if sweights < mincriterion */

    sweights = REAL(GET_SLOT(GET_SLOT(fitmem, PL2_expcovinfSym), 

                             PL2_sumweightsSym))[0];

    REAL(VECTOR_ELT(node, S3_SUMWEIGHTS))[0] = sweights;

    /* compute the prediction of this node */

    dprediction = REAL(S3get_prediction(node));

    /* <FIXME> feed raw numeric values OR dummy encoded factors as y 

       Problem: what happens for survival times ? */

    C_prediction(REAL(predy), nobs, ncol(predy), REAL(weights), 

                     sweights, dprediction);

    /* </FIXME> */

    teststat = REAL(S3get_teststat(node));

    pvalue = REAL(S3get_criterion(node));

    mask_pvalue(pvalue, variables_to_ignore, ninputs);

    /* try the two out of ninputs best inputs variables */

    /* <FIXME> be more flexible and add a parameter controlling

               the number of inputs tried </FIXME> */

    for (j = 0; j < 2; j++) {

        smax = C_max(pvalue, ninputs);

        REAL(S3get_maxcriterion(node))[0] = smax;

        /* if the global null hypothesis was rejected */

        if (smax > mincriterion && !TERMINAL) {

            /* the input variable with largest association to the response */

            jselect = C_whichmax(pvalue, teststat, ninputs) + 1;

            /* get the raw numeric values or the codings of a factor */

            x = get_variable(inputs, jselect);

            if (has_missings(inputs, jselect)) {

                expcovinf = GET_SLOT(get_varmemory(fitmem, jselect), 

                                    PL2_expcovinfSym);

                thisweights = C_tempweights(jselect, weights, fitmem, inputs);

            } else {

                expcovinf = GET_SLOT(fitmem, PL2_expcovinfSym);

                thisweights = REAL(weights);

            }

            /* <FIXME> handle ordered factors separatly??? </FIXME> */

            if (!is_nominal(inputs, jselect)) {

                /* search for a split in a ordered variable x */

                split = S3get_primarysplit(node);

                /* check if the n-vector of splitstatistics 

                   should be returned for each primary split */

                if (get_savesplitstats(tgctrl)) {

                    C_init_orderedsplit(split, nobs);

                    splitstat = REAL(S3get_splitstatistics(split));

                } else {

                    C_init_orderedsplit(split, 0);

                    splitstat = REAL(get_splitstatistics(fitmem));

                }

                C_split(REAL(x), 1, REAL(testy), q, thisweights, nobs,

                        INTEGER(get_ordering(inputs, jselect)), splitctrl, 

                        GET_SLOT(fitmem, PL2_linexpcov2sampleSym),

                        expcovinf, REAL(S3get_splitpoint(split)), &maxstat,

                        splitstat);

                S3set_variableID(split, jselect);

             } else {

                 /* search of a set of levels (split) in a numeric variable x */

                 split = S3get_primarysplit(node);

                /* check if the n-vector of splitstatistics 

                   should be returned for each primary split */

                if (get_savesplitstats(tgctrl)) {

                    C_init_nominalsplit(split, 

                        LENGTH(get_levels(inputs, jselect)), 

                        nobs);

                    splitstat = REAL(S3get_splitstatistics(split));

                } else {

                    C_init_nominalsplit(split, 

                        LENGTH(get_levels(inputs, jselect)), 

                        0);

                    splitstat = REAL(get_splitstatistics(fitmem));

                }

                 linexpcov = get_varmemory(fitmem, jselect);

                 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);

                 C_splitcategorical(INTEGER(x), 

                                    LENGTH(get_levels(inputs, jselect)), 

                                    REAL(testy), q, thisweights, 

                                    nobs, standstat, splitctrl, 

                                    GET_SLOT(fitmem, PL2_linexpcov2sampleSym),

                                    expcovinf, &cutpoint, 

                                    INTEGER(S3get_splitpoint(split)),

                                    &maxstat, splitstat);

                 /* compute which levels of a factor are available in this node 

                    (for printing) later on. A real `table' for this node would

                    induce too much overhead here. Maybe later. */

                 itable = INTEGER(S3get_table(split));

                 dxtransf = REAL(get_transformation(inputs, jselect));

                 for (k = 0; k < LENGTH(get_levels(inputs, jselect)); k++) {

                     itable[k] = 0;

                     for (i = 0; i < nobs; i++) {

                         if (dxtransf[k * nobs + i] * thisweights[i] > 0) {

                             itable[k] = 1;

                             continue;

                         }

                     }

                 }

                 Free(standstat);

            }

            if (maxstat == 0) {

                if (j == 1) {          

                    S3set_nodeterminal(node);

                } else {

                    /* do not look at jselect in next iteration */

                    pvalue[jselect - 1] = R_NegInf;

                }

            } else {

                S3set_variableID(split, jselect);

                break;

            }

        } else {

            S3set_nodeterminal(node);

            break;

        }

    }

}       

/**

    R-interface to C_Node

    *\param learnsample an object of class `LearningSample'

    *\param weights case weights

    *\param fitmem an object of class `TreeFitMemory'

    *\param controls an object of class `TreeControl'

*/

SEXP R_Node(SEXP learnsample, SEXP weights, SEXP fitmem, SEXP controls) {

     SEXP ans;

     PROTECT(ans = allocVector(VECSXP, NODE_LENGTH));

     C_init_node(ans, get_nobs(learnsample), get_ninputs(learnsample), 

                 get_maxsurrogate(get_splitctrl(controls)),

                 ncol(get_predict_trafo(GET_SLOT(learnsample, PL2_responsesSym))));

     C_Node(ans, learnsample, weights, fitmem, controls, 0, 1, NULL);

     UNPROTECT(1);

     return(ans);

}