Changeset 36180 for trunk/psLib

Timestamp:

Oct 3, 2013, 12:17:43 PM (13 years ago)

Author:

watersc1

Message:

Fix to mixture model code for cases where NAN values appear. This also fixes an unrelated calculation error that sometimes comes up in those cases.

File:

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

trunk/psLib/src/math/psMixtureModels.c

@@ -198 +198 @@
     psVectorInit(modes,NAN);
     for (int i = 0; i < N - Ncensored; i++) {
-      modes->data.F32[i - offsets->data.S32[i]] = tempModes->data.F32[i];
+      modes->data.F32[i + offsets->data.S32[i]] = tempModes->data.F32[i];
     }
     psFree(tempModes);
@@ -418 +418 @@
     }
     if (!P) {
-      P = psImageAlloc(N,m,PS_TYPE_F32);
+      P = psImageAlloc(m,N,PS_TYPE_F32);
     }
     psVectorInit(modes,NAN);
@@ -424 +424 @@
 
     for (int i = 0; i < N - Ncensored; i++) {
-      modes->data.F32[i - offsets->data.S32[i]] = tempModes->data.F32[i];
+      modes->data.F32[i + offsets->data.S32[i]] = tempModes->data.F32[i];
       for (int j = 0; j < dim; j++) {
-        P->data.F32[i - offsets->data.S32[i]][j] = tempP->data.F32[i][j];
+        P->data.F32[i + offsets->data.S32[i]][j] = tempP->data.F32[i][j];
       }
     }
@@ -511 +511 @@
   }
   for (i = 0; i < m; i++) {
-    pi->data.F32[i] = counts->data.F32[i] / N;
+    //pi->data.F32[i] = counts->data.F32[i] / (1.0 * N);
+    pi->data.F32[i] = 0.5;
   }
 
@@ -539 +540 @@
         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] =
@@ -545 +550 @@
 
         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);
         norm_K += P->data.F32[k][i];
       }
@@ -612 +622 @@
         }
       }
-    }    
+    }
+    //    printf("%d %f %f %f %f\n",*iterations,logL,oldlogL,means->data.F32[0][0],means->data.F32[1][0]);
+    *iterations = *iterations + 1;
   } while ((logL - oldlogL > 0)&&
            (*iterations < MAX_ITERATIONS));