| 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 | } |