/*=================================================================
 * Algorithme NIPALS
 *
 *============================================================*/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <math.h>
#include "mex.h"
#include "lapack.h"

/* ================================================== */
int isnul(double v) {
  /* isnul: retourne 1 si v ~= 0 
Permet la comparaison a zero des doubles avec une précision près */
  double eps=1e-37;
  if (fabs(v) < eps)
	   return (1);
  else
    return(0);
  } /* fin isnul */

/* ================================================== */
double ssq(int n, double *QQ) {
    /* ssq: calcul la somme des carrés des elts du vecteur QQ
INPUT
n : nombre d'elts de QQ
QQ: vecteur (n)
RETURN
la somme des carrés des elts (scalaire)
---------------------------------------------------- */
    int i;
    double som=0.0;
    for (i=0;i<n;i++) {
        som += (QQ[i]*QQ[i]);
    }
    return(som);
} /* fin ssq */
/* ================================================== */
void calcexpl(int N, int nth, double *Y, double *TT, double *CC,
        double * explY1c) {
    /* calcexpl: calcul des taux d'explication du tableau Y
INPUT
N : nombre de particules
nth: noombre de composantes
Y: réponses du modèle, centrées, réduites (N)
TT: (N X nth)
CC: (nth)
OUTPUT
explY1c: taux d'explication (nth)
---------------------------------------------------- */
    int i, h;
    double denoY, numY, indivi, explY1;
    denoY=ssq(N, Y);
    
    for (h=0; h<nth; h++) {
        numY=0.0;
        for (i=0; i<N; i++) {
            /* indivi=TT[i,h]*CC[h];*/
            indivi=TT[h*N+i]*CC[h];
            indivi=indivi*indivi;
            numY+=indivi;
        } /* fin i */
        explY1= (numY/denoY)*100;
        if (h==0) {
            explY1c[h]= explY1;
        } else {
            explY1c[h]=explY1c[h-1]+explY1;
        }
    } /* fin h */
    
} /* fin calcexpl */

/* ================================================== */
void  multMM(int n, int h, int k, double *TT, double *QQ, double *outp) {
    /* mumtMM:  multiplication d'une matrice par une autre
     * INPUT
     * n : nombre de lignes de la 1ier matrice
     * h: nombre de colonnes  de la 1ier matrice= nombre de lignes de la seconde
     * TT: matrice n X h
     * QQ: matrice h X k
     * OUTPUT
     * outp: matrice n X k,  préalablement allouée
     * ---------------------------------------------------- */
    int i,j, l;
    double a;
    
    for (i=0; i<n; i++) {
        for (j=0;j<k;j++) {
            /* outp[i, j]= somme(l=1:h)TT[i, l]*QQ[l,j] */
            a=0;
            for (l=0; l<h;l++) {
                a += (TT[l*n+i]*QQ[j*h+l]);
            }
            outp[j*n+i]=a;
        } /* fin j */
    } /* fin i */
    
} /* fin fonction multMM */
/* --------------------------------------------------- */
void multMt2(int n, int h,  double *TT, double *outp) {
    /* multMt2: multiplication de la transposée d'une matrice par elle-meme
     * INPUT
     * n : nombre de lignes
     * h: nombre de colonnes a considérer
     * nth: nombre de colonnes totales
     * TT: matrice n X h
     * OUTPUT
     * outp: pointeur sur la matrice en sortie h X h, préalablement allouée
     * ---------------------------------------------------- */
    int i,j, l;
    double a;
    
    for (i=0; i<h; i++) {
        for (j=0;j<h;j++) {
            /* outp[i, j]= somme(l=1:n)TT[l, i]*TT[l,j] */
            a=0;
            for (l=0; l<n;l++) {
                a += (TT[i*n+l]*TT[j*n+l]);
            }
            outp[j*h+i]=a;
        } /* fin j */
    } /* fin i */
    
} /* fin fonction multM2 */
/* --------------------------------------------------- */
int inverser(int h, double * A, double * outp) {
    /* Calcul  l'inverse d'une matrice carrée par LAPACK
     * INPUT
     * h: nombre de lignes et de colonnes de la matrice en entree
     * A: matrice en entree
     * OUTPUT
     * outp: matrice en sortie, allouée avant l'appel (h X h)
     * RETURN
     * le code d'erreur des sp lapack
     * ------------------------------------------------------- */
    
    /* DGE works in-place, so we copy the input first. */
    double *work;
    mwSignedIndex *iPivot;
    mxArray *Wwork, *mxPivot;
    mwSignedIndex  dims[2];
    ptrdiff_t info, lwork, lh;
    lh=(ptrdiff_t) h;
    info=0;
    memcpy(outp, A, h*h*sizeof(double));
    
    dims[0]= h+1;
    dims[1]=1;
    mxPivot = mxCreateNumericArray(2, dims, mxINT32_CLASS, mxREAL);
    if ( mxPivot==NULL) {
        mexPrintf("Fatal error: cannot allocate mxPivot\n");
        exit(EXIT_FAILURE);
    }
    
    iPivot = (mwSignedIndex*)mxGetData(mxPivot);
    dgetrf(&lh,&lh, outp,&lh,(ptrdiff_t *)iPivot,&info);
    if (info ==0) {
        Wwork = mxCreateDoubleMatrix(h, h, mxREAL);
        if ( Wwork==NULL) {
            mexPrintf("Fatal error: cannot allocate Wwork\n");
            exit(EXIT_FAILURE);
        }
        
        work=mxGetPr(Wwork);
        lwork=h*h;
        dgetri(&lh,outp,&lh,(ptrdiff_t *)iPivot,work,&lwork,&info);
        mxDestroyArray(Wwork);
    }
    mxDestroyArray(mxPivot);
    
    return(info);
    
}
/* ================================================== */
void centrereducpoly(double * Poly, int N, int nmono) {
    /* Centrer et réduire le tableau Poly
     *INPUT/OUTPUT
     * Poly: le tableau (nmono X N).
     * note: Poly[m,i] s'ecrit: Poly[i*nmono+m]
     * En sortie, il est centré et réduit
     * INPUT
     * N: nombre de particules
     * nmono: nombre de monomes
     * -------------------------------------------------------- */
    int i,j;
    double moyj, etj;
    
    /* calcul de la moyenne */
    for (j=0; j< nmono; j++) {
        moyj= 0.0;
        
        for (i=0; i<N; i++) {
            moyj += Poly[i*nmono+j];
            
        }
        moyj /= N;
        
        
        
        /* centrage des donnes */
        etj=0.0;
        for (i=0; i<N; i++) {
            Poly[i*nmono+j] = (Poly[i*nmono+j] - moyj);
            etj += (Poly[i*nmono+j]*Poly[i*nmono+j]);
            
        }
        
        etj = sqrt(etj / (N-1));
        
        /* Reduire les donnees */
        for (i=0; i<N; i++) {
            Poly[i*nmono+j] /= etj;
        }
        
    } /* fin j */
    
} /* fin centrereducpoly */
/*  ------------------------------------------------------- */
int carac(int N,  int nth,
        double *Y,
        double *TT, double *PP, double *CC,
        double *Q2Y1, double *Q2cumY1, double *explY1c ) {
    /* carac: calcul les Q2, Q2 cumulés et les taux d'explication
     * du tableau Y
     * INPUT
     * N:  nombre de particules = nombre de lignes de Y
     * nth: nombre de composantes
     * Y: vecteur des reponses (N).ATT: Sera centré-réduit en SORTIE
     * TT : (N X nth)
     * TT[h*N+i]  concerne la particule i et la composante h
     * PP: (nmono x Nbre de composantes voulues)
     * CC: (1 X Nbre de composantes voulues)
     * OUTPUT
     * Q2Y1: les Q2, vecteur (nth)
     * Q2cumY1: les Q2 cumulés, en % vecteur (nth)
     * explY1c: les taux d'explication, vecteur (nth)
     * RETURN
     * 0 ou un code d'erreur
     * ---------------------------------------------- */
    
    int h,i, itt, iyc1, kh, indh, indkh, flag, code=0;
    double a;
    double  rss, rpress, rpress2, hat, rusue1, prod;
    mxArray *TTwork, *invAwork, *Ayc1, *tempo;
    double *pTTwork, *pinvAwork, *yc1, *ptempo;
    
    
    /* Centrer, reduire la reponse */
    centrereducpoly(Y,  N, 1);
    
    /* on alloue dans ce sens yc1, pour ne pas avoir a allouer pour chaque valeur de h */
    Ayc1= mxCreateDoubleMatrix(nth, N,  mxREAL);
    if (Ayc1==NULL) {
        mexPrintf("Fatal error: cannot allocate Ayc1\n");
        exit(EXIT_FAILURE);
    }
    yc1 =  mxGetPr(Ayc1);
    /* yc1[i,j] est yc1[j*nth+i] , donc yc1[a,b] du papier est yc1[a*nth+b] */
    for(h=1; h<=nth; h++) {
        indh=h-1; /* les indices commencent a 0 */
        for (i=0;i<N;i++) {
            itt=indh*N+i;
            iyc1 = i*nth+indh;
            a=CC[indh] * TT[itt];
            if (h>1) {
                yc1[iyc1] = yc1[i*nth+(indh-1)]+ a; /* =yc1[i,h-1]+a du papier */
            }      else
                yc1[iyc1] =a;
        } /* fin i */
        
    } /* fin h */
    /* --- Fin calcul des Yi chapeaux --- */
    
    /* --- Calcul des Q2 a partir des th --- */
    
    for(h=1; h<=nth; h++) {
        indh=h-1; /* les indices commencent a 0 */
        rss=rpress2=0;
        TTwork = mxCreateDoubleMatrix(h, h, mxREAL);
        if (TTwork==NULL) {
            mexPrintf("Fatal error: cannot allocate TTwork\n");
            exit(EXIT_FAILURE);
        }
        pTTwork = mxGetPr(TTwork);
        multMt2(N, h, TT,pTTwork);
        invAwork= mxCreateDoubleMatrix(h, h, mxREAL);
        if (invAwork==NULL) {
            mexPrintf("Fatal error: cannot allocate invAwork\n");
            exit(EXIT_FAILURE);
        }
        
        pinvAwork = mxGetPr(invAwork);
        code=inverser(h, pTTwork, pinvAwork);
        if (code !=0) {
            return(code);
        } /* fin code !=0 */
        mxDestroyArray(TTwork);
        tempo=  mxCreateDoubleMatrix(N, h,  mxREAL);
        if (tempo==NULL) {
            mexPrintf("Fatal error: cannot allocate tempo\n");
            exit(EXIT_FAILURE);
        }
        ptempo= mxGetPr(tempo);
        multMM( N, h, h, TT, pinvAwork, ptempo);
        mxDestroyArray(invAwork);
        for (i=0; i<N; i++) {
            /* rusue1 = yc1[i,h]- y[i];  */
            rusue1 = yc1[i*nth+indh] - Y[i];
            hat=0;
            for (kh=1; kh<=h; kh++) {
                indkh=kh-1;
                hat += (ptempo[N*indkh+i]*TT[N*indkh+i]);
            } /* fin kh */
            rpress= rusue1 /(1-hat);
            rpress2 += (rpress*rpress);
            if (h>1) {
                /* rss += (yc1[i, h-1]-Y[i])**2; */
                a=  yc1[i*nth+indh-1] -Y[i];
                rss += (a*a);
            } /* fin h>0 */
        } /* fin i */
        mxDestroyArray(tempo);
        
        if (h==1) {
            rss= N-1;
        }
        Q2Y1[indh]= 1.0 - (rpress2/rss);
        
    } /* fin h */
    mxDestroyArray(Ayc1);
    
    /* --------- Construction des Q2 cumulés ----- */
    flag=0;
    for (h=0; h<nth; h++) {
        if (Q2Y1[h] <=0) flag++;
        Q2cumY1[h] = 0.0;
    }
    if (flag < nth) {
        prod=1;
        for (h=0; h<nth; h++) {
            if (Q2Y1[h] >0) {
                a= 1-Q2Y1[h];
                prod *= a;
                Q2cumY1[h]=1-prod;
                if ( (h>0) && (Q2cumY1[h] <=0)) {
                    Q2cumY1[h] = Q2cumY1[h-1];
                }
                /* on met les Q2cum, qui sont [0,1] entre [0,100] */
                Q2cumY1[h]=Q2cumY1[h]*100.0;
                
            } /* fin (Q2Y1[h] >0) */
        } /* fin h */
        for (h=1; h<nth; h++) {
	  if (Q2cumY1[h] <=0) {
	    Q2cumY1[h] = Q2cumY1[h-1];
	  }
	}
        
    } /* fin (flag < nth) */
    /* --------- Calcul des taux d'explication du tableau Y   ----- */
    calcexpl( N,  nth, Y, TT, CC, explY1c);
    return(code);
} /* fin fonction carac */

/* ================================================== */
double multA(double A, double *X, int i) {
    /* Multiplier un élément du tableau en entrée par une valeur.
     *  Fonction appelée par generPoly
     * INPUT
     * A: valeur scalaire
     * X: tableau
     * i: indice de l'élément de X à multiplier par A
     * OUTPUT
     * A*X[i]
     *-------------------------------------------------- */
    return (A*X[i]);
} /* fin multA */
/* ================================================== */
void generPoly(double *X, int p, int degre, double *res)  {
    /* Generation d'un polynome complet de degré 2 ou 3 ou 4
     * avec toutes les interactions et sans constantes
     * INPUT
     * X: vecteur des parametres (p valeurs)
     * degre: degré du polynome, doit être entre 2 et 4
     * p: le nombre de parametres
     * OUTPUT
     * res: tableau alloué avant l'appel
     * TEST
     * p=4; degre=2; N=2;
     * X=[1:4; 4 3 2 1]; XT=transpose(X);
     * a=generPoly(p, N, degre, XT);
     * On doit obtenir:
     * 1 1 2 3 4 2 4 6 8 3 9 12 4 16
     * 4 16 12 8 4 3 9 6 3 2 4 2 1 1
     * Sortie du test:
     * Les termes du polynome dans l'ordre, si p=4, degre=2:
     * X1, X1², X1.X2 X1.X3 X1.X4 X2 X2² X2.X3, X2.X4 X3 X3.X3, etc...
     *Les termes du polynome dans l'ordre, si p=4, degre=3:
     * X1 X1.X1 X1.X1.X1 X1².X2 X1².X3  X1².X4
     *   X1.X2 X1.X2.X2 X1.X2.X3 X1.X2.X4
     *     X1.X3 X1.X3² X1.X3.X4
     *       X1.X4 X1.X4²,
     * X2 X2² X2².X3 X2².X4 X2.X3 X2.X3.X3 X2.X3.X4 , etc...
     * ---------------------------------------------------------- */
    
    double un=1;
    int ires=0; /* indice dans le tableau resultat */
    int i, j,l,m;
    /* ..................................... */
    for (i=0; i<p; i++) {
        res[ires++]=multA(un, X, i);

        for (j=i; j<p; j++) {
            res[ires++]=multA(X[i],  X, j);
            if (degre >2) {
                for (l=j; l<p; l++) {
                    res[ires++]=multA(X[i]*X[j], X, l);
                    if (degre >3) {
                        for (m=l; m<p; m++) {
                            res[ires++]=multA(X[i]*X[j]*X[l], X, m);
                        }
                    } /* fin (degre >3) */
                }
            } /* fin (degre >2) */
        }
    }
    
    /* ..................................... */
    
} /* fin generPoly */

/* ================================================== */

double moy(double *X, int n) {
    /* Calculer la moyenne de X
     * INPUT
     * X: vecteur (n)
     * n: longueur de X
     * OUTPUT
     * la moyenne de X
     *-------------------------------------------------- */
    int i;
    double moyx=0.0;
    for (i=0; i<n; i++) {
        moyx += X[i];
    }
    moyx /= n;
    return(moyx);
} /* fin moy */
/* ================================================== */

double sd(double *X, int n) {
    /* Calculer l'ecart-type de X
     * INPUT
     * X: vecteur (n)
     * n: longueur de X
     * OUTPUT
     * l'écart-type de X
     * NOTE: on renvoie la racine carrée de la somme des carrés des éléments divisée par (n-1)
     *-------------------------------------------------------- */
    int i;
    double etx=0.0;
    for (i=0; i<n; i++) {
        etx += ( X[i]* X[i]);
    }
    etx = sqrt(etx / (n-1));
    return(etx);
}
/* ================================================== */

double correlation(double *X, double *Y, int n) {
    /* Calculer le coefficient de correlation de X et Y de longueur n
     * INPUT
     * X: vecteur (n)
     * Y:  vecteur (n)
     * n: nombre d'éléments de X et de Y
     * OUTPUT
     * le coefficient de correlation de X et Y
     * TEST
     * X= [2.9080 0.8252 1.3790]; Y= [-1.0582 -0.4686 -0.2725];
     * n=3; a=correlation(X,Y,n)
     ** a comparer à:
     * [r,p] =corrcoef(X,Y);
     * r(1,2) doit etre égal à -0.8768
     * --------------------------------------------------------- */
    
    int i;
    double moyx=0.0, moyy=0.0, etx=0.0, ety=0.0, a, b, vcorr=0.0;
    
    /* ..................................... */
    /* calcul des moyennes */
    for (i=0; i<n; i++) {
        moyx += X[i];
        moyy += Y[i];
    }
    moyx /= n;
    moyy /= n;
    
    /* ..................................... */
    /* covariance et écarts-types*/
    for (i=0; i<n; i++) {
        a= X[i] - moyx;
        etx += (a*a);
        b= Y[i] - moyy;
        ety += (b*b);
        vcorr += (a*b);
    }
    vcorr=vcorr/(n-1);
    etx = sqrt(etx / (n-1));
    ety = sqrt(ety / (n-1));
    
    /* ..................................... */
    /* covariance/ produit des écarts-types= coefficient de correlation */
    vcorr = vcorr / (etx*ety);
    
    
    return(vcorr);
    
} /* fin correlation */

/* ================================================== */
void nipals1( int nth, int nmono, int N,
        double *E, double *F,
        double *TT, double *WW,
        double *uh, double *PP, double *CC,
        double *wh ) {
    /* Algorithme NIPALS
     * INPUT
     * nth: nombre de composantes voulues
     * nmono: nombre de monomes
     * N: nombre de particules
     * INPUT-OUTPUT
     * E: tableau des monomes centré, réduit (N x nombre de monomes)
     * Valeurs rangées par colonnes:
     *  E[j*N+i] concerne le monome j et la particule i
     * F: réponse centrée, réduite à un temps donné (N)
     * NOTE: E,F sont modifiés en sortie
     * OUTPUT
     * TT: (N X Nbre de composantes voulues)
     * WW: (nmono x Nbre de composantes voulues)
     * Les valeurs sont rangées par colonnes:
     * TT[h*N+i]  concerne la particule i et la composante h
     * WW[h*nmono+j] concerne le monome j et la composante h
     * WORKING
     * uh: (N )
     * PP: (nmono x Nbre de composantes voulues)
     * CC: (1 X Nbre de composantes voulues)
     * wh : (nmono)
     * ---------------------------------------------- */
    
    int i, j, h, itt;
    double whj, uhj, thj, chj, chj2,  phj, lasqrt;
    
    /* ..................................... */
    /* Pour mémoire: exemples d'allocation en C/matlab:
     * TheString = mxCalloc(yLength, sizeof(char));
     * mxCreateDoubleMatrix(colLen, rowLen, mxREAL);
     * Call the Function built-in "conv":
     * mexCallMATLAB(1,&result,2,arguments,"conv");
     */
    /* ..................................... */
    for (i=0; i< N; i++) {
        uh[i]=F[i];
    }
    
    for (h=0; h<nth; h++) {
        
        /* ..... Calcul des wh ..... */
        uhj=0;
        for (i=0; i< N; i++) {
            uhj += (uh[i]* uh[i]);
        }
        for (j=0; j < nmono; j++) {
            whj=0;
            for (i=0; i< N; i++) {
                whj += (E[j*N+i] * uh[i]);
            }
            whj = whj/ uhj;
            wh[j]= whj;
            
        } /* fin j */
        
        
        lasqrt=0.0;
        for (j=0; j < nmono; j++) {
            lasqrt += (wh[j]*wh[j]);
        }
        lasqrt =sqrt(lasqrt); /*  sqrt est dans math.h */
        
        
        
        whj=0.0;
        for (j=0; j < nmono; j++) {
            
            wh[j] =  wh[j] / lasqrt;
            
            WW[h*nmono+j]= wh[j];
            whj += (wh[j]*wh[j]);
        }

        
        /* ..... Calcul des th ..... */
        for (i=0; i< N; i++) {
            thj=0;
            for (j=0; j < nmono; j++) {
                thj += (E[j*N+i] * wh[j]);
            }

            thj = thj / whj;
            TT[h*N+i] = thj;
        } /* fin i */
        /* .....  Calcul des ch  ..... */
        chj=0; thj=0;
        for (i=0; i< N; i++) {
            itt=h*N+i;
            chj += (F[i] * TT[itt]);
            thj += (TT[itt] *TT[itt]);
        }
        
        chj = chj/thj;
        CC[h] = chj;
         
        /* .....  Calcul des uh  ..... */
        chj2 = CC[h] * CC[h];
        
        
        for (i=0; i< N; i++) {
            uhj= F[i] * CC[h];
            uhj = uhj/chj2;
            uh[i] = uhj;
        } /* fin i */
        
        
        /* .....  Calcul des PP  ..... */
        for (j=0; j < nmono; j++) {
            phj=thj=0;
            for (i=0; i< N; i++) {
                itt=h*N+i;
                phj += (E[j*N+i] * TT[itt]);
                thj += (TT[itt] * TT[itt]);
            }
            phj = phj/thj;
            PP[h*nmono+j] = phj;
        }
        
        /* ..... Modifier E et F ..... */
        for (i=0; i< N; i++) {
            for (j=0; j < nmono; j++) {
                E[j*N+i] = E[j*N+i] - (TT[h*N+i] * PP[h*nmono+j]);
            }
            F[i] = F[i] - (TT[h*N+i] * CC[h]);
        }
         
        
    } /* fin h */
     
} /* fin nipals1 */

/* ================================================== */

void vippls(
        int nmono, int nth, int N,
        double *F, double *T, double *W,
        double *VIP,
        double *coryt, double * redonk) {
    /* INPUT
     * nmono: nombre de monomes
     * nth: nombre de composantes
     * N: nombre de particules
     * F: réponse centrée réduite (N)
     * T: (N X nombre de composantes)
     * W (nombre de monomes X nombre de composantes)
     * OUTPUT
     * VIP (nombre de monomes)
     * WORKING
     * coryt: (N)
     * redonk: (nombre de composantes)
     ** ------------------------------------------ */
    
    int i, iw, k, j;
    double  vcork, totw=0.0, redont =0.0;
    
    
    /* ..... Calcul du carré des correlations de T et F .... */
    for (k=0; k<nth; k++) {
        for (i=0; i< N; i++) {
            coryt[i]= T[k*N+i];
        }
        
        vcork = correlation( F, coryt, N);
        redonk[k] = (vcork * vcork);
        redont += redonk[k];
    } /* fin k */
    
    /* ..... Calcul des VIP ..... */
    for (i=0; i<nmono; i++) {
        totw=0.0;
        for (j=0; j< nth; j++) {
            iw=j*nmono+i;
            totw += (redonk[j] * (W[iw] * W[iw]));
        }
        VIP[i] = sqrt(((double)nmono / redont) * totw);
        
    } /* fin i */
    
    
} /* fin vippls */

/* ================================================== */
void sipls0(double *vipin, int nmono,
        unsigned int *indic, int np,
        double *tsivip,
        double * sivip) {
    /* INPUT
     * vipin: VIP de chaque monome (nmono)
     * nmono: nombre de monomes
     * indic: tableau de 0 et 1 (nmono X np)
     * Les valeurs sont rangées paramètre par paramètre:
     *  l'indication du paramètre j dans le monome i est dans indic[j*nmono+i]
     * np: nombre de paramètres
     * OUTPUT
     * tsivip: indices SIVIP totaux en absolu et en % (np X 2)
     * WORKING
     * sivip: (nmono)
     * -------------------------------------------- */
    
    int i, j, nonnul=0;
    double som;
    
    /* ..... Compter le nbre de vipin non nuls ..... */
    for (i=0; i< nmono; i++) {
      if (isnul(vipin[i]) ==0) { /* (vipin[i] != 0) */
            nonnul++; /* comparaison a zero d'un double a eps pres */
        }
    }
    
    /* ..... Calcul du carré des vipin divisé par le nbre de non nuls ..... */
    for (i=0; i< nmono; i++) {
        sivip[i] = (vipin[i]*vipin[i])/nonnul;
    }
    
    /* ..... SIVIP totaux en absolu .... */
    
    som=0.0;
    for (j=0; j<np; j++) {
        tsivip[j]=0.0;
        for (i=0; i<nmono; i++) {
            tsivip[j] += (sivip[i]*indic[j*nmono+i]);
        }
        som += (tsivip[j]);
    }
    
    /* ..... SIVIP totaux en pourcentage dans la 2ieme colonne de tsivip .... */
    for (j=0; j<np; j++) {
        tsivip[np+j]= (tsivip[j]/som) *100.0;
    }
    
} /* fin sipls0 */


/* ================================================== */

void calcIST(double *X, double *Y, int np, int N, int ntemps, int degre,
        int nmono,  int nth, unsigned int *indic,
        double *IST, double *PIST, double * Q2Y1, double *Q2cumY1,
        double *explY1c,
        double *Poly, double *E, double *F,double *TT, double *WW, double *uh,
        double *PP, double *CC, double *wh,
        double *sivip, double *tsivip, double *vip, double *coryt, double *redonk) {
    /* PROGRAMME PRINCIPAL de calcul des np IST de t à tmax (ntemps temps)
     * INPUT
     * X: valeurs des np paramètres sur les N particules (N X np)
     *      X[np*i +j] concerne la particule i, et le parametre j,
     *      c.a.d toutes les valeurs qui concernent une particule sont contigues
     * Y: valeurs du modèle calculées avec les valeurs des paramètres X (N, ntemps)
     *     Y[N*itemps+i] concerne la particule i au temps itemps
     *      Toutes les valeurs qui concernent un temps sont contigues
     * np: nombre de parametres
     * N: nombre de particules
     * ntemps: nombre de temps de calcul des IST
     * degre: degre du polynome (2 a 4)
     * nmono: nombre de monomes
     * nth: nombre de composantes
     * indic: tableau de 0 et 1 decrivant la présence de chaque paramètre dans chaque monome;
     * Les valeurs sont rangées paramètre par paramètre:
     *  l'indication du paramètre j dans le monome i est dans indic[j*nmono+i]
     * OUTPUT
     * IST (np, ntemps) en absolu
     * PIST (np, ntemps) en pourcentage
     * Q2Y1 (nth, ntemps): les Q2
     * Q2cumY1 (nth, ntemps): les Q2 cumulés
     * explY1c (nth, ntemps): les taux d'explication de Y
     * WORKING
     * Poly: tableau (nbre de monomes X N):  les monomes centrés-réduits
     * E: =Poly en entrée de nipals1 et modifié en sortie de nipals1
     * F: vecteur (N)= Y centré-réduit, modifié en sortie de nipals1
     * TT: (N X Nbre de composantes voulues)
     * WW: (nmono x Nbre de composantes voulues)
     * uh: (N )
     * PP, WWS: (nmono x Nbre de composantes voulues)
     * CC: (1 X Nbre de composantes voulues)
     * wh : (nmono)
     * sivip: vecteur (nbre de monomes) destiné à contenir les sivip
     * tsivip: tableau (np X 2) destiné à contenir les sivip totaux a un temps donné:
     *  1iere colonne: en absolu, 2ieme colonne: en pourcentage
     * vip (nbre de monomes)
     * coryt: (N)
     * redonk: (nombre de composantes)
     * -------------------------------------------- */
    int i, j, ii, itemps, code;
    double moyY, etY;
    
    /* ..................................... */
    
    
    /* ...... Generation des monomes pour chaque ligne de X ..... */
    /* ATTENTION: les valeurs dans Poly pour une particule
     * sont rangées à la queue leu. Poly est un tableau (nmono,N):
     * i.e le monome m pour la particule i est Poly[i*nmono+m]
     * En ce qui concerne X: X[np*i +j] concerne la particule i, et le parametre j
     */
    
    for (i=0; i<N; i++) {
        generPoly((X+ (np*i)), np, degre, Poly+(nmono*i));
        
    }
    
    
    /* ..... Centrer et réduire le tableau des monomes ..... */
    centrereducpoly( Poly,  N, nmono);
    
    
    /* ..... Boucle sur le temps ..... */
    for (itemps=0; itemps < ntemps; itemps++) {
        
        /* ..... Centrer et réduire le vecteur réponse dans F..... */
        
        moyY= moy((Y+ (itemps*N)), N); /* moyenne des Y au temps "itemps" */
        
        for (i=0; i<N; i++) {
            
            F[i] = (Y[itemps*N+i]- moyY);
        }
        
        etY= sd(F, N);
        for (i=0; i<N; i++) {
            F[i] /= etY;
        }
        
        /* ..... Appel de nipals1:   renvoie T et W ..... */
        /* remettre Poly dans E car nipals1 a modifié E au temps précédent */
        
        ii=0;
        for (i=0; i<N; i++) {
            for (j=0; j<nmono; j++) {
                E[j*N+i]= Poly[ii++];
            }
            
        }
        
        
        nipals1(  nth,  nmono,  N,
                E, F, TT, WW,  uh, PP,  CC,  wh );
        
        /* ..... Appel de vippls ..... */
        
        /* a besoin de la réponse non centrée-réduite, de TT et WW; renvoie vip */
        vippls(nmono, nth,  N,
                (Y+ (itemps*N)), TT, WW, vip, coryt, redonk);
        /* VERIF AU CASOU
         * pour verifier sum(vip**2) doit etre egale au nombre de monomes
         *
         * etY=0;
         * for (i=0; i<nmono; i++) {
         * etY +=(vip[i]*vip[i]);
         * }
         * mexPrintf(" somme des vip2 %g\n", etY);
         * FIN verif */
        
        
        /* ..... Calcul des Q2 ..... */
        /* carac va modifier Y (centrer-reduction) */
        code=carac( N, nth, (Y+itemps*N),
                TT, PP, CC,
                (Q2Y1+itemps*nth), (Q2cumY1+itemps*nth),
                (explY1c+itemps*nth));
        /* VERIF AU CASOU
        if (code !=0) {
            mexPrintf("calcIST: Probleme dans carac %d\n", code);
        }
	*/
        
        /* ..... Appel de sipls0: renvoie tsivip ..... */
        
        sipls0(vip,  nmono, indic, np,  tsivip, sivip);
        
        for (i=0; i<np; i++) {
            IST[np*itemps +i] = tsivip[i]; /*tisvp en absolu */
            PIST[np*itemps +i] = tsivip[np+i]; /*tisvp en pourcentage */
            
        }
        
    } /* ------------- fin itemps ------------------*/
    
    
} /* fin calcIST */


/* ================================================== */
/* The gateway function:
 * INPUT
 * nlhs: le nombre de sorties
 * nrhs: le nombre d'arguments d'entrée
 * prhs: les arguments d'entrée, qui sont:
 *   np: nombre de parametres
 *   N: nombre de particules
 *   ntemps: nombre de temps de calcul des IST
 *   degre: degre du polynome (2 a 4)
 *   nth: nombre de composantes voulues
 *   nmono: nombre de monomes
 *   X: valeurs des np paramètres sur les N particules (N X np)
 *   Y: valeurs du modèle calculées avec les valeurs des paramètres X (N, ntemps)
 *   indic: tableau de 0 et 1 decrivant la présence de chaque paramètre dans chaque monome
 *   Poly: tableau (nbre de monomes X N)  destiné à contenir les monomes
 *   E: tableau (nbre de monomes X N) de sauvegarde de Poly
 *   F: vecteur (N) de sauvegarde de Y
 *   TT: (N X Nbre de composantes voulues)
 *   WW: (nmono x Nbre de composantes voulues)
 *   uh: (N )
 *   PP: (nmono x Nbre de composantes voulues)
 *   CC: (1 X Nbre de composantes voulues)
 *   wh : (nmono)
 *   sivip: vecteur (nbre de monomes) destiné à contenir les sivip
 *   tsivip: tableau (np X 2) destiné à contenir les sivip totaux a un temps donné
 *   vip (nbre de monomes)
 *   coryt: (N)
 *   redonk: (nombre de composantes)
 *   np: nombre de parametres
 *   N: nombre de particules
 *   ntemps: nombre de temps de calcul des IST
 *   degre: degre du polynome (2 a 4)
 * OUTPUT
 * plhs: les sorties, qui sont:
 *   les IST (np, ntemps): indices totaux en absolu
 *   les PIST (np,ntemps): indices totaux en pourcentage
 *   Q2Y1 (nth, ntemps)
 *   Q2cumY1 (nth, ntemps): les Q2 cumulés
 *   explY1c (nth, ntemps): les taux d'explication de Y
 * ------------------------------------------------- */
void mexFunction( int nlhs, mxArray *plhs[],
        int nrhs, const mxArray *prhs[] )
{
    int  i, N, np, ntemps, degre, nmonomes,  nth;
    double *X, *Y, *Poly, *sivip, *tsivip, *IST, *PIST,
            *Q2Y1, *Q2cumY1, *explY1c,
            *vip, *coryt, *redonk;
    double *E, *F, *TT,  *WW, *uh, *PP, *CC, *wh;
    unsigned int *indic;
    mxArray *a;
    
    
    /* ..... Create a pointer to the  input   .... */
    i=0;
    /* les scalaires */
    np=(int ) mxGetScalar(prhs[i++]);
    N=(int ) mxGetScalar(prhs[i++]);
    ntemps = (int ) mxGetScalar(prhs[i++]);
    degre=(int) mxGetScalar(prhs[i++]);
    nth=(int) mxGetScalar(prhs[i++]);
    nmonomes= (int) mxGetScalar(prhs[i++]);
    
    /* les tableaux */
    X = mxGetPr(prhs[i++]);
    Y=   mxGetPr(prhs[i++]);
    indic = (unsigned int *)mxGetPr(prhs[i++]);
    /* le tableau de travail, destiné à contenir les monomes */
    Poly= mxGetPr(prhs[i++]);
    /* tableaux pour nipals1 */
    E= mxGetPr(prhs[i++]);
    F= mxGetPr(prhs[i++]);
    TT= mxGetPr(prhs[i++]);
    WW= mxGetPr(prhs[i++]);
    uh= mxGetPr(prhs[i++]);
    PP= mxGetPr(prhs[i++]);
    CC= mxGetPr(prhs[i++]);
    wh= mxGetPr(prhs[i++]);
    
    /* le tableau de travail, destiné à contenir les sivip */
    sivip = mxGetPr(prhs[i++]);
    /* le tableau de travail, destiné à contenir les tsivip */
    tsivip = mxGetPr(prhs[i++]);
    /* autres tableaux de travail */
    vip = mxGetPr(prhs[i++]);
    coryt = mxGetPr(prhs[i++]);
    redonk = mxGetPr(prhs[i++]);
    
    /* ..................................... */
    /* ...... Create the output matrix ..... */
    
    plhs[0] = mxCreateDoubleMatrix((mwSize)(np),(mwSize)(ntemps),mxREAL);
    plhs[1] = mxCreateDoubleMatrix((mwSize)(np),(mwSize)(ntemps),mxREAL);
    for (i=2; i<=4; i++) {
        a = mxCreateDoubleMatrix((mwSize)(nth), (mwSize)(ntemps),mxREAL);
        if (a==NULL) {
            mexPrintf("Allocation impossible\n");
            exit(EXIT_FAILURE);
        } else {
            plhs[i] =a;
            
        }
        
    }
    
    /* get a pointer to the real data in the output matrix */
    IST = mxGetPr(plhs[0]);
    PIST = mxGetPr(plhs[1]);
    Q2Y1 = mxGetPr(plhs[2]);
    Q2cumY1 = mxGetPr(plhs[3]);
    explY1c =     mxGetPr(plhs[4]);
    /* ...... Call the computational routine  .....*/
    
    calcIST(X, Y,   np, N, ntemps, degre, nmonomes,  nth, indic,
            IST, PIST,  Q2Y1, Q2cumY1,  explY1c,
            Poly, E, F, TT, WW, uh, PP, CC, wh,
            sivip, tsivip, vip, coryt, redonk) ;
    
    /* ..................................... */
    
} /* fin mexFunction */




