Commit | Line | Data |
---|---|---|
8e92c49c | 1 | #include <stdlib.h> |
1d3c1faa | 2 | #include <omp.h> |
8e92c49c BA |
3 | #include <gsl/gsl_linalg.h> |
4 | #include "EMGrank.h" | |
5 | #include "utils.h" | |
1d3c1faa BA |
6 | |
7 | // TODO: comment on constructionModelesLassoRank purpose | |
09ab3c16 | 8 | void constructionModelesLassoRank_core( |
3ec579a0 | 9 | // IN parameters |
afa07d41 BA |
10 | const float* Pi,// parametre initial des proportions |
11 | const float* Rho, // parametre initial de variance renormalisé | |
3ec579a0 BA |
12 | int mini, // nombre minimal d'itérations dans l'algorithme EM |
13 | int maxi, // nombre maximal d'itérations dans l'algorithme EM | |
afa07d41 BA |
14 | const float* X,// régresseurs |
15 | const float* Y,// réponse | |
16 | float tau, // seuil pour accepter la convergence | |
3ec579a0 BA |
17 | const int* A1, // matrice des coefficients des parametres selectionnes |
18 | int rangmin, //rang minimum autorisé | |
19 | int rangmax, //rang maximum autorisé | |
1d3c1faa | 20 | // OUT parameters (all pointers, to be modified) |
afa07d41 BA |
21 | float* phi,// estimateur ainsi calculé par le Lasso |
22 | float* lvraisemblance,// estimateur ainsi calculé par le Lasso | |
3ec579a0 BA |
23 | // additional size parameters |
24 | int n,// taille de l'echantillon | |
25 | int p,// nombre de covariables | |
26 | int m,// taille de Y (multivarié) | |
27 | int k,// nombre de composantes | |
28 | int L)// taille de glambda | |
1d3c1faa BA |
29 | { |
30 | //On cherche les rangs possiblement intéressants | |
3ec579a0 BA |
31 | int deltaRank = rangmax-rangmin+1; |
32 | int Size = (int)pow(deltaRank,k); | |
33 | int* Rank = (int*)malloc(Size*k*sizeof(int)); | |
34 | for (int r=0; r<k; r++) | |
35 | { | |
1d3c1faa BA |
36 | //On veut le tableau de toutes les combinaisons de rangs possibles |
37 | //Dans la première colonne : on répète (rangmax-rangmin)^(k-1) chaque chiffre : ca remplit la colonne | |
38 | //Dans la deuxieme : on répète (rangmax-rangmin)^(k-2) chaque chiffre, et on fait ca (rangmax-rangmin)^2 fois | |
39 | //... | |
40 | //Dans la dernière, on répète chaque chiffre une fois, et on fait ca (rangmin-rangmax)^(k-1) fois. | |
3ec579a0 BA |
41 | int indexInRank = 0; |
42 | int value = 0; | |
1d3c1faa BA |
43 | while (indexInRank < Size) |
44 | { | |
3ec579a0 BA |
45 | for (int u=0; u<pow(deltaRank,k-r-1); u++) |
46 | Rank[mi(indexInRank++,r,Size,k)] = rangmin + value; | |
1d3c1faa BA |
47 | value = (value+1) % deltaRank; |
48 | } | |
49 | } | |
3ec579a0 | 50 | |
1d3c1faa | 51 | //Initialize phi to zero, because unactive variables won't be assigned |
3ec579a0 | 52 | for (int i=0; i<p*m*k*L*Size; i++) |
1d3c1faa | 53 | phi[i] = 0.0; |
3ec579a0 | 54 | |
1d3c1faa | 55 | //initiate parallel section |
3ec579a0 | 56 | int lambdaIndex; |
1d3c1faa BA |
57 | omp_set_num_threads(OMP_NUM_THREADS); |
58 | #pragma omp parallel default(shared) private(lambdaIndex) | |
59 | { | |
60 | #pragma omp for schedule(dynamic,CHUNK_SIZE) nowait | |
61 | for (lambdaIndex=0; lambdaIndex<L; lambdaIndex++) | |
62 | { | |
63 | //On ne garde que les colonnes actives : active sera l'ensemble des variables informatives | |
3ec579a0 BA |
64 | int* active = (int*)malloc(p*sizeof(int)); |
65 | int longueurActive = 0; | |
66 | for (int j=0; j<p; j++) | |
1d3c1faa | 67 | { |
3ec579a0 BA |
68 | if (A1[mi(j,lambdaIndex,p,L)] != 0) |
69 | active[longueurActive++] = A1[mi(j,lambdaIndex,p,L)] - 1; | |
1d3c1faa | 70 | } |
3ec579a0 | 71 | |
1d3c1faa BA |
72 | if (longueurActive == 0) |
73 | continue; | |
3ec579a0 | 74 | |
1d3c1faa | 75 | //from now on, longueurActive > 0 |
afa07d41 BA |
76 | float* phiLambda = (float*)malloc(longueurActive*m*k*sizeof(float)); |
77 | float LLF; | |
3ec579a0 | 78 | for (int j=0; j<Size; j++) |
1d3c1faa BA |
79 | { |
80 | //[phiLambda,LLF] = EMGrank(Pi(:,lambdaIndex),Rho(:,:,:,lambdaIndex),mini,maxi,X(:,active),Y,tau,Rank(j,:)); | |
3ec579a0 BA |
81 | int* rank = (int*)malloc(k*sizeof(int)); |
82 | for (int r=0; r<k; r++) | |
83 | rank[r] = Rank[mi(j,r,Size,k)]; | |
afa07d41 | 84 | float* Xactive = (float*)malloc(n*longueurActive*sizeof(float)); |
3ec579a0 | 85 | for (int i=0; i<n; i++) |
1d3c1faa | 86 | { |
3ec579a0 BA |
87 | for (int jj=0; jj<longueurActive; jj++) |
88 | Xactive[mi(i,jj,n,longueurActive)] = X[mi(i,active[jj],n,p)]; | |
1d3c1faa | 89 | } |
afa07d41 | 90 | float* PiLambda = (float*)malloc(k*sizeof(float)); |
3ec579a0 BA |
91 | for (int r=0; r<k; r++) |
92 | PiLambda[r] = Pi[mi(r,lambdaIndex,k,L)]; | |
afa07d41 | 93 | float* RhoLambda = (float*)malloc(m*m*k*sizeof(float)); |
3ec579a0 | 94 | for (int u=0; u<m; u++) |
1d3c1faa | 95 | { |
3ec579a0 | 96 | for (int v=0; v<m; v++) |
1d3c1faa | 97 | { |
3ec579a0 | 98 | for (int r=0; r<k; r++) |
09ab3c16 | 99 | RhoLambda[ai(u,v,r,m,m,k)] = Rho[ai4(u,v,r,lambdaIndex,m,m,k,L)]; |
1d3c1faa BA |
100 | } |
101 | } | |
09ab3c16 | 102 | EMGrank_core(PiLambda,RhoLambda,mini,maxi,Xactive,Y,tau,rank, |
1d3c1faa BA |
103 | phiLambda,&LLF, |
104 | n,longueurActive,m,k); | |
105 | free(rank); | |
106 | free(Xactive); | |
107 | free(PiLambda); | |
108 | free(RhoLambda); | |
109 | //lvraisemblance((lambdaIndex-1)*Size+j,:) = [LLF, dot(Rank(j,:), length(active)-Rank(j,:)+m)]; | |
3ec579a0 | 110 | lvraisemblance[mi(lambdaIndex*Size+j,0,L*Size,2)] = LLF; |
1d3c1faa | 111 | //dot(Rank(j,:), length(active)-Rank(j,:)+m) |
afa07d41 | 112 | float dotProduct = 0.0; |
3ec579a0 BA |
113 | for (int r=0; r<k; r++) |
114 | dotProduct += Rank[mi(j,r,Size,k)] * (longueurActive-Rank[mi(j,r,Size,k)]+m); | |
115 | lvraisemblance[mi(lambdaIndex*Size+j,1,Size*L,2)] = dotProduct; | |
1d3c1faa | 116 | //phi(active,:,:,(lambdaIndex-1)*Size+j) = phiLambda; |
3ec579a0 | 117 | for (int jj=0; jj<longueurActive; jj++) |
1d3c1faa | 118 | { |
3ec579a0 | 119 | for (int mm=0; mm<m; mm++) |
1d3c1faa | 120 | { |
3ec579a0 BA |
121 | for (int r=0; r<k; r++) |
122 | phi[ai5(active[jj],mm,r,lambdaIndex,j,p,m,k,L,Size)] = phiLambda[jj*m*k+mm*k+r]; | |
1d3c1faa BA |
123 | } |
124 | } | |
125 | } | |
126 | free(active); | |
127 | free(phiLambda); | |
128 | } | |
129 | } | |
130 | free(Rank); | |
131 | } |