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