fix arrays reading in C, add type Real for double or float
[valse.git] / src / sources / selectiontotale.c
CommitLineData
8e92c49c 1#include <stdlib.h>
1d3c1faa 2#include <omp.h>
aa8df014 3#include <math.h>
8e92c49c
BA
4#include "EMGLLF.h"
5#include "utils.h"
1d3c1faa
BA
6
7// Main job on raw inputs (after transformation from mxArray)
09ab3c16 8void selectiontotale_core(
8e92c49c 9 // IN parameters
9ff729fb
BA
10 const Real* phiInit, // parametre initial de moyenne renormalisé
11 const Real* rhoInit, // parametre initial de variance renormalisé
12 const Real* piInit,// parametre initial des proportions
13 const Real* gamInit, // paramètre initial des probabilités a posteriori de chaque échantillon
8e92c49c
BA
14 int mini, // nombre minimal d'itérations dans lambdaIndex'algorithme EM
15 int maxi, // nombre maximal d'itérations dans lambdaIndex'algorithme EM
9ff729fb
BA
16 Real gamma, // valeur de gamma : puissance des proportions dans la pénalisation pour un Lasso adaptatif
17 const Real* glambda, // valeur des paramètres de régularisation du Lasso
18 const Real* X,// régresseurs
19 const Real* Y,// réponse
20 Real seuil, // seuil pour prendre en compte une variable
21 Real tau, // seuil pour accepter la convergence
1d3c1faa 22 // OUT parameters (all pointers, to be modified)
8e92c49c
BA
23 int* A1, // matrice des coefficients des parametres selectionnes
24 int* A2, // matrice des coefficients des parametres non selectionnes
9ff729fb
BA
25 Real* Rho,// estimateur ainsi calculé par le Lasso
26 Real* Pi,// estimateur ainsi calculé par le Lasso
1d3c1faa 27 // additional size parameters
8e92c49c
BA
28 int n,// taille de lambdaIndex'echantillon
29 int p,// nombre de covariables
30 int m,// taille de Y (multivarié)
31 int k,// nombre de composantes
32 int L) // taille de glambda
1d3c1faa
BA
33{
34 // Fill outputs with zeros: they might not be assigned
35 for (int u=0; u<p*(m+1)*L; u++)
36 {
37 A1[u] = 0;
38 A2[u] = 0;
39 }
40 for (int u=0; u<m*m*k*L; u++)
41 Rho[u] = 0.0;
42 for (int u=0; u<k*L; u++)
43 Pi[u] = 0.0;
8e92c49c 44
1d3c1faa 45 //initiate parallel section
8e92c49c 46 int lambdaIndex;
1d3c1faa
BA
47 omp_set_num_threads(OMP_NUM_THREADS);
48 #pragma omp parallel default(shared) private(lambdaIndex)
49 {
50 #pragma omp for schedule(dynamic,CHUNK_SIZE) nowait
51 for (lambdaIndex=0; lambdaIndex<L; lambdaIndex++)
52 {
53 //allocate output variables
9ff729fb
BA
54 Real* phi = (Real*)malloc(p*m*k*sizeof(Real));
55 Real* rho = (Real*)malloc(m*m*k*sizeof(Real));
56 Real* pi = (Real*)malloc(k*sizeof(Real));
57 Real* LLF = (Real*)malloc(maxi*sizeof(Real));
58 Real* S = (Real*)malloc(p*m*k*sizeof(Real));
aa8df014 59 EMGLLF_core(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,glambda[lambdaIndex],X,Y,tau,
1d3c1faa
BA
60 phi,rho,pi,LLF,S,
61 n,p,m,k);
62 free(LLF);
63 free(S);
8e92c49c 64
1d3c1faa 65 //Si un des coefficients est supérieur au seuil, on garde cette variable
8e92c49c
BA
66 int* selectedVariables = (int*)calloc(p*m,sizeof(int));
67 int* discardedVariables = (int*)calloc(p*m,sizeof(int));
1d3c1faa 68 int atLeastOneSelectedVariable = 0;
8e92c49c 69 for (int j=0; j<p; j++)
1d3c1faa 70 {
8e92c49c
BA
71 int cpt = 0;
72 int cpt2 = 0;
73 for (int jj=0; jj<m; jj++)
1d3c1faa 74 {
9ff729fb 75 Real maxPhi = 0.0;
8e92c49c 76 for (int r=0; r<k; r++)
1d3c1faa 77 {
8e92c49c
BA
78 if (fabs(phi[ai(j,jj,r,p,m,k)]) > maxPhi)
79 maxPhi = fabs(phi[ai(j,jj,r,p,m,k)]);
1d3c1faa
BA
80 }
81 if (maxPhi > seuil)
82 {
8e92c49c 83 selectedVariables[mi(j,cpt,p,m)] = jj+1;
1d3c1faa
BA
84 atLeastOneSelectedVariable = 1;
85 cpt++;
86 }
87 else
88 {
8e92c49c 89 discardedVariables[mi(j,cpt2,p,m)] = jj+1;
1d3c1faa
BA
90 cpt2++;
91 }
92 }
93 }
94 free(phi);
8e92c49c 95
1d3c1faa
BA
96 if (atLeastOneSelectedVariable)
97 {
8e92c49c
BA
98 int* vec = (int*)malloc(p*sizeof(int));
99 int vecSize = 0;
100 for (int j=0; j<p; j++)
1d3c1faa 101 {
8e92c49c 102 if (selectedVariables[mi(j,0,p,m)] != 0)
1d3c1faa
BA
103 vec[vecSize++] = j;
104 }
105
106 //A1
8e92c49c
BA
107 for (int j=0; j<p; j++)
108 A1[ai(j,0,lambdaIndex,p,m+1,L)] = (j < vecSize ? vec[j]+1 : 0);
109 for (int j=0; j<vecSize; j++)
1d3c1faa 110 {
8e92c49c
BA
111 for (int jj=1; jj<=m; jj++)
112 A1[ai(j,jj,lambdaIndex,p,m+1,L)] = selectedVariables[mi(vec[j],jj-1,p,m)];
1d3c1faa
BA
113 }
114 //A2
8e92c49c
BA
115 for (int j=0; j<p; j++)
116 A2[ai(j,0,lambdaIndex,p,m+1,L)] = j+1;
117 for (int j=0; j<p; j++)
1d3c1faa 118 {
8e92c49c
BA
119 for (int jj=1; jj<=m; jj++)
120 A2[ai(j,jj,lambdaIndex,p,m+1,L)] = discardedVariables[mi(j,jj-1,p,m)];
1d3c1faa
BA
121 }
122 //Rho
8e92c49c 123 for (int j=0; j<m; j++)
1d3c1faa 124 {
8e92c49c 125 for (int jj=0; jj<m; jj++)
1d3c1faa 126 {
8e92c49c
BA
127 for (int r=0; r<k; r++)
128 Rho[ai4(j,jj,r,lambdaIndex,m,m,k,L)] = rho[ai(j,jj,r,m,m,k)];
1d3c1faa
BA
129 }
130 }
131 //Pi
8e92c49c
BA
132 for (int r=0; r<k; r++)
133 Pi[mi(r,lambdaIndex,k,L)] = pi[r];
1d3c1faa
BA
134 free(vec);
135 }
8e92c49c 136
1d3c1faa
BA
137 free(selectedVariables);
138 free(discardedVariables);
139 free(rho);
140 free(pi);
141 }
142 }
143}