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