#include <math.h>
#include <stdio.h>
#include <stdlib.h>
/* Needed for building mex file:                                           */
#include <mex.h>
#include <matrix.h>
 

/*  Copyright (C) 2009 Cesare Magri                                        */

/*
 * -------
 * LICENSE
 * -------
 * This software is distributed free under the condition that:
 *
 * 1. it shall not be incorporated in software that is subsequently sold;
 *
 * 2. the authorship of the software shall be acknowledged and the following
 *    article shall be properly cited in any publication that uses results
 *    generated by the software:

 *      Magri C, Whittingstall K, Singh V, Logothetis NK, Panzeri S: A
 *      toolbox for the fast information analysis of multiple-site LFP, EEG
 *      and spike train recordings. BMC Neuroscience 2009 10(1):81;
 * 3.  this notice shall remain in place in each source file.
 */


/* Global variables:                                                       */
static bool initialized = false;
static mxArray *ClogC_mxArray = NULL;


/* ========================================================================*/
/* VECTORIAL KRONECKER'S PRODUCT                                           */
/* ========================================================================*/
void vecKronProd(
                 double  *xPtr,      /* input: pointer to array x          */
                 mwSize   xLen,      /* input: length of x                 */
                 double  *yPtr,      /* input: pointer to array y          */
                 mwSize   yLen,      /* input: length of y                 */
                 double **K_ptr_ptr, /* output                             */
                 mwSize*  KLen_ptr   /* output                             */
                )
{   
    double *K_ptr;
    mwIndex i, j, indx;
    
    /* Where to store the final array:                                     */
    K_ptr = *K_ptr_ptr;
    
    /* Length of final array:                                              */
    *KLen_ptr = xLen*yLen;
    
    K_ptr = mxCalloc(*KLen_ptr, sizeof(double));
    
    /* In the main funciton we map bins into responses as follows:         */
    /*     r = R[0,t,s] + R[1,t,s]*Nb[0] + R[2,t,s]*Nb[0]*Nb[1] + ...      */
    /* this explains why we first loop over j and then over i in the       */
    /* following nested loop:                                              */
    indx = 0;
    for(j=0; j<yLen; j++)
        for(i=0; i<xLen; i++)
            K_ptr[indx++] = xPtr[i] * yPtr[j];
    
    /* Storing result in the right place:                                  */
    *K_ptr_ptr = K_ptr;
}



/* ========================================================================*/
/* SPLIT PRODUCTS                                                          */
/* ========================================================================*/
void splitProds(
                double *X_ptr, /* pointer to data matrix                   */
                mwSize M,      /* number of lines of X                     */
                mwSize N,      /* number of columns of X                   */
                mwSize* n_ptr, /* pointer to array which specifies how many*/
                               /* elements in each row of X are really     */
                               /* occupied by data                         */
                
                double** K_ptr_ptr,
                mwSize*  KLen_ptr
               )
{      
    double *K_ptr;
    
    /* If the number of columns of X is greater than two the fucntion calls*/
    /* itself recursively. The following variable are thus used to break   */
    /* the problem in two parts.                                           */
    double *K1_ptr, *K2_ptr;
    mwSize  K1Len, K2Len;
    mwSize  N1, N2;
    
    K_ptr = *K_ptr_ptr;
    
    /* If N=2 compute the product directly                                 */
    if(N==2) {
        vecKronProd(&X_ptr[0], n_ptr[0], &X_ptr[M], n_ptr[1], &K_ptr, KLen_ptr);

    /* Otherwise split X in two parts and work on each half independently  */
    } else {
        /* First haf ------------------------------------------------------*/
        N1 = N/2; /* note: floored int ratio */
        if(N1==1) {
            K1_ptr = X_ptr;
            K1Len = n_ptr[0];
        } else {
            splitProds(X_ptr, M, N1, n_ptr, &K1_ptr, &K1Len);
        }

        /* Second half ----------------------------------------------------*/
        N2 = N - N1;
        splitProds(&X_ptr[N1*M], M, N2, &n_ptr[N1], &K2_ptr, &K2Len);

        /* Final product --------------------------------------------------*/
        vecKronProd(K1_ptr, K1Len, K2_ptr, K2Len, &K_ptr, KLen_ptr);

        if(N1>1) /* Freeing K1 when N1=1 would erase part of X */
            mxFree(K1_ptr);

        mxFree(K2_ptr);

    }

    *K_ptr_ptr = K_ptr;
}



/* ========================================================================*/
/* CREATE ClogC                                                            */
/* ========================================================================*/
double *create_ClogC(mwSize Ns, mwSize Nt)
{
    double *ClogC;
    
    ClogC_mxArray = mxCreateDoubleMatrix(Ns*Nt+100, 1, mxREAL);
    mexMakeArrayPersistent(ClogC_mxArray);
    ClogC = mxGetPr(ClogC_mxArray);
    
    return ClogC;
}



/* ========================================================================*/
/* CLEANUP FUNCTION                                                        */
/* ========================================================================*/
void cleanup(void) {
    mexPrintf("MEX-file is terminating, destroying MEX arrays\n");
    mxDestroyArray(ClogC_mxArray);
}



/* ========================================================================*/
/* INCREMENT                                                               */
/* ========================================================================*/
double increment(double C, double *ClogC)
{
    double out;
    mwSize indx;
    
    indx = (mwIndex) C;

    if(ClogC[indx+1] == 0)
        ClogC[indx] = (C) * log(C);

    out = ClogC[indx] - ClogC[indx-1];
        
    return out;
}



/* ========================================================================*/
/* ENTROPY COMPUTATION (GATEWAY FUNCTION)                                  */
/* ========================================================================*/
void mexFunction(int nlhs, mxArray *plhs[],
                 int nrhs, const mxArray *prhs[])
{

    /* Input --------------------------------------------------------------*/
    
    double  *R;         /* response matrix                                 */
    double  *Nt;        /* number of trials per stimulus                   */
    double   biasCorrNum;

    bool     doHR;      /* compute H(R)       (flag)                       */
    bool     doHRS;     /* compute H(R|S)     (flag)                       */
    bool     doHlR;     /* compute H_lin(R)   (flag)                       */
    bool     doHlRS;    /* compute H_lin(R|S) (flag)                       */
    bool     doHiR;     /* compute H_ind(R)   (flag)                       */
    bool     doHiRS;    /* compute H_ind(R|S) (flag)                       */
    bool     doChiR;    /* compute Chi(R)     (flag)                       */
    bool     doHshR;    /* compute H_sh(R)    (flag)                       */
    bool     doHshRS;   /* compute H_sh(R|S)  (flag)                       */
    
    /* Internal variables -------------------------------------------------*/
    
    mxArray *PrToInputMxArray;  /* generic pointer to mxArray, used to read*/
                                /* several input                           */
    
    /* Other useful flags:                                                 */
    bool     doMap;     /* map multi-dim response into scalar (flag)       */
    bool     doCr;      /* compute C(r)       (flag)                       */
    bool     doPir;     /* compute P_ind(r)   (flag)                       */
    bool     doCirs;    /* compute C_ind(r|s) (flag)                       */
    bool     doCrcs;    /* compute C(r_c|s)   (flag)                       */
    bool     doSh;      /* do the shuffling   (flag)                       */
    
    /* Dimensions ---------------------------------------------------------*/
    mwSize   L;         /* number of cells                                 */
    mwSize   Ns;        /* number of stimuli                               */
    mwSize  *Nb;        /* number of bins per cell                         */
    mwSize   Nr;        /* number of possible responses                    */
    
    mwSize   maxNb;     /* max number of bins: maxNb=max(Nb)               */
    mwSize  *cumprodNb; /* cumulative product of Nb:                       */
                        /*     cumprodNb(n) = Nb(n)*Nb(n-1)*...*Nb(1)      */
    mwSize  *base;      /* base used for mapping each response from an     */
                        /* L-dimensional array to single number:           */
                        /*     base = cumprodNb ./ Nb(1)                   */
    
    /* Storing useful dimension-related quantities:                        */
    mwSize   maxNt;             /* max number of trials: maxNt=max(Nt)     */
    double   inv_Nt;            /* inverse of Nt: inv_Nt=1/Nt              */
    double   NtpowL;            /* Nt power L: NtpowL=Nt^L                 */
    double   inv_NtpowL;        /* inverse of NtpowL: inv_NtpowL=1/(Nt^L)  */
    double   inv_Ntpow_Lminus1; /* inverse of Nt power (L - 1):            */
                                /*   inv_Ntpow_Lminus1 = 1/(Nt^(L-1))      */
    
    double   totNt;     /* total numner of trials: totNt=sum(Nt)           */
    double   inv_totNt; /* inverse of totNt: inv_totNt=1/totNt             */
    
    /* Other useful quantities                                             */
    double   log_2 = log(2);

    /* Indexes                                                             */
    mwIndex  c;         /* cell index                                      */
    mwIndex  t;         /* trial index                                     */
    mwIndex  s;         /* stimulus index                                  */
    mwIndex  b;         /* bin index                                       */
    mwIndex  r;         /* response index                                  */

    mwIndex  i;         /* counter                                         */
    mwIndex  j;         /* counter                                         */
    mwIndex  indx;      /* used used to store index to elements of a matrix*/

    /* Variables used for computing cell-shuffled quantities:              */
    double  *Rsh;       /* shuffled response matrix                        */
    mwIndex  rsh;       /* shuffled response index                         */
    mwIndex  bsh;       /* shuffled bin index;                             */
    mwIndex  tsh;       /* shuffled trial index                            */
    mwIndex  indxSh;    /* used used to store index to an element of a     */
                        /* shuffled matrix                                 */

    /* Count and probability matrices:                                     */
    double  *Cr;        /* pointer to C(r) values                          */
    double  *Crs;       /* pointer to C(r|s) values                        */
    double  *Crc;       /* pointer to C(rc) values                         */
    double  *Crcs;      /* pointer to C(rc|s) values                       */
    double  *Pir;       /* pointer to C_ind(r) values                      */
    double   Cirs;      /* pointer to C_ind(rc|s) values                   */
    double  *Cshr;      /* pointer to C_sh(r) values                       */
    double  *Cshrs;     /* pointer to C_sh(r|s) values                     */
    
    /* C*log(C) function:                                                  */
    double  *ClogC;
    mwSize   ClogCLen;
    
    /* Variables needed for calls to splitProd                             */
    double  *X_ptr;
    mwSize   M;
    mwSize   N1;
    mwSize   N2;
    mwSize  *n_ptr;
    
    double  *K1_ptr;
    double  *K2_ptr;
    
    mwSize   K1Len;
    mwSize   K2Len;
    
    /* Bias Values                                                         */
    double  *biasPr;
    double   biasHlRS;
    double   biasHlR;
    
    /* Output -------------------------------------------------------------*/
    
    double  *HR;
    double  *HRS;
    double  *HlR;
    double  *HlRS;
    double  *HiR;
    double  *HiRS;
    double  *ChiR;
    double  *HshR;
    double  *HshRS;
    
    /* Reading input ------------------------------------------------------*/
    R  = mxGetPr(prhs[0]);
    
    PrToInputMxArray = mxGetField(prhs[1], 0, "Nt");
    Nt = mxGetPr(PrToInputMxArray);
    
    PrToInputMxArray = mxGetField(prhs[1], 0, "biasCorrNum");
    biasCorrNum = *mxGetPr(PrToInputMxArray);

    PrToInputMxArray = mxGetField(prhs[1], 0, "doHR");
    doHR = *mxGetLogicals(PrToInputMxArray);
    
    PrToInputMxArray = mxGetField(prhs[1], 0, "doHRS");
    doHRS = *mxGetLogicals(PrToInputMxArray);
    
    PrToInputMxArray = mxGetField(prhs[1], 0, "doHlR");
    doHlR = *mxGetLogicals(PrToInputMxArray);
    
    PrToInputMxArray = mxGetField(prhs[1], 0, "doHlRS");
    doHlRS = *mxGetLogicals(PrToInputMxArray);
    
    PrToInputMxArray = mxGetField(prhs[1], 0, "doHiR");
    doHiR = *mxGetLogicals(PrToInputMxArray);
    
    PrToInputMxArray = mxGetField(prhs[1], 0, "doHiRS");
    doHiRS = *mxGetLogicals(PrToInputMxArray);
    
    PrToInputMxArray = mxGetField(prhs[1], 0, "doChiR");
    doChiR = *mxGetLogicals(PrToInputMxArray);
    
    PrToInputMxArray = mxGetField(prhs[1], 0, "doHshR");
    doHshR = *mxGetLogicals(PrToInputMxArray);
    
    PrToInputMxArray = mxGetField(prhs[1], 0, "doHshRS");
    doHshRS = *mxGetLogicals(PrToInputMxArray);
    
    /* Other useful flags:                                                 */
    doMap   = doHR || doHRS || doChiR;
    doCr    = doHR || doChiR;
    doPir   = doHiR || doChiR;
    doCirs  = doPir || doHiRS;
    doCrcs  = doHlRS || doCirs;
    doSh    = doHshR || doHshRS;

    /* Computing dimensions -----------------------------------------------*/
    
    L     = mxGetDimensions(prhs[0])[0];
    maxNt = mxGetDimensions(prhs[0])[1];
    Ns    = mxGetDimensions(prhs[0])[2];
    
    Nb         = mxCalloc(L, sizeof(mwSize));
    base       = mxCalloc(L, sizeof(mwSize));
    cumprodNb  = mxCalloc(L, sizeof(mwSize));
    
    if(doSh)
        Rsh = mxCalloc(L * maxNt * Ns, sizeof(double));
    
    maxNb = 0;
    Nr    = 1;
    for(c=0; c<L; c++) {
            
        for(s=0; s<Ns; s++) {
            
            for(t=0; t<Nt[s]; t++) {

                indx = c + t*L + s*L*maxNt;

                /* Find max(R(c,:,:))                                      */
                if(R[indx] > Nb[c]) {
                    Nb[c] = R[indx];
                }

                /* Shuffling                                               */
                /* ---------                                               */
                /* The shuffling routine is based on the following simple  */
                /* algorithm for permuting an array (i.e., any ordering    */
                /* of the array is equally probable or any element has     */
                /* equal chance of being in any position):                 */
                /*                                                         */
                /*     For i=1 to n-1:                                     */
                /*     - Let j=r(i+1)                                      */
                /*     - Swap a[i] and a[j]                                */
                /*                                                         */
                /* where r[n] is a random number generated between 0 and   */
                /* n-1.                                                    */

                if(doSh) {
                    /* Generating random number between 0 and t            */
                    tsh = rand() % (t+1);
                    indxSh = c + tsh*L + s*L*maxNt;

                    /* Swapping R[c, t, s] and R[c, tsh, s]                */
                    Rsh[indx]   = Rsh[indxSh];
                    Rsh[indxSh] = R[indx];
                }
            }
        }
        
        /* The response is assumed to start from zero, the number of bins  */
        /* must thus be increased by one                                   */
        Nb[c]++;
    
        /* Finding max number of bins for cell c                           */
        if(Nb[c] > maxNb) {
            maxNb = Nb[c];
        }
        
        /* Building cumulative product                                     */
        if(c==0) {
            cumprodNb[c] = Nb[c];
            base[c]      = 1;
        } else {
            cumprodNb[c] = cumprodNb[c-1] * Nb[c];
            base[c]      = cumprodNb[c-1];
        }
        
        /* Number of possible responses                                    */
        Nr *= Nb[c];
    }
    
    
    /* Allocating memory --------------------------------------------------*/
    if(doCr)    Cr    = mxCalloc(Nr             , sizeof(double));
    if(doHRS)   Crs   = mxCalloc(Nr * Ns        , sizeof(double));
    if(doHlR)   Crc   = mxCalloc(maxNb * L      , sizeof(double));
    if(doCrcs)  Crcs  = mxCalloc(maxNb * L * Ns , sizeof(double));
    if(doPir)   Pir   = mxCalloc(cumprodNb[L-1] , sizeof(double));
    if(doHshR)  Cshr  = mxCalloc(Nr             , sizeof(double));
    if(doHshRS) Cshrs = mxCalloc(Nr * Ns        , sizeof(double));
   

    /* HR                                                                  */
    plhs[0] = mxCreateDoubleScalar(0);
    HR = mxGetPr(plhs[0]);
    
    /* HRS                                                                 */
    if(doHRS) plhs[1] = mxCreateDoubleMatrix(Ns,1,mxREAL);
    else      plhs[1] = mxCreateDoubleScalar(0);
    HRS = mxGetPr(plhs[1]);

    /* HlR                                                                 */
    plhs[2] = mxCreateDoubleScalar(0);
    HlR = mxGetPr(plhs[2]);
    
    /* HlRS                                                                */
    if(doHlRS) plhs[3] = mxCreateDoubleMatrix(Ns,1,mxREAL);    
    else       plhs[3] = mxCreateDoubleScalar(0);
    HlRS = mxGetPr(plhs[3]);
    
    /* HiR                                                                 */
    plhs[4] = mxCreateDoubleScalar(0);
    HiR = mxGetPr(plhs[4]);
    
    /* HiRS                                                                */
    if(doHiRS) plhs[5] = mxCreateDoubleMatrix(Ns,1,mxREAL);
    else       plhs[5] = mxCreateDoubleScalar(0);
    HiRS = mxGetPr(plhs[5]);

    /* ChiR                                                                */
    plhs[6] = mxCreateDoubleScalar(0);
    ChiR = mxGetPr(plhs[6]);
    
    /* HshR                                                                */
    plhs[7] = mxCreateDoubleScalar(0);
    HshR = mxGetPr(plhs[7]);
    
    /* HshRS                                                               */
    if(doHshRS) plhs[8] = mxCreateDoubleMatrix(Ns,1,mxREAL);
    else        plhs[8] = mxCreateDoubleScalar(0);
    HshRS = mxGetPr(plhs[8]);
    
    biasPr = mxCalloc(1, sizeof(double));
    
    /* Initializing ClogC -------------------------------------------------*/
    
    /* If ClogC has not been already initialized:                          */
    if (!initialized) {
        /* Create persistent array and register its cleanup:               */
        ClogC = create_ClogC(Ns, maxNt);
        mexAtExit(cleanup);
        initialized = true;
    }
    else {
        /* If the size of ClogC is too small for the new input matrix we   */
        /* re-initialize it.                                               */
        ClogCLen = mxGetNumberOfElements(ClogC_mxArray);
        
        if(ClogCLen<Ns*maxNt) {
            ClogC = mxGetPr(ClogC_mxArray);
            mxFree(ClogC);
            ClogC = create_ClogC(Ns, maxNt);
        }
        else
            ClogC = mxGetPr(ClogC_mxArray);
    }
   
    /* Computing entropy --------------------------------------------------*/
    totNt = 0;
    for(s=0; s<Ns; s++) {
        
        totNt += Nt[s];
        inv_Nt = 1 / Nt[s];
        
            NtpowL        = pow(Nt[s],L);
        inv_NtpowL        = 1 / NtpowL;
        inv_Ntpow_Lminus1 = 1 / pow(Nt[s],L-1);
        
        biasHlRS = 0;
        biasHlR  = 0;

        for(t=0;t<Nt[s];t++) {
            
            r = 0;
            rsh = 0;
            
            for(c=0;c<L;c++) {
                
                b = R[c + t*L + s*L*maxNt];
                if(doMap)
                    r += b * base[c];
                
                if(doSh) {
                    bsh = Rsh[c + t*L + s*L*maxNt];            
                    rsh += bsh * base[c];
                }
            
                if(doCrcs) {
                    indx = b + c*maxNb + s*maxNb*L;
                    ++Crcs[indx];
                    
                    if(doHlRS) {
                        HlRS[s] += increment(Crcs[indx], ClogC);
                        
                        /* Panzeri and Treves '96 bias correction          */
                        if(biasCorrNum==2 && t==Nt[s]-1) {
                            panzeri_treves_96(&Crcs[c*maxNb + s*L*maxNb], Nb[c], Nt[s], biasPr);
                            biasHlRS += *biasPr;
                        }
                    }
                }

                if(doHlR) {
                    indx = b + c*maxNb;
                    ++Crc[indx];
                    
                    *HlR += increment(Crc[indx], ClogC);
                    
                     /* Panzeri and Treves '96 bias correction             */
                    if(biasCorrNum==2 && t==Nt[s]-1) {
                        panzeri_treves_96(&Crc[c*maxNb], Nb[c], totNt, biasPr);
                        biasHlR += *biasPr;
                    }
                }
            }
            
            if(doCr) {
                ++Cr[r];
                                
                if(doHR) {
                    HR[0] += increment(Cr[r], ClogC);
                }
            }
            
            if(doHRS) {
                indx = r + s*Nr;
                ++Crs[indx];
                
                HRS[s] += increment(Crs[indx], ClogC);
            }
            
            if(doHshR) {
                ++Cshr[rsh];
                
                *HshR += increment(Cshr[rsh], ClogC);
            }
            
            if(doHshRS) {
                indx = rsh + s*Nr;
                ++Cshrs[indx];
                
                HshRS[s] += increment(Cshrs[indx], ClogC);
            }
        }
        
        if(doCirs) {

            /* This part is essentially identical to part of splitProds.   */
            /* However, since we want to perform the last product in the   */
            /* main function we need to repeat this piece of code (this has*/
            /* the advantage of greatly reducing the number of variables   */
            /* that need to be passed to teh subfunctions).                */
            
            /* If L=2 we just do the product directly                      */
            if(L==2) {
                
                K1_ptr = &Crcs[s*maxNb*L];
                K2_ptr = &Crcs[s*maxNb*L + maxNb];
                
                K1Len = Nb[0];
                K2Len = Nb[1];
            } else {

                /* First half                                              */
                N1 = L/2; /* floored ration                                */

                if(N1 == 1) {
                    /* If N1=1 we it must be L=3, thus we split the        */
                    /* products as                                         */
                    /*      P1 * (P2 * P3)                                 */

                    /* In this case for K1 we point directly to the        */
                    /* beginning of Crcs and its length is Nb[0]           */
                    K1_ptr = &Crcs[s*maxNb*L];
                    K1Len = Nb[0];
                } else {
                    X_ptr = &Crcs[s*maxNb*L];
                    M = maxNb;
    
                    n_ptr = &Nb[0];
                    
                    splitProds(X_ptr, M, N1, n_ptr, &K1_ptr, &K1Len);
                }

                X_ptr = &Crcs[s*maxNb*L + N1*maxNb];
                M = maxNb;
                N2 = L - N1;
                n_ptr = &Nb[N1];
    
                splitProds(X_ptr, M, N2, n_ptr, &K2_ptr, &K2Len);
            }

            /* We perform the last Kronecker product here                  */
            indx = 0;
            for(j=0; j<K2Len; j++) {
                for(i=0; i<K1Len; i++) {
                    
                    Cirs = K1_ptr[i] * K2_ptr[j];
                    
                    if(Cirs>0) {
                        if(doHiRS)
                            HiRS[s] += Cirs * log(Cirs);

                        if(doPir)
                            Pir[indx] += inv_Ntpow_Lminus1 * Cirs;
                    }

                    if(s==Ns-1) {
                        if(doHiR && Pir[indx]>0)
                            *HiR += Pir[indx] * log(Pir[indx]);

                        if(doChiR && Pir[indx]>0)
                            *ChiR +=  Cr[indx] * log(Pir[indx]);
                        
                    }
                    
                    ++indx;
                }
            }
            
            /* Freeing K1 and K2 when L=2 would erase Crcs                 */
            if(L> 2) {
                /* Freeing K1 when N1=1 would erase part of Crcs           */
                if(N1>1) {
                    mxFree(K1_ptr);
                }
                mxFree(K2_ptr);
            }
            
        }
        
        if(doHRS) {
            HRS[s] = (log(Nt[s]) - inv_Nt * HRS[s]) / log_2;
            
            /* Bias correction                                             */
            if(biasCorrNum==2) {
                panzeri_treves_96(&Crs[s*Nr], Nr, Nt[s], biasPr);
                HRS[s] += *biasPr;
            }
        }
        
        if(doHshRS) {
            HshRS[s] = (log(Nt[s]) - inv_Nt * HshRS[s]) / log_2;
            
            /* Bias correction                                             */
            if(biasCorrNum==2) {
                panzeri_treves_96(&Cshrs[s*Nr], Nr, Nt[s], biasPr);
                HshRS[s] += *biasPr;
            }
        }
        
        if(doHlRS) {
            HlRS[s] = (log(Nt[s]) * L - inv_Nt * HlRS[s]) / log_2;
            HlRS[s] += biasHlRS;
        }
        
        if(doHiRS)
            HiRS[s] = (log(NtpowL) - inv_NtpowL * HiRS[s]) / log_2;

    }
    
    inv_totNt  = 1 / totNt;
    
    if(doHR) {
        HR[0] = (log(totNt) - inv_totNt * HR[0]) / log_2;

        /* Bias correction                                                 */
        if(biasCorrNum==2) {
            panzeri_treves_96(Cr, Nr, totNt, biasPr);
            HR[0] += *biasPr;
        }
    }

    if(doHshR) {
        HshR[0] = (log(totNt) - inv_totNt * HshR[0]) / log_2;
        
        /* Bias correction                                                 */
        if(biasCorrNum==2) {
            panzeri_treves_96(Cshr, Nr, totNt, biasPr);
            HshR[0] += *biasPr;
        }
    }

    if(doHlR) {
        HlR[0]  = (log(totNt) * L - inv_totNt * HlR[0]) / log_2;
        HlR[0] += biasHlR;
    }

    if(doHiR)
        HiR[0]  = (log(totNt) - inv_totNt * HiR[0]) / log_2;

    if(doChiR)
        ChiR[0] = (log(totNt) - inv_totNt * ChiR[0]) / log_2;

    /* Freeing memory                                                      */
    if(doHR)    mxFree(Cr);
    if(doHRS)   mxFree(Crs);
    if(doHlR)   mxFree(Crc);
    if(doCrcs)  mxFree(Crcs);
    if(doPir)   mxFree(Pir);
    if(doHshR)  mxFree(Cshr);
    if(doHshRS) mxFree(Cshrs);
}

