Changeset 36680 for branches/eam_branches/ps2-tc3-20130727/psLib/src/math/psMixtureModels.c

Timestamp:

Apr 21, 2014, 5:42:34 AM (12 years ago)

Author:

eugene

Message:

merge changes from trunk

Location:

branches/eam_branches/ps2-tc3-20130727

Files:

• . (modified) (1 prop)
• psLib/src/math/psMixtureModels.c (modified) (12 diffs, 1 prop)

branches/eam_branches/ps2-tc3-20130727/psLib/src/math/psMixtureModels.c

@@ -120 +120 @@
   *Ncensored = 0;
 
+  int k = 0;
   for (int i = 0; i < in->numRows; i++) {
     int isCensored = 0;
@@ -130 +131 @@
       *Ncensored = *Ncensored + 1;
     }
-    offsets->data.S32[i] = *Ncensored;
-  }
+    else {
+      offsets->data.S32[k] = i - k;
+      k++;
+    }
+    //    offsets->data.S32[i] = isCensored;
+  }
+
   psImage *out;
   if (*Ncensored == 0) {
@@ -198 +204 @@
     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);
@@ -320 +326 @@
         }
       }
-      i = modes->data.F32[k];
-      counts->data.F32[i] += 1.0;
+      if (modes->data.F32[k] != -1) {
+        i = modes->data.F32[k];
+        counts->data.F32[i] += 1.0;
+      }
     }
 
@@ -335 +343 @@
     for (k = 0; k < N; k++) {
       i = modes->data.F32[k];
-      for (j = 0; j < dim; j++) {
-        means->data.F32[i][j] += D->data.F32[k][j];
+      if (i != -1) {
+        for (j = 0; j < dim; j++) {
+          means->data.F32[i][j] += D->data.F32[k][j];
+        }
       }
     }
@@ -403 +413 @@
     success = psMMGMMUncensored(Dcensored,
                                 dim,
-                                N,
+                                N - Ncensored,
                                 tempModes,
                                 means,
@@ -418 +428 @@
     }
     if (!P) {
-      P = psImageAlloc(N,m,PS_TYPE_F32);
+      P = psImageAlloc(m,N,PS_TYPE_F32);
     }
     psVectorInit(modes,NAN);
@@ -424 +434 @@
 
     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 +521 @@
   }
   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 +550 @@
         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 +560 @@
 
         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 +632 @@
         }
       }
-    }    
+    }
+    //    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));