6 // Main job on raw inputs (after transformation from mxArray)
9 const double* phiInit
, // parametre initial de moyenne renormalisé
10 const double* rhoInit
, // parametre initial de variance renormalisé
11 const double* piInit
,// parametre initial des proportions
12 const double* 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 double gamma
, // valeur de gamma : puissance des proportions dans la pénalisation pour un Lasso adaptatif
16 const double* glambda
, // valeur des paramètres de régularisation du Lasso
17 const double* X
,// régresseurs
18 const double* Y
,// réponse
19 double seuil
, // seuil pour prendre en compte une variable
20 double 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 double* Rho
,// estimateur ainsi calculé par le Lasso
25 double* Pi
,// estimateur ainsi calculé par le Lasso
26 // additional size parameters
27 int n
,// taille de lambdaIndex'echantillon
28 int p
,// nombre de covariables
29 int m
,// taille de Y (multivarié)
30 int k
,// nombre de composantes
31 int 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 double* phi
= (double*)malloc(p
*m
*k
*sizeof(double));
54 double* rho
= (double*)malloc(m
*m
*k
*sizeof(double));
55 double* pi
= (double*)malloc(k
*sizeof(double));
56 double* LLF
= (double*)malloc(maxi
*sizeof(double));
57 double* S
= (double*)malloc(p
*m
*k
*sizeof(double));
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 int* selectedVariables
= (int*)calloc(p
*m
,sizeof(int));
66 int* discardedVariables
= (int*)calloc(p
*m
,sizeof(int));
67 int atLeastOneSelectedVariable
= 0;
68 for (int j
=0; j
<p
; j
++)
72 for (int jj
=0; jj
<m
; jj
++)
75 for (int r
=0; r
<k
; r
++)
77 if (fabs(phi
[ai(j
,jj
,r
,p
,m
,k
)]) > maxPhi
)
78 maxPhi
= fabs(phi
[ai(j
,jj
,r
,p
,m
,k
)]);
82 selectedVariables
[mi(j
,cpt
,p
,m
)] = jj
+1;
83 atLeastOneSelectedVariable
= 1;
88 discardedVariables
[mi(j
,cpt2
,p
,m
)] = jj
+1;
95 if (atLeastOneSelectedVariable
)
97 int* vec
= (int*)malloc(p
*sizeof(int));
99 for (int j
=0; j
<p
; j
++)
101 if (selectedVariables
[mi(j
,0,p
,m
)] != 0)
106 for (int j
=0; j
<p
; j
++)
107 A1
[ai(j
,0,lambdaIndex
,p
,m
+1,L
)] = (j
< vecSize
? vec
[j
]+1 : 0);
108 for (int j
=0; j
<vecSize
; j
++)
110 for (int jj
=1; jj
<=m
; jj
++)
111 A1
[ai(j
,jj
,lambdaIndex
,p
,m
+1,L
)] = selectedVariables
[mi(vec
[j
],jj
-1,p
,m
)];
114 for (int j
=0; j
<p
; j
++)
115 A2
[ai(j
,0,lambdaIndex
,p
,m
+1,L
)] = j
+1;
116 for (int j
=0; j
<p
; j
++)
118 for (int jj
=1; jj
<=m
; jj
++)
119 A2
[ai(j
,jj
,lambdaIndex
,p
,m
+1,L
)] = discardedVariables
[mi(j
,jj
-1,p
,m
)];
122 for (int j
=0; j
<m
; j
++)
124 for (int jj
=0; jj
<m
; jj
++)
126 for (int r
=0; r
<k
; r
++)
127 Rho
[ai4(j
,jj
,r
,lambdaIndex
,m
,m
,k
,L
)] = rho
[ai(j
,jj
,r
,m
,m
,k
)];
131 for (int r
=0; r
<k
; r
++)
132 Pi
[mi(r
,lambdaIndex
,k
,L
)] = pi
[r
];
136 free(selectedVariables
);
137 free(discardedVariables
);