C-part of tests ready [TODO: R part]
[valse.git] / src / sources / selectiontotale.c
1 #include <stdlib.h>
2 #include <omp.h>
3 #include <math.h>
4 #include "EMGLLF.h"
5 #include "utils.h"
6
7 // Main job on raw inputs (after transformation from mxArray)
8 void selectiontotale_core(
9 // IN parameters
10 const float* phiInit, // parametre initial de moyenne renormalisé
11 const float* rhoInit, // parametre initial de variance renormalisé
12 const float* piInit,// parametre initial des proportions
13 const float* gamInit, // paramètre initial des probabilités a posteriori de chaque échantillon
14 int mini, // nombre minimal d'itérations dans lambdaIndex'algorithme EM
15 int maxi, // nombre maximal d'itérations dans lambdaIndex'algorithme EM
16 float gamma, // valeur de gamma : puissance des proportions dans la pénalisation pour un Lasso adaptatif
17 const float* glambda, // valeur des paramètres de régularisation du Lasso
18 const float* X,// régresseurs
19 const float* Y,// réponse
20 float seuil, // seuil pour prendre en compte une variable
21 float tau, // seuil pour accepter la convergence
22 // OUT parameters (all pointers, to be modified)
23 int* A1, // matrice des coefficients des parametres selectionnes
24 int* A2, // matrice des coefficients des parametres non selectionnes
25 float* Rho,// estimateur ainsi calculé par le Lasso
26 float* Pi,// estimateur ainsi calculé par le Lasso
27 // additional size parameters
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
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;
44
45 //initiate parallel section
46 int lambdaIndex;
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
54 float* phi = (float*)malloc(p*m*k*sizeof(float));
55 float* rho = (float*)malloc(m*m*k*sizeof(float));
56 float* pi = (float*)malloc(k*sizeof(float));
57 float* LLF = (float*)malloc(maxi*sizeof(float));
58 float* S = (float*)malloc(p*m*k*sizeof(float));
59 EMGLLF_core(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,glambda[lambdaIndex],X,Y,tau,
60 phi,rho,pi,LLF,S,
61 n,p,m,k);
62 free(LLF);
63 free(S);
64
65 //Si un des coefficients est supérieur au seuil, on garde cette variable
66 int* selectedVariables = (int*)calloc(p*m,sizeof(int));
67 int* discardedVariables = (int*)calloc(p*m,sizeof(int));
68 int atLeastOneSelectedVariable = 0;
69 for (int j=0; j<p; j++)
70 {
71 int cpt = 0;
72 int cpt2 = 0;
73 for (int jj=0; jj<m; jj++)
74 {
75 float maxPhi = 0.0;
76 for (int r=0; r<k; r++)
77 {
78 if (fabs(phi[ai(j,jj,r,p,m,k)]) > maxPhi)
79 maxPhi = fabs(phi[ai(j,jj,r,p,m,k)]);
80 }
81 if (maxPhi > seuil)
82 {
83 selectedVariables[mi(j,cpt,p,m)] = jj+1;
84 atLeastOneSelectedVariable = 1;
85 cpt++;
86 }
87 else
88 {
89 discardedVariables[mi(j,cpt2,p,m)] = jj+1;
90 cpt2++;
91 }
92 }
93 }
94 free(phi);
95
96 if (atLeastOneSelectedVariable)
97 {
98 int* vec = (int*)malloc(p*sizeof(int));
99 int vecSize = 0;
100 for (int j=0; j<p; j++)
101 {
102 if (selectedVariables[mi(j,0,p,m)] != 0)
103 vec[vecSize++] = j;
104 }
105
106 //A1
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++)
110 {
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)];
113 }
114 //A2
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++)
118 {
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)];
121 }
122 //Rho
123 for (int j=0; j<m; j++)
124 {
125 for (int jj=0; jj<m; jj++)
126 {
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)];
129 }
130 }
131 //Pi
132 for (int r=0; r<k; r++)
133 Pi[mi(r,lambdaIndex,k,L)] = pi[r];
134 free(vec);
135 }
136
137 free(selectedVariables);
138 free(discardedVariables);
139 free(rho);
140 free(pi);
141 }
142 }
143 }