1 #include "selectiontotale.h"
4 #include "omp_num_threads.h"
6 // Main job on raw inputs (after transformation from mxArray)
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
33 // Fill outputs with zeros: they might not be assigned
34 for (int u
=0; u
<p
*(m
+1)*L
; u
++)
39 for (int u
=0; u
<m
*m
*k
*L
; u
++)
41 for (int u
=0; u
<k
*L
; u
++)
44 //initiate parallel section
46 omp_set_num_threads(OMP_NUM_THREADS
);
47 #pragma omp parallel default(shared) private(lambdaIndex)
49 #pragma omp for schedule(dynamic,CHUNK_SIZE) nowait
50 for (lambdaIndex
=0; lambdaIndex
<L
; lambdaIndex
++)
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
,
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
++)
72 for (mwSize jj
=0; jj
<m
; jj
++)
75 for (mwSize r
=0; r
<k
; r
++)
77 if (fabs(phi
[j
*m
*k
+jj
*k
+r
]) > maxPhi
)
78 maxPhi
= fabs(phi
[j
*m
*k
+jj
*k
+r
]);
82 selectedVariables
[j
*m
+cpt
] = jj
+1;
83 atLeastOneSelectedVariable
= 1;
88 discardedVariables
[j
*m
+cpt2
] = jj
+1;
95 if (atLeastOneSelectedVariable
)
97 Int
* vec
= (Int
*)malloc(p
*sizeof(Int
));
99 for (mwSize j
=0; j
<p
; j
++)
101 if (selectedVariables
[j
*m
+0] != 0)
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
++)
110 for (mwSize jj
=1; jj
<=m
; jj
++)
111 A1
[j
*(m
+1)*L
+jj
*L
+lambdaIndex
] = selectedVariables
[vec
[j
]*m
+jj
-1];
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
++)
118 for (mwSize jj
=1; jj
<=m
; jj
++)
119 A2
[j
*(m
+1)*L
+jj
*L
+lambdaIndex
] = discardedVariables
[j
*m
+jj
-1];
122 for (mwSize j
=0; j
<m
; j
++)
124 for (mwSize jj
=0; jj
<m
; jj
++)
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
];
131 for (mwSize r
=0; r
<k
; r
++)
132 Pi
[r
*L
+lambdaIndex
] = pi
[r
];
136 free(selectedVariables
);
137 free(discardedVariables
);