X-Git-Url: https://git.auder.net/?a=blobdiff_plain;f=pkg%2Fsrc%2Fsources%2FEMGLLF.c;h=d9f92c049745c1304b1852df9f509c18a841a308;hb=23b9fb13bc6e82d7ca43bfb83aa85b6cd69c52c0;hp=a57a37962ca2000f09117d2a362f1b9e1d9582b5;hpb=2e813ad23c796bbed3d5ba685b8fa002bdc6689d;p=valse.git

diff --git a/pkg/src/sources/EMGLLF.c b/pkg/src/sources/EMGLLF.c
index a57a379..d9f92c0 100644
--- a/pkg/src/sources/EMGLLF.c
+++ b/pkg/src/sources/EMGLLF.c
@@ -39,6 +39,7 @@ void EMGLLF_core(
 
 	//Other local variables: same as in R
 	Real* gam = (Real*)malloc(n*k*sizeof(Real));
+	Real* logGam = (Real*)malloc(k*sizeof(Real));
 	copyArray(gamInit, gam, n*k);
 	Real* Gram2 = (Real*)malloc(p*p*k*sizeof(Real));
 	Real* ps2 = (Real*)malloc(p*m*k*sizeof(Real));
@@ -47,7 +48,6 @@ void EMGLLF_core(
 	Real* Y2 = (Real*)malloc(n*m*k*sizeof(Real));
 	*llh = -INFINITY;
 	Real* pi2 = (Real*)malloc(k*sizeof(Real));
-	const Real EPS = 1e-15;
 	// Additional (not at this place, in R file)
 	Real* gam2 = (Real*)malloc(k*sizeof(Real));
 	Real* sqNorm2 = (Real*)malloc(k*sizeof(Real));
@@ -61,7 +61,7 @@ void EMGLLF_core(
 	Real* Rho = (Real*)malloc(m*m*k*sizeof(Real));
 	Real* Pi = (Real*)malloc(k*sizeof(Real));
 
-	for (int ite=0; ite<maxi; ite++)
+	for (int ite=1; ite<=maxi; ite++)
 	{
 		copyArray(phi, Phi, p*m*k);
 		copyArray(rho, Rho, m*m*k);
@@ -300,19 +300,26 @@ void EMGLLF_core(
 					sqNorm2[r] += (YiRhoR[u]-XiPhiR[u]) * (YiRhoR[u]-XiPhiR[u]);
 			}
 
-			Real sumGamI = 0.;
+			// Update gam[,]; use log to avoid numerical problems
+			Real maxLogGam = -INFINITY;
 			for (int r=0; r<k; r++)
 			{
-				gam[mi(i,r,n,k)] = pi[r] * exp(-.5*sqNorm2[r]) * detRho[r];
-				sumGamI += gam[mi(i,r,n,k)];
+				logGam[r] = log(pi[r]) - .5 * sqNorm2[r] + log(detRho[r]);
+				if (maxLogGam < logGam[r])
+					maxLogGam = logGam[r];
 			}
-
-			sumLogLLH += log(sumGamI) - log(gaussConstM);
-			if (sumGamI > EPS) //else: gam[i,] is already ~=0
+			Real norm_fact = 0.;
+			for (int r=0; r<k; r++)
 			{
-				for (int r=0; r<k; r++)
-					gam[mi(i,r,n,k)] /= sumGamI;
+				logGam[r] = logGam[r] - maxLogGam; //adjust without changing proportions
+				gam[mi(i,r,n,k)] = exp(logGam[r]); //gam[i, ] <- exp(logGam)
+				norm_fact += gam[mi(i,r,n,k)]; //norm_fact <- sum(gam[i, ])
 			}
+			// gam[i, ] <- gam[i, ] / norm_fact
+			for (int r=0; r<k; r++)
+				gam[mi(i,r,n,k)] /= norm_fact;
+
+			sumLogLLH += log(norm_fact) - log(gaussConstM);
 		}
 
 		//sumPen = sum(pi^gamma * b)
@@ -320,9 +327,9 @@ void EMGLLF_core(
 		for (int r=0; r<k; r++)
 			sumPen += pow(pi[r],gamma) * b[r];
 		Real last_llh = *llh;
-		//llh = -sumLogLLH/n + lambda*sumPen
-		*llh = -invN * sumLogLLH + lambda * sumPen;
-		Real dist = ite==0 ? *llh : (*llh - last_llh) / (1. + fabs(*llh));
+		//llh = -sumLogLLH/n #+ lambda*sumPen
+		*llh = -invN * sumLogLLH; //+ lambda * sumPen;
+		Real dist = ( ite==1 ? *llh : (*llh - last_llh) / (1. + fabs(*llh)) );
 
 		//Dist1 = max( abs(phi-Phi) / (1+abs(phi)) )
 		Real Dist1 = 0.;
@@ -394,11 +401,13 @@ void EMGLLF_core(
 	//free memory
 	free(b);
 	free(gam);
+	free(logGam);
 	free(Phi);
 	free(Rho);
 	free(Pi);
 	free(Gram2);
 	free(ps2);
+	free(detRho);
 	gsl_matrix_free(matrix);
 	gsl_permutation_free(permutation);
 	free(XiPhiR);
@@ -408,4 +417,4 @@ void EMGLLF_core(
 	free(X2);
 	free(Y2);
 	free(sqNorm2);
-}
+}