Commit | Line | Data |
---|---|---|
8e92c49c BA |
1 | #include "utils.h" |
2 | #include <stdlib.h> | |
d6d71630 | 3 | #include <math.h> |
1d3c1faa BA |
4 | #include <gsl/gsl_linalg.h> |
5 | ||
b42f0f40 | 6 | // TODO: don't recompute indexes ai(...) and mi(...) when possible |
2e813ad2 | 7 | void EMGLLF_core( |
1d3c1faa | 8 | // IN parameters |
9ff729fb BA |
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 | |
8e92c49c BA |
13 | int mini, // nombre minimal d'itérations dans l'algorithme EM |
14 | int maxi, // nombre maximal d'itérations dans l'algorithme EM | |
9ff729fb BA |
15 | Real gamma, // puissance des proportions dans la pénalisation pour un Lasso adaptatif |
16 | Real lambda, // valeur du paramètre de régularisation du Lasso | |
17 | const Real* X, // régresseurs | |
18 | const Real* Y, // réponse | |
19 | Real tau, // seuil pour accepter la convergence | |
1d3c1faa | 20 | // OUT parameters (all pointers, to be modified) |
9ff729fb BA |
21 | Real* phi, // parametre de moyenne renormalisé, calculé par l'EM |
22 | Real* rho, // parametre de variance renormalisé, calculé par l'EM | |
23 | Real* pi, // parametre des proportions renormalisé, calculé par l'EM | |
d6d71630 BA |
24 | Real* llh, // (derniere) log vraisemblance associée à cet échantillon, |
25 | // pour les valeurs estimées des paramètres | |
9ff729fb | 26 | Real* S, |
8be79c46 | 27 | int* affec, |
1d3c1faa | 28 | // additional size parameters |
8e92c49c BA |
29 | int n, // nombre d'echantillons |
30 | int p, // nombre de covariables | |
31 | int m, // taille de Y (multivarié) | |
32 | int k) // nombre de composantes dans le mélange | |
1d3c1faa BA |
33 | { |
34 | //Initialize outputs | |
35 | copyArray(phiInit, phi, p*m*k); | |
36 | copyArray(rhoInit, rho, m*m*k); | |
37 | copyArray(piInit, pi, k); | |
1d3c1faa | 38 | //S is already allocated, and doesn't need to be 'zeroed' |
4cab944a | 39 | |
b42f0f40 | 40 | //Other local variables: same as in R |
9ff729fb | 41 | Real* gam = (Real*)malloc(n*k*sizeof(Real)); |
1d3c1faa | 42 | copyArray(gamInit, gam, n*k); |
b42f0f40 BA |
43 | Real* Gram2 = (Real*)malloc(p*p*k*sizeof(Real)); |
44 | Real* ps2 = (Real*)malloc(p*m*k*sizeof(Real)); | |
9ff729fb | 45 | Real* b = (Real*)malloc(k*sizeof(Real)); |
b42f0f40 BA |
46 | Real* X2 = (Real*)malloc(n*p*k*sizeof(Real)); |
47 | Real* Y2 = (Real*)malloc(n*m*k*sizeof(Real)); | |
d6d71630 | 48 | *llh = -INFINITY; |
9ff729fb | 49 | Real* pi2 = (Real*)malloc(k*sizeof(Real)); |
b42f0f40 BA |
50 | const Real EPS = 1e-15; |
51 | // Additional (not at this place, in R file) | |
52 | Real* gam2 = (Real*)malloc(k*sizeof(Real)); | |
ef67d338 | 53 | Real* sqNorm2 = (Real*)malloc(k*sizeof(Real)); |
d6d71630 | 54 | Real* detRho = (Real*)malloc(k*sizeof(Real)); |
1d3c1faa BA |
55 | gsl_matrix* matrix = gsl_matrix_alloc(m, m); |
56 | gsl_permutation* permutation = gsl_permutation_alloc(m); | |
9ff729fb BA |
57 | Real* YiRhoR = (Real*)malloc(m*sizeof(Real)); |
58 | Real* XiPhiR = (Real*)malloc(m*sizeof(Real)); | |
ef67d338 | 59 | const Real gaussConstM = pow(2.*M_PI,m/2.); |
b42f0f40 BA |
60 | Real* Phi = (Real*)malloc(p*m*k*sizeof(Real)); |
61 | Real* Rho = (Real*)malloc(m*m*k*sizeof(Real)); | |
62 | Real* Pi = (Real*)malloc(k*sizeof(Real)); | |
4cab944a | 63 | |
d6d71630 | 64 | for (int ite=0; ite<maxi; ite++) |
1d3c1faa BA |
65 | { |
66 | copyArray(phi, Phi, p*m*k); | |
67 | copyArray(rho, Rho, m*m*k); | |
68 | copyArray(pi, Pi, k); | |
4cab944a BA |
69 | |
70 | // Calculs associés a Y et X | |
71 | for (int r=0; r<k; r++) | |
1d3c1faa | 72 | { |
4cab944a | 73 | for (int mm=0; mm<m; mm++) |
1d3c1faa | 74 | { |
ef67d338 | 75 | //Y2[,mm,r] = sqrt(gam[,r]) * Y[,mm] |
4cab944a | 76 | for (int u=0; u<n; u++) |
435cb841 | 77 | Y2[ai(u,mm,r,n,m,k)] = sqrt(gam[mi(u,r,n,k)]) * Y[mi(u,mm,n,m)]; |
1d3c1faa | 78 | } |
4cab944a | 79 | for (int i=0; i<n; i++) |
1d3c1faa | 80 | { |
ef67d338 | 81 | //X2[i,,r] = sqrt(gam[i,r]) * X[i,] |
4cab944a | 82 | for (int u=0; u<p; u++) |
e39bc178 | 83 | X2[ai(i,u,r,n,p,k)] = sqrt(gam[mi(i,r,n,k)]) * X[mi(i,u,n,p)]; |
1d3c1faa | 84 | } |
4cab944a | 85 | for (int mm=0; mm<m; mm++) |
1d3c1faa | 86 | { |
ef67d338 | 87 | //ps2[,mm,r] = crossprod(X2[,,r],Y2[,mm,r]) |
4cab944a | 88 | for (int u=0; u<p; u++) |
1d3c1faa | 89 | { |
9ff729fb | 90 | Real dotProduct = 0.; |
4cab944a | 91 | for (int v=0; v<n; v++) |
46a2e676 | 92 | dotProduct += X2[ai(v,u,r,n,p,k)] * Y2[ai(v,mm,r,n,m,k)]; |
e39bc178 | 93 | ps2[ai(u,mm,r,p,m,k)] = dotProduct; |
1d3c1faa BA |
94 | } |
95 | } | |
4cab944a | 96 | for (int j=0; j<p; j++) |
1d3c1faa | 97 | { |
4cab944a | 98 | for (int s=0; s<p; s++) |
1d3c1faa | 99 | { |
ef67d338 | 100 | //Gram2[j,s,r] = crossprod(X2[,j,r], X2[,s,r]) |
9ff729fb | 101 | Real dotProduct = 0.; |
4cab944a BA |
102 | for (int u=0; u<n; u++) |
103 | dotProduct += X2[ai(u,j,r,n,p,k)] * X2[ai(u,s,r,n,p,k)]; | |
104 | Gram2[ai(j,s,r,p,p,k)] = dotProduct; | |
1d3c1faa BA |
105 | } |
106 | } | |
107 | } | |
108 | ||
109 | ///////////// | |
110 | // Etape M // | |
111 | ///////////// | |
4cab944a | 112 | |
1d3c1faa | 113 | // Pour pi |
4cab944a | 114 | for (int r=0; r<k; r++) |
1d3c1faa | 115 | { |
ef67d338 | 116 | //b[r] = sum(abs(phi[,,r])) |
9ff729fb | 117 | Real sumAbsPhi = 0.; |
4cab944a BA |
118 | for (int u=0; u<p; u++) |
119 | for (int v=0; v<m; v++) | |
120 | sumAbsPhi += fabs(phi[ai(u,v,r,p,m,k)]); | |
1d3c1faa BA |
121 | b[r] = sumAbsPhi; |
122 | } | |
ef67d338 | 123 | //gam2 = colSums(gam) |
4cab944a | 124 | for (int u=0; u<k; u++) |
1d3c1faa | 125 | { |
9ff729fb | 126 | Real sumOnColumn = 0.; |
4cab944a BA |
127 | for (int v=0; v<n; v++) |
128 | sumOnColumn += gam[mi(v,u,n,k)]; | |
1d3c1faa BA |
129 | gam2[u] = sumOnColumn; |
130 | } | |
ef67d338 | 131 | //a = sum(gam %*% log(pi)) |
9ff729fb | 132 | Real a = 0.; |
4cab944a | 133 | for (int u=0; u<n; u++) |
1d3c1faa | 134 | { |
9ff729fb | 135 | Real dotProduct = 0.; |
4cab944a BA |
136 | for (int v=0; v<k; v++) |
137 | dotProduct += gam[mi(u,v,n,k)] * log(pi[v]); | |
1d3c1faa BA |
138 | a += dotProduct; |
139 | } | |
4cab944a | 140 | |
1d3c1faa | 141 | //tant que les proportions sont negatives |
d6d71630 BA |
142 | int kk = 0, |
143 | pi2AllPositive = 0; | |
9ff729fb | 144 | Real invN = 1./n; |
1d3c1faa BA |
145 | while (!pi2AllPositive) |
146 | { | |
ef67d338 BA |
147 | //pi2 = pi + 0.1^kk * ((1/n)*gam2 - pi) |
148 | Real pow_01_kk = pow(0.1,kk); | |
4cab944a | 149 | for (int r=0; r<k; r++) |
ef67d338 BA |
150 | pi2[r] = pi[r] + pow_01_kk * (invN*gam2[r] - pi[r]); |
151 | //pi2AllPositive = all(pi2 >= 0) | |
1d3c1faa | 152 | pi2AllPositive = 1; |
4cab944a | 153 | for (int r=0; r<k; r++) |
1d3c1faa BA |
154 | { |
155 | if (pi2[r] < 0) | |
156 | { | |
157 | pi2AllPositive = 0; | |
158 | break; | |
159 | } | |
160 | } | |
161 | kk++; | |
162 | } | |
4cab944a | 163 | |
435cb841 | 164 | //sum(pi^gamma * b) |
9ff729fb | 165 | Real piPowGammaDotB = 0.; |
4cab944a | 166 | for (int v=0; v<k; v++) |
1d3c1faa | 167 | piPowGammaDotB += pow(pi[v],gamma) * b[v]; |
435cb841 | 168 | //sum(pi2^gamma * b) |
9ff729fb | 169 | Real pi2PowGammaDotB = 0.; |
4cab944a | 170 | for (int v=0; v<k; v++) |
1d3c1faa | 171 | pi2PowGammaDotB += pow(pi2[v],gamma) * b[v]; |
435cb841 BA |
172 | //sum(gam2 * log(pi2)) |
173 | Real gam2DotLogPi2 = 0.; | |
4cab944a | 174 | for (int v=0; v<k; v++) |
435cb841 BA |
175 | gam2DotLogPi2 += gam2[v] * log(pi2[v]); |
176 | ||
ef67d338 | 177 | //t(m) la plus grande valeur dans la grille O.1^k tel que ce soit décroissante ou constante |
435cb841 | 178 | while (-invN*a + lambda*piPowGammaDotB < -invN*gam2DotLogPi2 + lambda*pi2PowGammaDotB |
8e92c49c | 179 | && kk<1000) |
1d3c1faa | 180 | { |
ef67d338 BA |
181 | Real pow_01_kk = pow(0.1,kk); |
182 | //pi2 = pi + 0.1^kk * (1/n*gam2 - pi) | |
4cab944a | 183 | for (int v=0; v<k; v++) |
ef67d338 | 184 | pi2[v] = pi[v] + pow_01_kk * (invN*gam2[v] - pi[v]); |
435cb841 | 185 | //pi2 was updated, so we recompute pi2PowGammaDotB and gam2DotLogPi2 |
4cab944a BA |
186 | pi2PowGammaDotB = 0.; |
187 | for (int v=0; v<k; v++) | |
1d3c1faa | 188 | pi2PowGammaDotB += pow(pi2[v],gamma) * b[v]; |
435cb841 | 189 | gam2DotLogPi2 = 0.; |
4cab944a | 190 | for (int v=0; v<k; v++) |
435cb841 | 191 | gam2DotLogPi2 += gam2[v] * log(pi2[v]); |
1d3c1faa BA |
192 | kk++; |
193 | } | |
9ff729fb | 194 | Real t = pow(0.1,kk); |
ef67d338 | 195 | //sum(pi + t*(pi2-pi)) |
9ff729fb | 196 | Real sumPiPlusTbyDiff = 0.; |
4cab944a | 197 | for (int v=0; v<k; v++) |
1d3c1faa | 198 | sumPiPlusTbyDiff += (pi[v] + t*(pi2[v] - pi[v])); |
ef67d338 | 199 | //pi = (pi + t*(pi2-pi)) / sum(pi + t*(pi2-pi)) |
4cab944a | 200 | for (int v=0; v<k; v++) |
1d3c1faa | 201 | pi[v] = (pi[v] + t*(pi2[v] - pi[v])) / sumPiPlusTbyDiff; |
4cab944a | 202 | |
1d3c1faa | 203 | //Pour phi et rho |
4cab944a | 204 | for (int r=0; r<k; r++) |
1d3c1faa | 205 | { |
4cab944a | 206 | for (int mm=0; mm<m; mm++) |
1d3c1faa | 207 | { |
d6d71630 BA |
208 | Real ps = 0., |
209 | nY2 = 0.; | |
210 | // Compute ps, and nY2 = sum(Y2[,mm,r]^2) | |
4cab944a | 211 | for (int i=0; i<n; i++) |
1d3c1faa | 212 | { |
435cb841 | 213 | //< X2[i,,r] , phi[,mm,r] > |
ef67d338 | 214 | Real dotProduct = 0.; |
4cab944a | 215 | for (int u=0; u<p; u++) |
a2d68d1d | 216 | dotProduct += X2[ai(i,u,r,n,p,k)] * phi[ai(u,mm,r,p,m,k)]; |
d6d71630 BA |
217 | //ps = ps + Y2[i,mm,r] * sum(X2[i,,r] * phi[,mm,r]) |
218 | ps += Y2[ai(i,mm,r,n,m,k)] * dotProduct; | |
219 | nY2 += Y2[ai(i,mm,r,n,m,k)] * Y2[ai(i,mm,r,n,m,k)]; | |
1d3c1faa | 220 | } |
d6d71630 BA |
221 | //rho[mm,mm,r] = (ps+sqrt(ps^2+4*nY2*gam2[r])) / (2*nY2) |
222 | rho[ai(mm,mm,r,m,m,k)] = (ps + sqrt(ps*ps + 4*nY2 * gam2[r])) / (2*nY2); | |
1d3c1faa BA |
223 | } |
224 | } | |
435cb841 | 225 | |
4cab944a | 226 | for (int r=0; r<k; r++) |
1d3c1faa | 227 | { |
4cab944a | 228 | for (int j=0; j<p; j++) |
1d3c1faa | 229 | { |
4cab944a | 230 | for (int mm=0; mm<m; mm++) |
1d3c1faa | 231 | { |
d6d71630 | 232 | //sum(phi[-j,mm,r] * Gram2[j,-j,r]) |
435cb841 | 233 | Real phiDotGram2 = 0.; |
b42f0f40 BA |
234 | for (int u=0; u<p; u++) |
235 | { | |
236 | if (u != j) | |
435cb841 | 237 | phiDotGram2 += phi[ai(u,mm,r,p,m,k)] * Gram2[ai(j,u,r,p,p,k)]; |
b42f0f40 | 238 | } |
435cb841 | 239 | //S[j,mm,r] = -rho[mm,mm,r]*ps2[j,mm,r] + sum(phi[-j,mm,r] * Gram2[j,-j,r]) |
d6d71630 BA |
240 | S[ai(j,mm,r,p,m,k)] = -rho[ai(mm,mm,r,m,m,k)] * ps2[ai(j,mm,r,p,m,k)] |
241 | + phiDotGram2; | |
435cb841 BA |
242 | Real pirPowGamma = pow(pi[r],gamma); |
243 | if (fabs(S[ai(j,mm,r,p,m,k)]) <= n*lambda*pirPowGamma) | |
244 | phi[ai(j,mm,r,p,m,k)] = 0.; | |
245 | else if (S[ai(j,mm,r,p,m,k)] > n*lambda*pirPowGamma) | |
ef67d338 | 246 | { |
435cb841 | 247 | phi[ai(j,mm,r,p,m,k)] = (n*lambda*pirPowGamma - S[ai(j,mm,r,p,m,k)]) |
4cab944a | 248 | / Gram2[ai(j,j,r,p,p,k)]; |
ef67d338 | 249 | } |
1d3c1faa | 250 | else |
ef67d338 | 251 | { |
435cb841 | 252 | phi[ai(j,mm,r,p,m,k)] = -(n*lambda*pirPowGamma + S[ai(j,mm,r,p,m,k)]) |
4cab944a | 253 | / Gram2[ai(j,j,r,p,p,k)]; |
ef67d338 | 254 | } |
1d3c1faa BA |
255 | } |
256 | } | |
257 | } | |
4cab944a | 258 | |
1d3c1faa BA |
259 | ///////////// |
260 | // Etape E // | |
261 | ///////////// | |
4cab944a | 262 | |
d6d71630 | 263 | // Precompute det(rho[,,r]) for r in 1...k |
2e813ad2 | 264 | int signum; |
d6d71630 BA |
265 | for (int r=0; r<k; r++) |
266 | { | |
267 | for (int u=0; u<m; u++) | |
268 | { | |
269 | for (int v=0; v<m; v++) | |
270 | matrix->data[u*m+v] = rho[ai(u,v,r,m,m,k)]; | |
271 | } | |
272 | gsl_linalg_LU_decomp(matrix, permutation, &signum); | |
273 | detRho[r] = gsl_linalg_LU_det(matrix, signum); | |
274 | } | |
275 | ||
d6d71630 | 276 | Real sumLogLLH = 0.; |
4cab944a | 277 | for (int i=0; i<n; i++) |
1d3c1faa | 278 | { |
4cab944a | 279 | for (int r=0; r<k; r++) |
1d3c1faa | 280 | { |
ef67d338 | 281 | //compute Y[i,]%*%rho[,,r] |
4cab944a | 282 | for (int u=0; u<m; u++) |
1d3c1faa | 283 | { |
b42f0f40 | 284 | YiRhoR[u] = 0.; |
4cab944a | 285 | for (int v=0; v<m; v++) |
aa8df014 | 286 | YiRhoR[u] += Y[mi(i,v,n,m)] * rho[ai(v,u,r,m,m,k)]; |
1d3c1faa | 287 | } |
4cab944a | 288 | |
435cb841 | 289 | //compute X[i,]%*%phi[,,r] |
4cab944a | 290 | for (int u=0; u<m; u++) |
1d3c1faa | 291 | { |
b42f0f40 | 292 | XiPhiR[u] = 0.; |
4cab944a BA |
293 | for (int v=0; v<p; v++) |
294 | XiPhiR[u] += X[mi(i,v,n,p)] * phi[ai(v,u,r,p,m,k)]; | |
1d3c1faa | 295 | } |
4cab944a | 296 | |
ef67d338 | 297 | //compute sq norm || Y(:,i)*rho(:,:,r)-X(i,:)*phi(:,:,r) ||_2^2 |
b42f0f40 | 298 | sqNorm2[r] = 0.; |
4cab944a | 299 | for (int u=0; u<m; u++) |
ef67d338 | 300 | sqNorm2[r] += (YiRhoR[u]-XiPhiR[u]) * (YiRhoR[u]-XiPhiR[u]); |
1d3c1faa | 301 | } |
ef67d338 | 302 | |
b42f0f40 | 303 | Real sumGamI = 0.; |
4cab944a | 304 | for (int r=0; r<k; r++) |
1d3c1faa | 305 | { |
d6d71630 BA |
306 | gam[mi(i,r,n,k)] = pi[r] * exp(-.5*sqNorm2[r]) * detRho[r]; |
307 | sumGamI += gam[mi(i,r,n,k)]; | |
1d3c1faa | 308 | } |
435cb841 | 309 | |
d6d71630 BA |
310 | sumLogLLH += log(sumGamI) - log(gaussConstM); |
311 | if (sumGamI > EPS) //else: gam[i,] is already ~=0 | |
1d3c1faa | 312 | { |
d6d71630 BA |
313 | for (int r=0; r<k; r++) |
314 | gam[mi(i,r,n,k)] /= sumGamI; | |
1d3c1faa BA |
315 | } |
316 | } | |
ef67d338 BA |
317 | |
318 | //sumPen = sum(pi^gamma * b) | |
b42f0f40 | 319 | Real sumPen = 0.; |
4cab944a | 320 | for (int r=0; r<k; r++) |
1d3c1faa | 321 | sumPen += pow(pi[r],gamma) * b[r]; |
d6d71630 BA |
322 | Real last_llh = *llh; |
323 | //llh = -sumLogLLH/n + lambda*sumPen | |
324 | *llh = -invN * sumLogLLH + lambda * sumPen; | |
325 | Real dist = ite==0 ? *llh : (*llh - last_llh) / (1. + fabs(*llh)); | |
ef67d338 BA |
326 | |
327 | //Dist1 = max( abs(phi-Phi) / (1+abs(phi)) ) | |
b42f0f40 | 328 | Real Dist1 = 0.; |
4cab944a | 329 | for (int u=0; u<p; u++) |
1d3c1faa | 330 | { |
4cab944a | 331 | for (int v=0; v<m; v++) |
1d3c1faa | 332 | { |
4cab944a | 333 | for (int w=0; w<k; w++) |
1d3c1faa | 334 | { |
b42f0f40 BA |
335 | Real tmpDist = fabs(phi[ai(u,v,w,p,m,k)]-Phi[ai(u,v,w,p,m,k)]) |
336 | / (1.+fabs(phi[ai(u,v,w,p,m,k)])); | |
1d3c1faa BA |
337 | if (tmpDist > Dist1) |
338 | Dist1 = tmpDist; | |
339 | } | |
340 | } | |
341 | } | |
ef67d338 | 342 | //Dist2 = max( (abs(rho-Rho)) / (1+abs(rho)) ) |
b42f0f40 | 343 | Real Dist2 = 0.; |
4cab944a | 344 | for (int u=0; u<m; u++) |
1d3c1faa | 345 | { |
4cab944a | 346 | for (int v=0; v<m; v++) |
1d3c1faa | 347 | { |
4cab944a | 348 | for (int w=0; w<k; w++) |
1d3c1faa | 349 | { |
b42f0f40 BA |
350 | Real tmpDist = fabs(rho[ai(u,v,w,m,m,k)]-Rho[ai(u,v,w,m,m,k)]) |
351 | / (1.+fabs(rho[ai(u,v,w,m,m,k)])); | |
1d3c1faa BA |
352 | if (tmpDist > Dist2) |
353 | Dist2 = tmpDist; | |
354 | } | |
355 | } | |
356 | } | |
ef67d338 | 357 | //Dist3 = max( (abs(pi-Pi)) / (1+abs(Pi))) |
b42f0f40 | 358 | Real Dist3 = 0.; |
4cab944a | 359 | for (int u=0; u<n; u++) |
1d3c1faa | 360 | { |
4cab944a | 361 | for (int v=0; v<k; v++) |
1d3c1faa | 362 | { |
b42f0f40 | 363 | Real tmpDist = fabs(pi[v]-Pi[v]) / (1.+fabs(pi[v])); |
1d3c1faa BA |
364 | if (tmpDist > Dist3) |
365 | Dist3 = tmpDist; | |
366 | } | |
367 | } | |
368 | //dist2=max([max(Dist1),max(Dist2),max(Dist3)]); | |
d6d71630 | 369 | Real dist2 = Dist1; |
1d3c1faa BA |
370 | if (Dist2 > dist2) |
371 | dist2 = Dist2; | |
372 | if (Dist3 > dist2) | |
373 | dist2 = Dist3; | |
ef67d338 | 374 | |
d6d71630 BA |
375 | if (ite >= mini && (dist >= tau || dist2 >= sqrt(tau))) |
376 | break; | |
1d3c1faa | 377 | } |
ef67d338 | 378 | |
8be79c46 BA |
379 | //affec = apply(gam, 1, which.max) |
380 | for (int i=0; i<n; i++) | |
381 | { | |
382 | Real rowMax = 0.; | |
383 | affec[i] = 0; | |
384 | for (int j=0; j<k; j++) | |
385 | { | |
386 | if (gam[mi(i,j,n,k)] > rowMax) | |
387 | { | |
388 | affec[i] = j+1; //R indices start at 1 | |
389 | rowMax = gam[mi(i,j,n,k)]; | |
390 | } | |
391 | } | |
392 | } | |
37e11bb0 | 393 | |
1d3c1faa BA |
394 | //free memory |
395 | free(b); | |
396 | free(gam); | |
1d3c1faa BA |
397 | free(Phi); |
398 | free(Rho); | |
399 | free(Pi); | |
1d3c1faa BA |
400 | free(Gram2); |
401 | free(ps2); | |
402 | gsl_matrix_free(matrix); | |
403 | gsl_permutation_free(permutation); | |
404 | free(XiPhiR); | |
405 | free(YiRhoR); | |
406 | free(gam2); | |
407 | free(pi2); | |
408 | free(X2); | |
409 | free(Y2); | |
ef67d338 | 410 | free(sqNorm2); |
1d3c1faa | 411 | } |