From: Benjamin Auder <benjamin.auder@somewhere>
Date: Mon, 20 Mar 2017 03:18:01 +0000 (+0100)
Subject: remove also NY21 as in R code
X-Git-Url: https://git.auder.net/variants/current/doc/scripts/%7B%7B%20asset%28%27mixstore/%7B%7B%20pkg.url%20%7D%7D?a=commitdiff_plain;h=37e11bb06399b30c3a26d5e21a80a921898053d6;p=valse.git

remove also NY21 as in R code
---

diff --git a/pkg/src/sources/EMGLLF.c b/pkg/src/sources/EMGLLF.c
index e41fe3c..a028919 100644
--- a/pkg/src/sources/EMGLLF.c
+++ b/pkg/src/sources/EMGLLF.c
@@ -54,7 +54,6 @@ void EMGLLF_core(
 	const Real EPS = 1e-15;
 	// Additional (not at this place, in R file)
 	Real* gam2 = (Real*)malloc(k*sizeof(Real));
-	Real* nY21 = (Real*)malloc(n*m*k*sizeof(Real));
 	Real* sqNorm2 = (Real*)malloc(k*sizeof(Real));
 	gsl_matrix* matrix = gsl_matrix_alloc(m, m);
 	gsl_permutation* permutation = gsl_permutation_alloc(m);
@@ -216,18 +215,17 @@ void EMGLLF_core(
 						dotProduct += X2[ai(i,u,r,n,p,k)] * phi[ai(u,mm,r,p,m,k)];
 					//ps1[i,mm,r] = Y2[i,mm,r] * sum(X2[i,,r] * phi[,mm,r])
 					ps1[ai(i,mm,r,n,m,k)] = Y2[ai(i,mm,r,n,m,k)] * dotProduct;
-					nY21[ai(i,mm,r,n,m,k)] = Y2[ai(i,mm,r,n,m,k)] * Y2[ai(i,mm,r,n,m,k)];
 				}
 				//ps[mm,r] = sum(ps1[,mm,r])
 				Real sumPs1 = 0.;
 				for (int u=0; u<n; u++)
 					sumPs1 += ps1[ai(u,mm,r,n,m,k)];
 				ps[mi(mm,r,m,k)] = sumPs1;
-				//nY2[mm,r] = sum(nY21[,mm,r])
-				Real sumNy21 = 0.;
+				//nY2[mm,r] = sum(nY2[,mm,r])
+				Real sumNy2 = 0.;
 				for (int u=0; u<n; u++)
-					sumNy21 += nY21[ai(u,mm,r,n,m,k)];
-				nY2[mi(mm,r,m,k)] = sumNy21;
+					sumNy2 += nY2[ai(u,mm,r,n,m,k)];
+				nY2[mi(mm,r,m,k)] = sumNy2;
 				//rho[mm,mm,r] = (ps[mm,r]+sqrt(ps[mm,r]^2+4*nY2[mm,r]*(gam2[r]))) / (2*nY2[mm,r])
 				rho[ai(mm,mm,r,m,m,k)] = ( ps[mi(mm,r,m,k)] + sqrt( ps[mi(mm,r,m,k)]*ps[mi(mm,r,m,k)]
 					+ 4*nY2[mi(mm,r,m,k)] * gam2[r] ) ) / (2*nY2[mi(mm,r,m,k)]);
@@ -380,6 +378,8 @@ void EMGLLF_core(
 		ite++;
 	}
 
+	//TODO: affec = apply(gam, 1,which.max) à traduire, et retourner affec aussi.
+
 	//free memory
 	free(b);
 	free(gam);
@@ -390,7 +390,6 @@ void EMGLLF_core(
 	free(ps);
 	free(nY2);
 	free(ps1);
-	free(nY21);
 	free(Gram2);
 	free(ps2);
 	gsl_matrix_free(matrix);