Commit | Line | Data |
---|---|---|
1d3c1faa | 1 | #include "EMGLLF.h" |
8e92c49c BA |
2 | #include "utils.h" |
3 | #include <stdlib.h> | |
1d3c1faa BA |
4 | #include <gsl/gsl_linalg.h> |
5 | #include <omp.h> | |
1d3c1faa BA |
6 | |
7 | // TODO: comment on constructionModelesLassoMLE purpose | |
09ab3c16 | 8 | void constructionModelesLassoMLE_core( |
3ec579a0 | 9 | // IN parameters |
afa07d41 BA |
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 | |
3ec579a0 BA |
14 | int mini,// nombre minimal d'itérations dans l'algorithme EM |
15 | int maxi,// nombre maximal d'itérations dans l'algorithme EM | |
afa07d41 BA |
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 | |
3ec579a0 BA |
22 | const int* A1, // matrice des coefficients des parametres selectionnes |
23 | const int* A2, // matrice des coefficients des parametres non selectionnes | |
1d3c1faa | 24 | // OUT parameters |
afa07d41 BA |
25 | float* phi,// estimateur ainsi calculé par le Lasso |
26 | float* rho,// estimateur ainsi calculé par le Lasso | |
27 | float* pi, // estimateur ainsi calculé par le Lasso | |
28 | float* lvraisemblance, // estimateur ainsi calculé par le Lasso | |
1d3c1faa | 29 | // additional size parameters |
3ec579a0 BA |
30 | int n, // taille de l'echantillon |
31 | int p, // nombre de covariables | |
32 | int m, // taille de Y (multivarié) | |
33 | int k, // nombre de composantes | |
34 | int L) // taille de glambda | |
1d3c1faa BA |
35 | { |
36 | //preparation: phi = 0 | |
3ec579a0 | 37 | for (int u=0; u<p*m*k*L; u++) |
1d3c1faa | 38 | phi[u] = 0.0; |
3ec579a0 | 39 | |
1d3c1faa | 40 | //initiate parallel section |
3ec579a0 | 41 | int lambdaIndex; |
1d3c1faa BA |
42 | omp_set_num_threads(OMP_NUM_THREADS); |
43 | #pragma omp parallel default(shared) private(lambdaIndex) | |
44 | { | |
45 | #pragma omp for schedule(dynamic,CHUNK_SIZE) nowait | |
46 | for (lambdaIndex=0; lambdaIndex<L; lambdaIndex++) | |
47 | { | |
48 | //~ a = A1(:,1,lambdaIndex); | |
49 | //~ a(a==0) = []; | |
3ec579a0 BA |
50 | int* a = (int*)malloc(p*sizeof(int)); |
51 | int lengthA = 0; | |
52 | for (int j=0; j<p; j++) | |
1d3c1faa | 53 | { |
3ec579a0 BA |
54 | if (A1[ai(j,0,lambdaIndex,p,m+1,L)] != 0) |
55 | a[lengthA++] = A1[ai(j,0,lambdaIndex,p,m+1,L)] - 1; | |
1d3c1faa BA |
56 | } |
57 | if (lengthA == 0) | |
58 | continue; | |
3ec579a0 | 59 | |
1d3c1faa | 60 | //Xa = X(:,a) |
afa07d41 | 61 | float* Xa = (float*)malloc(n*lengthA*sizeof(float)); |
3ec579a0 | 62 | for (int i=0; i<n; i++) |
1d3c1faa | 63 | { |
3ec579a0 BA |
64 | for (int j=0; j<lengthA; j++) |
65 | Xa[mi(i,j,n,lengthA)] = X[mi(i,a[j],n,p)]; | |
1d3c1faa | 66 | } |
3ec579a0 | 67 | |
1d3c1faa | 68 | //phia = phiInit(a,:,:) |
afa07d41 | 69 | float* phia = (float*)malloc(lengthA*m*k*sizeof(float)); |
3ec579a0 | 70 | for (int j=0; j<lengthA; j++) |
1d3c1faa | 71 | { |
3ec579a0 | 72 | for (int mm=0; mm<m; mm++) |
1d3c1faa | 73 | { |
3ec579a0 BA |
74 | for (int r=0; r<k; r++) |
75 | phia[ai(j,mm,r,lengthA,m,k)] = phiInit[ai(a[j],mm,r,p,m,k)]; | |
1d3c1faa BA |
76 | } |
77 | } | |
3ec579a0 | 78 | |
1d3c1faa BA |
79 | //[phiLambda,rhoLambda,piLambda,~,~] = EMGLLF(... |
80 | // phiInit(a,:,:),rhoInit,piInit,gamInit,mini,maxi,gamma,0,X(:,a),Y,tau); | |
afa07d41 BA |
81 | float* phiLambda = (float*)malloc(lengthA*m*k*sizeof(float)); |
82 | float* rhoLambda = (float*)malloc(m*m*k*sizeof(float)); | |
83 | float* piLambda = (float*)malloc(k*sizeof(float)); | |
84 | float* LLF = (float*)malloc((maxi+1)*sizeof(float)); | |
85 | float* S = (float*)malloc(lengthA*m*k*sizeof(float)); | |
09ab3c16 | 86 | EMGLLF_core(phia,rhoInit,piInit,gamInit,mini,maxi,gamma,0.0,Xa,Y,tau, |
1d3c1faa BA |
87 | phiLambda,rhoLambda,piLambda,LLF,S, |
88 | n,lengthA,m,k); | |
89 | free(Xa); | |
90 | free(phia); | |
91 | free(LLF); | |
92 | free(S); | |
3ec579a0 | 93 | |
1d3c1faa BA |
94 | //~ for j=1:length(a) |
95 | //~ phi(a(j),:,:,lambdaIndex) = phiLambda(j,:,:); | |
96 | //~ end | |
3ec579a0 | 97 | for (int j=0; j<lengthA; j++) |
1d3c1faa | 98 | { |
3ec579a0 | 99 | for (int mm=0; mm<m; mm++) |
1d3c1faa | 100 | { |
3ec579a0 BA |
101 | for (int r=0; r<k; r++) |
102 | phi[ai4(a[j],mm,r,lambdaIndex,p,m,k,L)] = phiLambda[ai(j,mm,r,p,m,k)]; | |
1d3c1faa BA |
103 | } |
104 | } | |
105 | free(phiLambda); | |
106 | //~ rho(:,:,:,lambdaIndex) = rhoLambda; | |
3ec579a0 | 107 | for (int u=0; u<m; u++) |
1d3c1faa | 108 | { |
3ec579a0 | 109 | for (int v=0; v<m; v++) |
1d3c1faa | 110 | { |
3ec579a0 BA |
111 | for (int r=0; r<k; r++) |
112 | rho[ai4(u,v,r,lambdaIndex,m,m,k,L)] = rhoLambda[ai(u,v,r,m,m,k)]; | |
1d3c1faa BA |
113 | } |
114 | } | |
115 | free(rhoLambda); | |
116 | //~ pi(:,lambdaIndex) = piLambda; | |
3ec579a0 BA |
117 | for (int r=0; r<k; r++) |
118 | pi[mi(r,lambdaIndex,k,L)] = piLambda[r]; | |
1d3c1faa | 119 | free(piLambda); |
3ec579a0 BA |
120 | |
121 | int dimension = 0; | |
122 | int* b = (int*)malloc(m*sizeof(int)); | |
123 | for (int j=0; j<p; j++) | |
1d3c1faa BA |
124 | { |
125 | //~ b = A2(j,2:end,lambdaIndex); | |
126 | //~ b(b==0) = []; | |
3ec579a0 BA |
127 | int lengthB = 0; |
128 | for (int mm=0; mm<m; mm++) | |
1d3c1faa | 129 | { |
3ec579a0 BA |
130 | if (A2[ai(j,mm+1,lambdaIndex,p,m+1,L)] != 0) |
131 | b[lengthB++] = A2[ai(j,mm+1,lambdaIndex,p,m+1,L)] - 1; | |
1d3c1faa BA |
132 | } |
133 | //~ if length(b) > 0 | |
134 | //~ phi(A2(j,1,lambdaIndex),b,:,lambdaIndex) = 0.0; | |
135 | //~ end | |
136 | if (lengthB > 0) | |
137 | { | |
3ec579a0 | 138 | for (int mm=0; mm<lengthB; mm++) |
1d3c1faa | 139 | { |
3ec579a0 | 140 | for (int r=0; r<k; r++) |
09ab3c16 | 141 | phi[ai4( A2[ai(j,0,lambdaIndex,p,m+1,L)]-1, b[mm], r, lambdaIndex, p, m, k, L)] = 0.; |
1d3c1faa BA |
142 | } |
143 | } | |
3ec579a0 | 144 | |
1d3c1faa BA |
145 | //~ c = A1(j,2:end,lambdaIndex); |
146 | //~ c(c==0) = []; | |
147 | //~ dimension = dimension + length(c); | |
3ec579a0 | 148 | for (int mm=0; mm<m; mm++) |
1d3c1faa | 149 | { |
3ec579a0 | 150 | if (A1[ai(j,mm+1,lambdaIndex,p,m+1,L)] != 0) |
1d3c1faa BA |
151 | dimension++; |
152 | } | |
153 | } | |
154 | free(b); | |
3ec579a0 | 155 | |
1d3c1faa | 156 | int signum; |
afa07d41 BA |
157 | float* densite = (float*)calloc(L*n,sizeof(float)); |
158 | float sumLogDensit = 0.0; | |
1d3c1faa BA |
159 | gsl_matrix* matrix = gsl_matrix_alloc(m, m); |
160 | gsl_permutation* permutation = gsl_permutation_alloc(m); | |
afa07d41 BA |
161 | float* YiRhoR = (float*)malloc(m*sizeof(float)); |
162 | float* XiPhiR = (float*)malloc(m*sizeof(float)); | |
3ec579a0 | 163 | for (int i=0; i<n; i++) |
1d3c1faa BA |
164 | { |
165 | //~ for r=1:k | |
166 | //~ delta = Y(i,:)*rho(:,:,r,lambdaIndex) - (X(i,a)*(phi(a,:,r,lambdaIndex))); | |
167 | //~ densite(i,lambdaIndex) = densite(i,lambdaIndex) +... | |
168 | //~ pi(r,lambdaIndex)*det(rho(:,:,r,lambdaIndex))/(sqrt(2*PI))^m*exp(-dot(delta,delta)/2.0); | |
169 | //~ end | |
3ec579a0 | 170 | for (int r=0; r<k; r++) |
1d3c1faa BA |
171 | { |
172 | //compute det(rho(:,:,r,lambdaIndex)) [TODO: avoid re-computations] | |
3ec579a0 | 173 | for (int u=0; u<m; u++) |
1d3c1faa | 174 | { |
3ec579a0 BA |
175 | for (int v=0; v<m; v++) |
176 | matrix->data[u*m+v] = rho[ai4(u,v,r,lambdaIndex,m,m,k,L)]; | |
1d3c1faa BA |
177 | } |
178 | gsl_linalg_LU_decomp(matrix, permutation, &signum); | |
afa07d41 | 179 | float detRhoR = gsl_linalg_LU_det(matrix, signum); |
3ec579a0 | 180 | |
1d3c1faa | 181 | //compute Y(i,:)*rho(:,:,r,lambdaIndex) |
3ec579a0 | 182 | for (int u=0; u<m; u++) |
1d3c1faa BA |
183 | { |
184 | YiRhoR[u] = 0.0; | |
3ec579a0 BA |
185 | for (int v=0; v<m; v++) |
186 | YiRhoR[u] += Y[mi(i,v,n,m)] * rho[ai4(v,u,r,lambdaIndex,m,m,k,L)]; | |
1d3c1faa | 187 | } |
3ec579a0 | 188 | |
1d3c1faa | 189 | //compute X(i,a)*phi(a,:,r,lambdaIndex) |
3ec579a0 | 190 | for (int u=0; u<m; u++) |
1d3c1faa BA |
191 | { |
192 | XiPhiR[u] = 0.0; | |
3ec579a0 BA |
193 | for (int v=0; v<lengthA; v++) |
194 | XiPhiR[u] += X[mi(i,a[v],n,p)] * phi[ai4(a[v],u,r,lambdaIndex,p,m,k,L)]; | |
1d3c1faa | 195 | } |
3ec579a0 BA |
196 | // On peut remplacer X par Xa dans ce dernier calcul, mais je ne sais pas si c'est intéressant ... |
197 | ||
1d3c1faa | 198 | // compute dotProduct < delta . delta > |
afa07d41 | 199 | float dotProduct = 0.0; |
3ec579a0 | 200 | for (int u=0; u<m; u++) |
1d3c1faa | 201 | dotProduct += (YiRhoR[u]-XiPhiR[u]) * (YiRhoR[u]-XiPhiR[u]); |
3ec579a0 BA |
202 | |
203 | densite[mi(lambdaIndex,i,L,n)] += (pi[mi(r,lambdaIndex,k,L)]*detRhoR/pow(sqrt(2.0*M_PI),m))*exp(-dotProduct/2.0); | |
204 | } | |
1d3c1faa BA |
205 | sumLogDensit += log(densite[lambdaIndex*n+i]); |
206 | } | |
3ec579a0 BA |
207 | lvraisemblance[mi(lambdaIndex,0,L,2)] = sumLogDensit; |
208 | lvraisemblance[mi(lambdaIndex,1,L,2)] = (dimension+m+1)*k-1; | |
209 | ||
1d3c1faa BA |
210 | free(a); |
211 | free(YiRhoR); | |
212 | free(XiPhiR); | |
213 | free(densite); | |
214 | gsl_matrix_free(matrix); | |
215 | gsl_permutation_free(permutation); | |
216 | } | |
217 | } | |
218 | } |