/*=================================================================
 *  Generer le tableau indic
 *  * COMPILATION sous matlab par:
 * mex SRC/SIMULATION_STRATEGIE/C/generIndic.c
 *============================================================*/

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

/*=================================================================*/
int presence2(int type, int i, int j) {
    /* retourne 1 ou 0 selon que le monome à 2 variables doit etre considéré selon la valeur de type */
    
    switch (type) {
        case 1:
            /* poly complet */
            return 1;
        case 2:
            /* sans carrés */
            if  (i==j) {return 0; } else { return 1; }
        case 3:
            /* sans interactions */
            if  (i!=j) {return 0; } else { return 1; }
    }
}


/*=================================================================*/
int presence3(int type, int i, int j, int l) {
    /* retourne 1 ou 0 selon que le monome à 3 variables doit etre considéré selon la valeur de type */
    
    switch (type) {
        case 1:
            /* poly complet */
            return 1;
        case 2:
            /* sans carrés */
            if  ((i==j) || (i==l) || (j==l)) {return 0; } else { return 1; }
        case 3:
            /* sans interactions */
            if  ((i!=j) || (i!=l) || (j!=l)) {return 0; } else { return 1; }
    }
}


/*=================================================================*/
int presence4(int type, int i, int j, int l, int m) {
    /* retourne 1 ou 0 selon que le monome à 4 variables doit etre considéré selon la valeur de type */
    
    switch (type) {
        case 1:
            /* poly complet */
            return 1;
        case 2:
            /* sans carrés */
            if  ((i==j) || (i==l) ||  (i==m) ||
                    (j==l) || (j==m) || (l==m)) {return 0; } else { return 1; }
        case 3:
            /* sans interactions */
            if  ((i!=j) || (i!=l) ||  (i!=m) ||
                    (j!=l) || (j!=m) || (l!=m)) {return 0; } else { return 1; }
            
    }
}


/*=================================================================*/
void generIndic(int nmono, int p, int degre, int type, unsigned int *indic) {
    /* Generer le tableau indic
     * ENTREES
     * nmono: nombre de monomes
     * p: nombre de paramètres
     * degre: degré du polynome, doit être entre 2 et 4
     * type: type du polynome:
     * 1: complet, 2: sans carrés, 3; sans interaction
     * ENTREE-SORTIE
     * indic: tableau de 0 et 1 (nmono X p) initialisé à 0
     * indic[i,j]=1, si le parametre j est présent dans le monome i
     * Les valeurs sont rangées en colonnes
     * indic[i,j] => indic[j*nmono+i]
     * ----------------------------------------------- */
    /* Principe:
on considère que les monomes sont dans l'ordre :
(p=3, degre=2)
X1, X1*X1, X1*X2, X1*X3, X2, X2*X2, X2*X3, X3, X3*X3
 (p=3, degre=3)
X1, X1*X1, X1*X1*X1, X1*X1*X2,X1*X1*X3, X1*X2, X1*X2*X2, X1*X2*X3,
X1*X3,  X1*X3*X3,
X2, X2*X2, X2*X2*X2,X2*X2*X3, X2*X3*X3,
X3, X3*X3,  X3*X3*X3
Dans les boucles qui suivent, i correspond à l'indice du 1ier membre
du monome, j au 2ieme, l au 3ieme
     */
    
    
    int val, ires=0; /* indice dans le tableau resultat */
    int i, j,l,m;
    
    for (i=0; i<p; i++) {
        
        indic[i* nmono+ires]=1;
        ires++;
        for (j=i; j<p; j++) {
            switch (degre) {
                case 2:
                    val=presence2(type, i,j);
                    break;
                case 3:
                       /* without 2 input interactions or without interactions*/
                    if ((type==2) || (type==4))
                        val=presence2(3,i,j);
                    else
                        val=presence2(1,i,j);
                    
                            break;
                case 4:
                    if ((type==2) ||  (type==5))
                        /* without 2 input interactions or without interactions*/
                        val=presence2(3,i,j);
                    else
                        val=presence2(1,i,j);
                    
                            break;
            } /* fin switch */
            
            
            indic[i* nmono+ires]=val;
            indic[j* nmono+ires]=val;
            ires++;
            if (degre >2) {
                for (l=j; l<p; l++) {
                    switch (degre) {
                        case 3:
                            if ( (type==3)|| (type==4))
                                /* without 3 input interactions or without interactions*/
                                val=presence3(3, i,j,l);
                            else
                                val=presence3(1, i,j,l);
                            break;
                        case 4:
                            if ( (type==3)|| (type==5))
                                /* without 3 input interactions or without interactions*/
                                val=presence3(3, i,j,l);
                            else
                                val=presence3(1, i,j,l);
                            break;
                    } /* fin switch */
                    
                    indic[i* nmono+ires]=val;
                    indic[j* nmono+ires]=val;
                    indic[l* nmono+ires]=val;
                    
                    ires++;
                    if (degre >3) {
                        for (m=l; m<p; m++) {
                            if ( (type==4)|| (type==5))
/* without 4 input interactions or without interactions*/
                                val=presence4(3, i,j,l,m);
                            else
                                val=presence4(1, i,j,l,m);
                            indic[i* nmono+ires]=val;
                            indic[j* nmono+ires]=val;
                            indic[l* nmono+ires]=val;
                            indic[m* nmono+ires]=val;
                            ires++;
                        }
                    } /* fin (degre >3) */
                }
            } /* fin (degre >2) */
        }
    }
    
    
} /* fin generIndic */

/* ================================================== */
/* The gateway function:
 * ENTREES
 * nlhs: le nombre de sorties
 * nrhs: le nombre d'arguments d'entrée
 * prhs: les arguments d'entrée
 * SORTIE
 * plhs: les sorties
 * prhs: les arguments d'entrée
 * TEST
 * np=4; nmono=14; degre=2; indic=zeros(nmono, np, 'uint32');
 * a=generIndic(nmono, np, degre, indic)
 * ------------------------------------------------- */
void mexFunction( int nlhs, mxArray *plhs[],
        int nrhs, const mxArray *prhs[] )
{
    int  nmono, np,degre, type;
    unsigned int *indic;
    nmono=(int ) mxGetScalar(prhs[0]);
    np=(int) mxGetScalar(prhs[1]);
    degre=(int) mxGetScalar(prhs[2]);
    type=(int) mxGetScalar(prhs[3]);
    
    
    
    /* create the output matrix */
    plhs[0] = mxCreateNumericMatrix((mwSize)nmono,(mwSize)np, mxUINT32_CLASS, mxREAL);
    
    /* get a pointer to the real data in the output matrix */
    indic = (unsigned int *)mxGetPr(plhs[0]);
    
    
    generIndic( nmono, np, degre, type, indic);
    
}
/* fin mexFunction */

