/**
    Conditional Distributions
    *\file Distributions.c
    *\author $Author$
    *\date $Date$
*/
                
#include "party.h"


/**
    Conditional asymptotic P-value of a quadratic form\n
    *\param tstat test statistic
    *\param df degree of freedom
*/
        
double C_quadformConditionalPvalue(const double tstat, const double df) {
    return(pchisq(tstat, df, 0, 0));
}


/**
    R-interface to C_quadformConditionalPvalue\n
    *\param tstat test statitstic
    *\param df degree of freedom
*/

SEXP R_quadformConditionalPvalue(SEXP tstat, SEXP df) {

    SEXP ans;
    
    PROTECT(ans = allocVector(REALSXP, 1));
    REAL(ans)[0] = C_quadformConditionalPvalue(REAL(tstat)[0], REAL(df)[0]);
    UNPROTECT(1);
    return(ans);
}


/**
    Conditional asymptotic P-value of a maxabs-type test statistic\n
    Basically the functionality from package `mvtnorm' \n
    *\param tstat test statitstic
    *\param Sigma covariance matrix
    *\param pq nrow(Sigma)
    *\param maxpts number of Monte-Carlo steps
    *\param releps relative error
    *\param abseps absolute error
    *\param tol tolerance
*/

double C_maxabsConditionalPvalue(const double tstat, const double *Sigma, 
    const int pq, int *maxpts, double *releps, double *abseps, double *tol) {

    int *n, *nu, *inform, i, j, *infin, sub, *index, nonzero, iz, jz;
    double *lower, *upper, *delta, *corr, *sd, *myerror,
           *prob, ans;

    /* univariate problem */
    if (pq == 1) 
        return(2*pnorm(fabs(tstat)*-1.0, 0.0, 1.0, 1, 0)); /* return P-value */
    
    n = Calloc(1, int);
    nu = Calloc(1, int);
    myerror = Calloc(1, double);
    prob = Calloc(1, double);
    nu[0] = 0;
    inform = Calloc(1, int);
    n[0] = pq;
        
    if (n[0] == 2)  
         corr = Calloc(1, double);
    else 
         corr = Calloc(n[0] + ((n[0] - 2) * (n[0] - 1))/2, double);
    
    sd = Calloc(n[0], double);
    lower = Calloc(n[0], double);
    upper = Calloc(n[0], double);
    infin = Calloc(n[0], int);
    delta = Calloc(n[0], double);
    index = Calloc(n[0], int);

    /* determine elements with non-zero variance */ 

    nonzero = 0;
    for (i = 0; i < n[0]; i++) {
        if (Sigma[i*n[0] + i] > tol[0]) {
            index[nonzero] = i;
            nonzero++;
        }
    }

    /* mvtdst assumes the unique elements of the triangular 
       covariance matrix to be passes as argument CORREL 
    */

    for (iz = 0; iz < nonzero; iz++) {

        /* handle elements with non-zero variance only */
        i = index[iz];

        /* standard deviations */
        sd[i] = sqrt(Sigma[i*n[0] + i]);
                
        /* always look at the two-sided problem */           
        lower[iz] = fabs(tstat) * -1.0;
        upper[iz] = fabs(tstat);
        infin[iz] = 2;
        delta[iz] = 0.0;
        
        /* set up vector of correlations, i.e., the upper 
           triangular part of the covariance matrix) */
        for (jz = 0; jz < iz; jz++) {
            j = index[jz]; 
            sub = (int) (jz + 1) + (double) ((iz - 1) * iz) / 2 - 1;
            if (sd[i] == 0.0 || sd[j] == 0.0) 
                corr[sub] = 0.0; 
            else 
                corr[sub] = Sigma[i*n[0] + j] / (sd[i] * sd[j]);
        }
    }
    n[0] = nonzero;
        
    /* call FORTRAN subroutine */
    F77_CALL(mvtdst)(n, nu, lower, upper, infin, corr, delta, 
                     maxpts, abseps, releps, myerror, prob, inform);
                         
    /* inform == 0 means: everything is OK */
    switch (inform[0]) {
        case 0: break;
        case 1: warning("cmvnorm: completion with ERROR > EPS"); break;
        case 2: warning("cmvnorm: N > 1000 or N < 1"); 
                prob[0] = 0.0; 
                break;
        case 3: warning("cmvnorm: correlation matrix not positive semi-definite"); 
                prob[0] = 0.0; 
                break;
        default: warning("cmvnorm: unknown problem in MVTDST");
                 prob[0] = 0.0;
    }
    ans = prob[0];
    Free(corr); Free(sd); Free(lower); Free(upper); 
    Free(infin); Free(delta); Free(myerror); Free(prob);
    Free(n); Free(nu); Free(inform); 
    return(1 - ans);  /* return P-value */
}


/**
   R-interface to C_maxabsConditionalPvalue \n
   *\param tstat test statitstic
   *\param Sigma covariance matrix
   *\param maxpts number of Monte-Carlo steps
   *\param releps relative error
   *\param abseps absolute error
   *\param tol tolerance
*/

SEXP R_maxabsConditionalPvalue(SEXP tstat, SEXP Sigma, SEXP maxpts, 
                               SEXP releps, SEXP abseps, SEXP tol) {
           
    SEXP ans;                    
    int pq;
    
    pq = nrow(Sigma);
   
    PROTECT(ans = allocVector(REALSXP, 1));
    REAL(ans)[0] = C_maxabsConditionalPvalue(REAL(tstat)[0], REAL(Sigma), pq, 
        INTEGER(maxpts), REAL(releps), REAL(abseps), REAL(tol));
    UNPROTECT(1);
    return(ans);
}


/**
    Monte-Carlo approximation to the conditional pvalues 
    *\param criterion vector of node criteria for each input 
    *\param learnsample an object of class `LearningSample'
    *\param weights case weights
    *\param fitmem an object of class `TreeFitMemory'
    *\param varctrl an object of class `VariableControl'
    *\param gtctrl an object of class `GlobalTestControl'
    *\param ans_pvalues return values; vector of adjusted pvalues 
*/

void C_MonteCarlo(double *criterion, SEXP learnsample, SEXP weights, 
                  SEXP fitmem, SEXP varctrl, SEXP gtctrl, double *ans_pvalues) {

    int ninputs, nobs, j, i, k;
    SEXP responses, inputs, y, x, xmem, expcovinf;
    double sweights, *stats, tmp = 0.0, smax, *dweights; 
    int m, *counts, b, B, *dummy, *permindex, *index, *permute;
    
    ninputs = get_ninputs(learnsample);
    nobs = get_nobs(learnsample);
    responses = GET_SLOT(learnsample, PL2_responsesSym);
    inputs = GET_SLOT(learnsample, PL2_inputsSym);
    dweights = REAL(weights);
    
    /* number of Monte-Carlo replications */
    B = get_nresample(gtctrl);
    
    /* y = get_transformation(responses, 1); */
    y = get_test_trafo(responses);
    
    expcovinf = GET_SLOT(fitmem, PL2_expcovinfSym);

    sweights = REAL(GET_SLOT(expcovinf, PL2_sumweightsSym))[0];
    m = (int) sweights;
    
    stats = Calloc(ninputs, double);
    counts = Calloc(ninputs, int);
    
    dummy = Calloc(m, int);
    permute = Calloc(m, int);
    index = Calloc(m, int);
    permindex = Calloc(m, int);
                
    /* expand weights, see appendix of 
       `Unbiased Recursive Partitioning: A Conditional Inference Framework' */
    j = 0;
    for (i = 0; i < nobs; i++) {
        for (k = 0; k < dweights[i]; k++) {
            index[j] = i;
            j++;
        }
    }

    for (b = 0; b < B; b++) {

        /* generate a admissible permutation */
        C_SampleNoReplace(dummy, m, m, permute);
        for (k = 0; k < m; k++) permindex[k] = index[permute[k]];

        /* for all input variables */
        for (j = 1; j <= ninputs; j++) {
            x = get_transformation(inputs, j);

            /* compute test statistic or pvalue for the permuted data */
            xmem = get_varmemory(fitmem, j);
            if (!has_missings(inputs, j)) {
                C_PermutedLinearStatistic(REAL(x), ncol(x), REAL(y), ncol(y), 
                    nobs, m, index, permindex, 
                    REAL(GET_SLOT(xmem, PL2_linearstatisticSym)));
            } else {
                error("cannot resample with missing values");
            }
            
            /* compute the criterion, i.e. something to be MAXIMISED */
            C_TeststatCriterion(xmem, varctrl, &tmp, &stats[j - 1]);
        }
        
        /* the maximum of the permuted test statistics / 1 - pvalues */
        smax = C_max(stats, ninputs);

        /* count the number of permuted > observed */
        for (j = 0; j < ninputs; j++) {
            if (smax > criterion[j]) counts[j]++;
        }
    }
    
    /* return adjusted pvalues */
    for (j = 0; j < ninputs; j++)
        ans_pvalues[j] = (double) counts[j] / B;
                
    /* <FIXME> we try to assess the linear statistics later on 
               (in C_Node, for categorical variables) 
               but have used this memory for resampling here */

    for (j = 1; j <= ninputs; j++) {
        x = get_transformation(inputs, j);
        /* re-compute linear statistics for unpermuted data */
        xmem = get_varmemory(fitmem, j);
        C_LinearStatistic(REAL(x), ncol(x), REAL(y), ncol(y), 
                      dweights, nobs, 
                      REAL(GET_SLOT(xmem, PL2_linearstatisticSym)));
    }
    /* </FIXME> */
    
    Free(stats); Free(counts); Free(dummy); Free(permute); 
    Free(index); Free(permindex);
}


/**
    R-interface to C_MonteCarlo \n
    *\param criterion vector of node criteria for each input
    *\param learnsample an object of class `LearningSample'
    *\param weights case weights
    *\param fitmem an object of class `TreeFitMemory'
    *\param varctrl an object of class `VariableControl'
    *\param gtctrl an object of class `GlobalTestControl'
*/

SEXP R_MonteCarlo(SEXP criterion, SEXP learnsample, SEXP weights, 
                  SEXP fitmem, SEXP varctrl, SEXP gtctrl) {
                  
     SEXP ans;
     
     GetRNGstate();
     
     PROTECT(ans = allocVector(REALSXP, get_ninputs(learnsample)));
     C_MonteCarlo(REAL(criterion), learnsample, weights, fitmem, varctrl, 
                  gtctrl, REAL(ans));
                  
     PutRNGstate();
                  
     UNPROTECT(1);
     return(ans);
}