From 435cb8419ebcdb624aa053f351c981133d58d6b6 Mon Sep 17 00:00:00 2001
From: Benjamin Auder <benjamin.auder@somewhere>
Date: Sun, 2 Apr 2017 11:11:35 +0200
Subject: [PATCH] Fix errors in EMGLLF.c :)

---
 pkg/src/sources/EMGLLF.c         | 55 +++++++++++++++++---------------
 test/generate_test_data/EMGLLF.R | 34 ++++++++------------
 2 files changed, 43 insertions(+), 46 deletions(-)

diff --git a/pkg/src/sources/EMGLLF.c b/pkg/src/sources/EMGLLF.c
index 7ce6ded..520d57b 100644
--- a/pkg/src/sources/EMGLLF.c
+++ b/pkg/src/sources/EMGLLF.c
@@ -77,7 +77,7 @@ void EMGLLF_core(
 			{
 				//Y2[,mm,r] = sqrt(gam[,r]) * Y[,mm]
 				for (int u=0; u<n; u++)
-					Y2[ai(u,mm,r,n,m,k)] = sqrt(gam[mi(u,r,n,k)]) * Y[mi(u,mm,m,n)];
+					Y2[ai(u,mm,r,n,m,k)] = sqrt(gam[mi(u,r,n,k)]) * Y[mi(u,mm,n,m)];
 			}
 			for (int i=0; i<n; i++)
 			{
@@ -164,33 +164,34 @@ void EMGLLF_core(
 			kk++;
 		}
 
-		//(pi.^gamma)*b
+		//sum(pi^gamma * b)
 		Real piPowGammaDotB = 0.;
 		for (int v=0; v<k; v++)
 			piPowGammaDotB += pow(pi[v],gamma) * b[v];
-		//(pi2.^gamma)*b
+		//sum(pi2^gamma * b)
 		Real pi2PowGammaDotB = 0.;
 		for (int v=0; v<k; v++)
 			pi2PowGammaDotB += pow(pi2[v],gamma) * b[v];
-		//transpose(gam2)*log(pi2)
-		Real prodGam2logPi2 = 0.;
+		//sum(gam2 * log(pi2))
+		Real gam2DotLogPi2 = 0.;
 		for (int v=0; v<k; v++)
-			prodGam2logPi2 += gam2[v] * log(pi2[v]);
+			gam2DotLogPi2 += gam2[v] * log(pi2[v]);
+
 		//t(m) la plus grande valeur dans la grille O.1^k tel que ce soit décroissante ou constante
-		while (-invN*a + lambda*piPowGammaDotB < -invN*prodGam2logPi2 + lambda*pi2PowGammaDotB
+		while (-invN*a + lambda*piPowGammaDotB < -invN*gam2DotLogPi2 + lambda*pi2PowGammaDotB
 			&& kk<1000)
 		{
 			Real pow_01_kk = pow(0.1,kk);
 			//pi2 = pi + 0.1^kk * (1/n*gam2 - pi)
 			for (int v=0; v<k; v++)
 				pi2[v] = pi[v] + pow_01_kk * (invN*gam2[v] - pi[v]);
-			//pi2 was updated, so we recompute pi2PowGammaDotB and prodGam2logPi2
+			//pi2 was updated, so we recompute pi2PowGammaDotB and gam2DotLogPi2
 			pi2PowGammaDotB = 0.;
 			for (int v=0; v<k; v++)
 				pi2PowGammaDotB += pow(pi2[v],gamma) * b[v];
-			prodGam2logPi2 = 0.;
+			gam2DotLogPi2 = 0.;
 			for (int v=0; v<k; v++)
-				prodGam2logPi2 += gam2[v] * log(pi2[v]);
+				gam2DotLogPi2 += gam2[v] * log(pi2[v]);
 			kk++;
 		}
 		Real t = pow(0.1,kk);
@@ -209,7 +210,7 @@ void EMGLLF_core(
 			{
 				for (int i=0; i<n; i++)
 				{
-					//< X2(i,:,r) , phi(:,mm,r) >
+					//< X2[i,,r] , phi[,mm,r] >
 					Real dotProduct = 0.;
 					for (int u=0; u<p; u++)
 						dotProduct += X2[ai(i,u,r,n,p,k)] * phi[ai(u,mm,r,p,m,k)];
@@ -221,16 +222,17 @@ void EMGLLF_core(
 				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(Y2[,mm,r])
-				Real sumNy2 = 0.;
+				//nY2[mm,r] = sum(Y2[,mm,r]^2)
+				Real sumY2 = 0.;
 				for (int u=0; u<n; u++)
-					sumNy2 += Y2[ai(u,mm,r,n,m,k)];
-				nY2[mi(mm,r,m,k)] = sumNy2;
+					sumY2 += Y2[ai(u,mm,r,n,m,k)] * Y2[ai(u,mm,r,n,m,k)];
+				nY2[mi(mm,r,m,k)] = sumY2;
 				//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)]);
 			}
 		}
+
 		for (int r=0; r<k; r++)
 		{
 			for (int j=0; j<p; j++)
@@ -238,25 +240,25 @@ void EMGLLF_core(
 				for (int mm=0; mm<m; mm++)
 				{
 					//sum(phi[-j,mm,r] * Gram2[j, setdiff(1:p,j),r])
-					Real dotPhiGram2 = 0.0;
+					Real phiDotGram2 = 0.;
 					for (int u=0; u<p; u++)
 					{
 						if (u != j)
-							dotPhiGram2 += phi[ai(u,mm,r,p,m,k)] * Gram2[ai(j,u,r,p,p,k)];
+							phiDotGram2 += phi[ai(u,mm,r,p,m,k)] * Gram2[ai(j,u,r,p,p,k)];
 					}
-					//S[j,mm,r] = -rho[mm,mm,r]*ps2[j,mm,r] + sum(phi[-j,mm,r] * Gram2[j, setdiff(1:p,j),r])
-					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;
-					Real pow_pir_gamma = pow(pi[r],gamma);
-					if (fabs(S[ai(j,mm,r,p,m,k)]) <= n*lambda*pow_pir_gamma)
-						phi[ai(j,mm,r,p,m,k)] = 0;
-					else if (S[ai(j,mm,r,p,m,k)] > n*lambda*pow_pir_gamma)
+					//S[j,mm,r] = -rho[mm,mm,r]*ps2[j,mm,r] + sum(phi[-j,mm,r] * Gram2[j,-j,r])
+					S[ai(j,mm,r,p,m,k)] = -rho[ai(mm,mm,r,m,m,k)] * ps2[ai(j,mm,r,p,m,k)] + phiDotGram2;
+					Real pirPowGamma = pow(pi[r],gamma);
+					if (fabs(S[ai(j,mm,r,p,m,k)]) <= n*lambda*pirPowGamma)
+						phi[ai(j,mm,r,p,m,k)] = 0.;
+					else if (S[ai(j,mm,r,p,m,k)] > n*lambda*pirPowGamma)
 					{
-						phi[ai(j,mm,r,p,m,k)] = (n*lambda*pow_pir_gamma - S[ai(j,mm,r,p,m,k)])
+						phi[ai(j,mm,r,p,m,k)] = (n*lambda*pirPowGamma - S[ai(j,mm,r,p,m,k)])
 							/ Gram2[ai(j,j,r,p,p,k)];
 					}
 					else
 					{
-						phi[ai(j,mm,r,p,m,k)] = -(n*lambda*pow_pir_gamma + S[ai(j,mm,r,p,m,k)])
+						phi[ai(j,mm,r,p,m,k)] = -(n*lambda*pirPowGamma + S[ai(j,mm,r,p,m,k)])
 							/ Gram2[ai(j,j,r,p,p,k)];
 					}
 				}
@@ -281,7 +283,7 @@ void EMGLLF_core(
 						YiRhoR[u] += Y[mi(i,v,n,m)] * rho[ai(v,u,r,m,m,k)];
 				}
 
-				//compute X(i,:)*phi(:,:,r)
+				//compute X[i,]%*%phi[,,r]
 				for (int u=0; u<m; u++)
 				{
 					XiPhiR[u] = 0.;
@@ -311,6 +313,7 @@ void EMGLLF_core(
 				sumLLF1 += Gam[mi(i,r,n,k)] / gaussConstM;
 				sumGamI += Gam[mi(i,r,n,k)];
 			}
+
 			sumLogLLF2 += log(sumLLF1);
 			for (int r=0; r<k; r++)
 			{
diff --git a/test/generate_test_data/EMGLLF.R b/test/generate_test_data/EMGLLF.R
index 37859d8..fce0a8f 100644
--- a/test/generate_test_data/EMGLLF.R
+++ b/test/generate_test_data/EMGLLF.R
@@ -50,17 +50,17 @@ EMGLLF_R = function(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,ta
           Gram2[j,s,r] = crossprod(X2[,j,r], X2[,s,r])
       }
     }
-    
+
     ##########
     #Etape M #
     ##########
     
     #pour pi
-    for (r in 1:k){
-      b[r] = sum(abs(phi[,,r]))}
+    for (r in 1:k)
+      b[r] = sum(abs(phi[,,r]))
     gam2 = colSums(gam)
     a = sum(gam %*% log(pi))
-    
+
     #tant que les props sont negatives
     kk = 0
     pi2AllPositive = FALSE
@@ -71,8 +71,7 @@ EMGLLF_R = function(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,ta
       kk = kk+1
     }
 
-#if (ite==2) browser()
-    #t[m] la plus grande valeur dans la grille O.1^k tel que ce soit décroissante ou constante
+    #t(m) la plus grande valeur dans la grille O.1^k tel que ce soit décroissante ou constante
     while( kk < 1000 && -a/n + lambda * sum(pi^gamma * b) <
 			-sum(gam2 * log(pi2))/n + lambda * sum(pi2^gamma * b) )
 		{
@@ -81,7 +80,7 @@ EMGLLF_R = function(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,ta
     }
     t = 0.1^kk
     pi = (pi + t*(pi2-pi)) / sum(pi + t*(pi2-pi))
-    
+
     #Pour phi et rho
     for (r in 1:k)
 		{
@@ -94,18 +93,17 @@ EMGLLF_R = function(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,ta
         ps[mm,r] = sum(ps1[,mm,r])
         nY2[mm,r] = sum(Y2[,mm,r]^2)
         rho[mm,mm,r] = (ps[mm,r]+sqrt(ps[mm,r]^2+4*nY2[mm,r]*gam2[r])) / (2*nY2[mm,r])
-      }
+			}
     }
+
     for (r in 1:k)
 		{
       for (j in 1:p)
 			{
         for (mm in 1:m)
 				{
-          S[j,mm,r] = -rho[mm,mm,r]*ps2[j,mm,r] + sum(phi[-j,mm,r] * Gram2[j, setdiff(1:p, j),r])
-#						(if(j>1) sum(phi[1:(j-1),mm,r] * Gram2[j,1:(j-1),r]) else 0) +
-#						(if(j<p) sum(phi[(j+1):p,mm,r] * Gram2[j,(j+1):p,r]) else 0)
-          if (abs(S[j,mm,r]) <= n*lambda*(pi[r]^gamma))
+          S[j,mm,r] = -rho[mm,mm,r]*ps2[j,mm,r] + sum(phi[-j,mm,r] * Gram2[j,-j,r])
+					if (abs(S[j,mm,r]) <= n*lambda*(pi[r]^gamma))
             phi[j,mm,r]=0
           else if(S[j,mm,r] > n*lambda*(pi[r]^gamma))
             phi[j,mm,r] = (n*lambda*(pi[r]^gamma)-S[j,mm,r]) / Gram2[j,j,r]
@@ -124,16 +122,14 @@ EMGLLF_R = function(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,ta
 		{
       #precompute sq norms to numerically adjust their values
       sqNorm2 = rep(0,k)
-      for (r in 1:k){
-        sqNorm2[r] = sum( (Y[i,]%*%rho[,,r]-X[i,]%*%phi[,,r])^2 )}
+      for (r in 1:k)
+        sqNorm2[r] = sum( (Y[i,]%*%rho[,,r]-X[i,]%*%phi[,,r])^2 )
 
-      #compute Gam(:,:)
+      #compute Gam[,]
       sumLLF1 = 0.0;
       for (r in 1:k)
 			{
-				#FIXME: numerical problems, because 0 < det(Rho[,,r] < EPS; what to do ?!
-        #       consequence: error in while() at line 77
-				Gam[i,r] = pi[r] * exp(-0.5*sqNorm2[r])* det(rho[,,r])
+				Gam[i,r] = pi[r] * exp(-0.5*sqNorm2[r]) * det(rho[,,r])
         sumLLF1 = sumLLF1 + Gam[i,r] / (2*base::pi)^(m/2)
       }
       sumLogLLF2 = sumLogLLF2 + log(sumLLF1)
@@ -146,9 +142,7 @@ EMGLLF_R = function(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,ta
 
     sumPen = sum(pi^gamma * b)
     LLF[ite] = -sumLogLLF2/n + lambda*sumPen
-
     dist = ifelse( ite == 1, LLF[ite], (LLF[ite]-LLF[ite-1]) / (1+abs(LLF[ite])) )
-
     Dist1 = max( (abs(phi-Phi)) / (1+abs(phi)) )
     Dist2 = max( (abs(rho-Rho)) / (1+abs(rho)) )
     Dist3 = max( (abs(pi-Pi)) / (1+abs(Pi)) )
-- 
2.44.0