Commit | Line | Data |
---|---|---|
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 | 8 | void selectiontotale_core( |
8e92c49c | 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 | |
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 | |
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 | |
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 | |
afa07d41 BA |
25 | float* Rho,// estimateur ainsi calculé par le Lasso |
26 | float* 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 | |
afa07d41 BA |
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)); | |
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 | { |
afa07d41 | 75 | float 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 | } |