/**
    Test statistics for conditional inference
    *\file TestStatistic.c
    *\author $Author$
    *\date $Date$
*/
                
#include "party.h"


/**
    Standardizes a statistic t of length pq with mean mu and covariance Sigma
    for variances > tol \n
    *\param t the vector of statistics
    *\param mu expectations
    *\param Sigma covariance matrix
    *\param pq dimension of t
    *\param tol tolerance for variances
    *\param ans return value; a pointer to a REALSXP-vector of length pq
*/

void C_standardize(const double *t, const double *mu, const double *Sigma, 
                   int pq, double tol, double *ans) {
                   
    int i;
    double sd;
    
    for (i = 0; i < pq; i++) {
        sd = Sigma[i*pq + i]; 
        if (sd > tol)
            ans[i] = (t[i] - mu[i])/sqrt(sd);
        else
            ans[i] = 0.0;
    }
}


/**
    Absolute values of standardized statistics
    *\param t the vector of statistics
    *\param mu expectations
    *\param Sigma covariance matrix
    *\param pq dimension of t
    *\param tol tolerance for variances
    *\param ans return value; a pointer to a REALSXP-vector of length pq
*/

void C_absstandardize(const double *t, const double *mu, const double *Sigma,
                      int pq, double tol, double *ans) {
                      
    C_standardize(t, mu, Sigma, pq, tol, ans);
    C_abs(ans, pq);
}


/**
    Maximum absolute values of standardized statistics
    *\param t the vector of statistics
    *\param mu expectations
    *\param Sigma covariance matrix
    *\param pq dimension of t
    *\param tol tolerance for variances
*/

double C_maxabsTestStatistic(const double *t, const double *mu, const double *Sigma,
                             int pq, double tol) {
                           
     double *mem, ans;
     
     mem = Calloc(pq, double);
     C_absstandardize(t, mu, Sigma, pq, tol, mem);
     ans = C_max(mem, pq);
     Free(mem);
     return(ans);
}


/**
    R-interface to C_maxabsTestStatistic
    *\param t the vector of statistics
    *\param mu expectations
    *\param Sigma covariance matrix
    *\param tol tolerance for variances
*/

SEXP R_maxabsTestStatistic(SEXP t, SEXP mu, SEXP Sigma, SEXP tol) {

     SEXP ans;
     int pq;
     
     pq = LENGTH(t);
     
     PROTECT(ans = allocVector(REALSXP, 1));
     REAL(ans)[0] = C_maxabsTestStatistic(REAL(t), REAL(mu), REAL(Sigma), pq, 
                                          REAL(tol)[0]);
     UNPROTECT(1);
     return(ans);
}


/**
    Quadratic form t(t - mu) SigmaPlus (t - mu) \n
    *\param t the vector of statistics
    *\param mu expectations
    *\param SigmaPlus Moore-Penrose inverse
    *\param pq dimension of t
*/
                                
double C_quadformTestStatistic(const double *t, const double *mu, 
                               const double *SigmaPlus, int pq) {

    int i, j;
    double quadform = 0.0, *tmmu, *tmmuSigmaPlus;

    tmmu = Calloc(pq, double);
    for (i = 0; i < pq; i++)
        tmmu[i] = t[i] - mu[i];
    
    tmmuSigmaPlus = Calloc(pq, double);
    for (i = 0; i < pq; i++)  {
        tmmuSigmaPlus[i] = 0.0;
        for (j = 0; j < pq; j++)
            tmmuSigmaPlus[i] += tmmu[j] * SigmaPlus[i * pq + j];
        quadform += tmmuSigmaPlus[i] * tmmu[i];
    }

    Free(tmmu); Free(tmmuSigmaPlus);
    return(quadform);
}


/**
    R-interface to C_quadformTestStatistic \n
    *\param t the vector of statistics
    *\param mu expectations
    *\param SigmaPlus Moore-Penrose inverse
*/
    
SEXP R_quadformTestStatistic(SEXP t, SEXP mu, SEXP SigmaPlus) {

    SEXP ans;
    int pq;
    
    pq = LENGTH(t);
    PROTECT(ans = allocVector(REALSXP, 1));
    REAL(ans)[0] = C_quadformTestStatistic(REAL(t), 
                                           REAL(mu), REAL(SigmaPlus), pq);
    UNPROTECT(1);
    return(ans);
}