Changeset 40729 for trunk/psLib

Timestamp:

May 10, 2019, 9:54:48 AM (7 years ago)

Author:

eugene

Message:

do not let logL drive to -inf (limit probability to >= 1-e6)

File:

• trunk/psLib/src/math/psMixtureModels.c (modified) (9 diffs)

trunk/psLib/src/math/psMixtureModels.c

@@ -42 +42 @@
 
 
-// #define TESTING
+# define TESTING 0
+# define TESTING_VERBOSE 0
 /////////////////////////////////////////////////////////////////////////////////
 // Internal functions                                                          //
@@ -363 +364 @@
     }
 
+# 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 +559 @@
         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 +646 @@
       }
     }
-    //    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 +688 @@
   int iterations;
   double V,V0;
-#ifdef TESTING
-  int i,j;
-#endif
+
   if (!modes) {
     modes = psVectorAlloc(N,PS_TYPE_F32);
@@ -680 +699 @@
       // 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 +714 @@
       // 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 +733 @@
     // 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 +748 @@
     // 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]);