Index: trunk/psLib/src/math/psMixtureModels.c
===================================================================
--- trunk/psLib/src/math/psMixtureModels.c	(revision 40551)
+++ trunk/psLib/src/math/psMixtureModels.c	(revision 40729)
@@ -42,5 +42,6 @@
 
 
-// #define TESTING
+# define TESTING 0
+# define TESTING_VERBOSE 0
 /////////////////////////////////////////////////////////////////////////////////
 // Internal functions                                                          //
@@ -363,7 +364,15 @@
     }
 
+# if (TESTING)
+    fprintf (stderr, "kmeans : %d %f %f %f\n",*iterations,diff,means->data.F32[0][0],means->data.F32[1][0]);
+# endif
+    
     *iterations += 1;
   } while ((*iterations < MAX_ITERATIONS)&&
 	   (diff > THRESHOLD));
+
+# if (TESTING)
+  fprintf (stderr, "psMMkmeans diff: %f\n", diff);
+# endif
 
   *V = diff;
@@ -550,19 +559,23 @@
 	for (j = 0; j < dim; j++) { // VV
 	  M->data.F32[j] = means->data.F32[i][j];
-	  //	  printf("%f %f %f\n",
-	  //		 X->data.F32[j],M->data.F32[j],
-	  //		 S->data.F32[j][j]);
-	}
-
-	
-	P->data.F32[k][i] =
-	  pi->data.F32[i] * psMMNDGaussian(X,M,S);
+#if (TESTING_VERBOSE)
+	  	  printf("%f %f %f\n", X->data.F32[j],M->data.F32[j], S->data.F32[j][j]);
+#endif
+	}
+
+	// if log-likelihood becomes too small, we risk numerical overflow 
+	float Prob = psMMNDGaussian(X,M,S);
+	if (Prob < 1e-6) { Prob = 1e-6; }
+
+	P->data.F32[k][i] = pi->data.F32[i] * Prob;
 
 	logL += log(P->data.F32[k][i]);
-	//	printf("%d %d %f %f %f %f\n",
-	//	       k,i,
-	//	       pi->data.F32[i],
-	//	       psMMNDGaussian(X,M,S),
-	//	       P->data.F32[k][i],logL);
+#if (TESTING_VERBOSE)
+	printf("%d %d %f %f %f %f\n",
+		       k,i,
+		       pi->data.F32[i],
+		       psMMNDGaussian(X,M,S),
+		       P->data.F32[k][i],logL);
+#endif
 	norm_K += P->data.F32[k][i];
       }
@@ -633,8 +646,16 @@
       }
     }
-    //    printf("%d %f %f %f %f\n",*iterations,logL,oldlogL,means->data.F32[0][0],means->data.F32[1][0]);
+
+# if (TESTING)
+    fprintf (stderr, "GMM: %d %f %f %f %f\n",*iterations,logL,oldlogL,means->data.F32[0][0],means->data.F32[1][0]);
+# endif
+
     *iterations = *iterations + 1;
   } while ((logL - oldlogL > 0)&&
 	   (*iterations < MAX_ITERATIONS));
+
+# if (TESTING)
+  fprintf (stderr, "psMMGMM logL: %f\n", logL);
+# endif
 
   psFree(counts);
@@ -667,7 +688,5 @@
   int iterations;
   double V,V0;
-#ifdef TESTING
-  int i,j;
-#endif
+
   if (!modes) {
     modes = psVectorAlloc(N,PS_TYPE_F32);
@@ -680,7 +699,8 @@
       // Error calculating k-means
     }
-#ifdef TESTING
-    for (i = 0 ; i < m -1; i++) {
-      for (j = 0; j < dim; j++) {
+
+#if (TESTING)
+    for (int i = 0 ; i < m -1; i++) {
+      for (int j = 0; j < dim; j++) {
 	fprintf(stderr,"KMEANS POST: (%d/%d) (%d/%d) %g\n",
 		i,m - 1,j,dim,means->data.F32[i][j]);
@@ -694,7 +714,8 @@
       // Error calculating GMM
     }
-#ifdef TESTING
-    for (i = 0 ; i < m - 1; i++) {
-      for (j = 0; j < dim; j++) {
+
+#if (TESTING)
+    for (int i = 0 ; i < m - 1; i++) {
+      for (int j = 0; j < dim; j++) {
 	fprintf(stderr,"GMM POST: (%d/%d) (%d/%d) %g\n",
 		i,m - 1,j,dim,means->data.F32[i][j]);
@@ -712,7 +733,8 @@
     // Error calculating k-means
   }
-#ifdef TESTING
-  for (i = 0 ; i < m; i++) {
-    for (j = 0; j < dim; j++) {
+
+#if (TESTING)
+  for (int i = 0 ; i < m; i++) {
+    for (int j = 0; j < dim; j++) {
       fprintf(stderr,"KMEANS POST: (%d/%d) (%d/%d) %g\n",
 	      i,m,j,dim,means->data.F32[i][j]);
@@ -726,7 +748,8 @@
     // Error calculating GMM
   }
-#ifdef TESTING
-  for (i = 0 ; i < m; i++) {
-    for (j = 0; j < dim; j++) {
+
+#if (TESTING)
+  for (int i = 0 ; i < m; i++) {
+    for (int j = 0; j < dim; j++) {
       fprintf(stderr,"GMM POST: (%d/%d) (%d/%d) %g\n",
 	      i,m,j,dim,means->data.F32[i][j]);