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


/*
 *   Copyright (C) 2010 Cesare Magri
 *   Version: 6a
 */

/*
 * -------
 * 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.
 *
 */

/* ----------
 * DISCLAIMER
 * ----------
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
 * IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
 * THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
 * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
 * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
 * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
 * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
 * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */


/* 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;
}



/* ========================================================================*/
/* DIVIDE AND CONQUER                                                      */
/* ========================================================================*/
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                   */
    int     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)                       */
    
    bool     testMode;   /* specifies if working in test condition  (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  *biasPtr;
    double   biasHlRS;
    double   biasHlR;
    
    /* Output -------------------------------------------------------------*/
    
    double  *HR;
    double  *HRS;
    double  *HlR;
    double  *HlRS;
    double  *HiR;
    double  *HiRS;
    double  *ChiR;
    double  *HshR;
    double  *HshRS;
    
    /* Variables used for calling user-defined bias corrections            */
    mwSize  bufLen;
    char   *biasCorrFuncName;
    
    mxArray *lengthNrInputToBiasCorrFunc[1];
    double  *lengthNrInputToBiasCorrFuncPr;
    
    mxArray *lengthNbInputToBiasCorrFunc[1];
    double  *lengthNbInputToBiasCorrFuncPr;
    
    mxArray *outputOfBiasCorrFunc[1];
    
    mwIndex k;
    
    
    /* Reading input ------------------------------------------------------*/
    R  = mxGetPr(prhs[0]);
    
    PrToInputMxArray = mxGetField(prhs[1], 0, "Nt");
    Nt = mxGetPr(PrToInputMxArray);
    
    PrToInputMxArray = mxGetField(prhs[1], 0, "biasCorrNum");
    biasCorrNum = (int) *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);
    
    PrToInputMxArray = mxGetField(prhs[1], 0, "testMode");
    testMode = *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];
    /* If the users inputs a response matrix with a single stimulus then   */
    /* the number of dimensions for the response matrix will be 2 and      */
    /* asking for the number of elements in the third dimension will return*/
    /* a meaningless value. Thus first we need to check that there is      */
    /* indeed a third dimension. We do this using MXGETNUMBEROFDIMENSIONS  */
    if(mxGetNumberOfDimensions(prhs[0])>2)
        Ns = mxGetDimensions(prhs[0])[2];
    else
        Ns = 1;
    
    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) {
                    /* In test-mode we do not want to shuffle              */
                    if(testMode) {
                        Rsh[indx] = R[indx];
                    }
                    else {
                    /* 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];
    }
    
    if(biasCorrNum==-1) {
        /* Reading function name                                           */
        PrToInputMxArray = mxGetField(prhs[1], 0, "biasCorrFuncName");
        bufLen = mxGetN(PrToInputMxArray)*sizeof(mxChar) + 1;
        biasCorrFuncName = mxMalloc(bufLen);
        mxGetString(PrToInputMxArray, biasCorrFuncName, bufLen);
        
        if(doHlR || doHlRS) {
            lengthNbInputToBiasCorrFunc[0] = mxCreateDoubleMatrix(maxNb,  1, mxREAL);
            lengthNbInputToBiasCorrFuncPr  = mxGetPr(lengthNbInputToBiasCorrFunc[0]);
        }

        if(doHR || doHRS || doHshR || doHshRS) {
            lengthNrInputToBiasCorrFunc[0] = mxCreateDoubleMatrix(Nr, 1, mxREAL);
            lengthNrInputToBiasCorrFuncPr  = mxGetPr(lengthNrInputToBiasCorrFunc[0]);
            
            for(k=0; k<Nr; k++)
                lengthNrInputToBiasCorrFuncPr[k] = Cr[k];
        }
    }
    
    /* 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]);
    
    biasPtr = 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;
    biasHlR  = 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;

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

                        if(biasCorrNum!=0 && t==Nt[s]-1) {
                            switch(biasCorrNum) {
                                /* User-defined bias correction routines   */
                                case -1:
                                    for(k=0; k<maxNb; k++)
                                        lengthNbInputToBiasCorrFuncPr[k] = Crcs[c*maxNb + s*L*maxNb + k];

                                    mexCallMATLAB(1, &outputOfBiasCorrFunc[0], 1, &lengthNbInputToBiasCorrFunc[0], biasCorrFuncName);
                                    biasPtr = mxGetPr(outputOfBiasCorrFunc[0]);
                                    break;
                                
                                /* Panzeri and Treves '96 bias correction  */
                                case 2:
                                    panzeri_treves_96(&Crcs[c*maxNb + s*L*maxNb], Nb[c], Nt[s], biasPtr);
                                    break;
                            }
                            
                            biasHlRS += *biasPtr;
                        }
                    }
                }

                if(doHlR) {
                    indx = b + c*maxNb;
                    ++Crc[indx];
                    
                    *HlR += increment(Crc[indx], ClogC);
                    
                    if(biasCorrNum!=0 && s==Ns-1 && t==Nt[s]-1) {
                        switch(biasCorrNum) {
                            /* User-defined bias correction routines       */
                            case -1:
                                for(k=0; k<maxNb; k++)
                                    lengthNbInputToBiasCorrFuncPr[k] = Crc[c*maxNb + k];

                                mexCallMATLAB(1, &outputOfBiasCorrFunc[0], 1, &lengthNbInputToBiasCorrFunc[0], biasCorrFuncName);
                                biasPtr = mxGetPr(outputOfBiasCorrFunc[0]);
                                break;
                            
                            /* Panzeri and Treves '96 bias correction      */
                            case 2:
                                panzeri_treves_96(&Crc[c*maxNb], Nb[c], totNt, biasPtr);
                                break;
                        }
                        biasHlR += *biasPtr;
                    }
                }
            } /* End of loop on c ---------------------------------------- */
            
            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);
            }
        } /* End of loop on t -------------------------------------------- */
        
        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 the 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 ratio                                 */
                
                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!=0) {
                switch(biasCorrNum) {
                    /* User-defined bias correction routines               */
                    case -1:
                        /* Copying values to be passed to user-defined     */
                        /* biascorr functions (see "Creating Potential     */
                        /* Memory Leaks" under "Memory Management Issues"  */
                        /* in MATLAB help for an explanation of why the    */
                        /* the following method is used)                   */
                        for(k=0; k<Nr; k++)
                            lengthNrInputToBiasCorrFuncPr[k] = Crs[s*Nr + k];
                                                
                        mexCallMATLAB(1, &outputOfBiasCorrFunc[0], 1, &lengthNrInputToBiasCorrFunc[0], biasCorrFuncName);
                        biasPtr = mxGetPr(outputOfBiasCorrFunc[0]);
                        
                        break;
                        
                    /* Panzeri and Treves '96 bias correction              */
                    case 2:
                        panzeri_treves_96(&Crs[s*Nr], Nr, Nt[s], biasPtr);
                        break;
                }
                
                HRS[s] += *biasPtr;
            }
        }
        
        if(doHshRS) {
            HshRS[s] = (log(Nt[s]) - inv_Nt * HshRS[s]) / log_2;
            
            /* Bias correction                                             */
            if(biasCorrNum!=0) {
                switch(biasCorrNum) {
                    /* User-defined bias correction routines               */
                    case -1:

                        /* Copying values to be passed to user-defined     */
                        /* biascorr functions (see "Creating Potential     */
                        /* Memory Leaks" under "Memory Management Issues"  */
                        /* in MATLAB help for an explanation of why the    */
                        /* the following method is used)                   */

                        for(k=0; k<Nr; k++)
                            lengthNrInputToBiasCorrFuncPr[k] = Cshrs[s*Nr + k];

                        mexCallMATLAB(1, &outputOfBiasCorrFunc[0], 1, &lengthNrInputToBiasCorrFunc[0], biasCorrFuncName);
                        biasPtr = mxGetPr(outputOfBiasCorrFunc[0]);
                        break;
                        
                    /* Panzeri and Treves '96 bias correction              */
                    case 2:
                        panzeri_treves_96(&Cshrs[s*Nr], Nr, Nt[s], biasPtr);
                        break;
                }
                HshRS[s] += *biasPtr;
            }
        }
        
        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;

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

        /* Bias correction                                                 */
        if(biasCorrNum!=0) {
            switch(biasCorrNum) {
                /* User-defined bias correction routines                   */
                case -1:
                    for(k=0; k<Nr; k++)
                        lengthNrInputToBiasCorrFuncPr[k] = Cr[k];

                    mexCallMATLAB(1, &outputOfBiasCorrFunc[0], 1, &lengthNrInputToBiasCorrFunc[0], biasCorrFuncName);
                    biasPtr = mxGetPr(outputOfBiasCorrFunc[0]);
                    break;

                /* Panzeri and Treves '96 bias correction                  */
                case 2:
                    panzeri_treves_96(Cr, Nr, totNt, biasPtr);
                    break;
            }
        *HR += *biasPtr;
        }
    }

    if(doHshR) {
        HshR[0] = (log(totNt) - inv_totNt * HshR[0]) / log_2;

        /* Bias correction                                                 */
        if(biasCorrNum!=0) {
            switch(biasCorrNum) {
                /* User-defined bias correction routines                   */
                case -1:
                    for(k=0; k<Nr; k++)
                        lengthNrInputToBiasCorrFuncPr[k] = Cshr[k];

                    mexCallMATLAB(1, &outputOfBiasCorrFunc[0], 1, &lengthNrInputToBiasCorrFunc[0], biasCorrFuncName);
                    biasPtr = mxGetPr(outputOfBiasCorrFunc[0]);
                    break;

                /* Panzeri and Treves '96 bias correction                  */
                case 2:
                    panzeri_treves_96(Cshr, Nr, totNt, biasPtr);
                    break;
            }
            *HshR += *biasPtr;
        }
    }

    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;

    if(biasCorrNum==-1) {
        if(doHlR || doHlRS)
            mxDestroyArray(lengthNbInputToBiasCorrFunc[0]);
        
        if(doHR || doHRS || doHshR || doHshRS)
            mxDestroyArray(lengthNrInputToBiasCorrFunc[0]);
        
        mxDestroyArray(outputOfBiasCorrFunc[0]);
    }

    /* 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);
}

