Changeset 35785 for trunk/psLib

Timestamp:

Jul 9, 2013, 6:36:20 PM (13 years ago)

Author:

watersc1

Message:

Fixed psMM errors in handling NAN values. I've added a wrapper function that censors these values before passing the valie data into the prior version of the functions. This is necessary as a mask won't work (there are matrix operations I don't want to bother sorting masking for).

Location:

trunk/psLib/src/math

Files:

• psMixtureModels.c (modified) (4 diffs)
• psMixtureModels.h (modified) (3 diffs)

trunk/psLib/src/math/psMixtureModels.c

@@ -112 +112 @@
 }
 
+psImage *psMMCensorData(psImage *in,
+                        long *Ncensored,
+                        psVector *offsets) {
+  // Allocate offset vector
+  offsets = psVectorRecycle(offsets,in->numRows,PS_TYPE_S32);
+
+  *Ncensored = 0;
+
+  for (int i = 0; i < in->numRows; i++) {
+    int isCensored = 0;
+    for (int j = 0; j < in->numCols; j++) {
+      if (!isfinite(in->data.F32[i][j])) {
+        isCensored = 1;
+      }
+    }
+    if (isCensored) {
+      *Ncensored = *Ncensored + 1;
+    }
+    offsets->data.S32[i] = *Ncensored;
+  }
+  psImage *out;
+  if (*Ncensored == 0) {
+    out = psMemIncrRefCounter(in);
+  }
+  else {
+    out = psImageAlloc(in->numCols,in->numRows - *Ncensored,PS_TYPE_F32);
+    for (int i = 0; i < out->numRows; i++) {
+      for (int j = 0; j < out->numCols; j++) {
+        out->data.F32[i][j] = in->data.F32[i + offsets->data.S32[i]][j];
+      }
+    }
+  }
+  return(out);
+}
 
 /////////////////////////////////////////////////////////////////////////////////
@@ -117 +151 @@
 /////////////////////////////////////////////////////////////////////////////////
 
-
 #define MAX_ITERATIONS 100
 #define THRESHOLD      1e-3
 
-
-bool psMMkmeans(psImage *D,
-                int dim,
-                long N,
-                psVector *modes,
-                psImage *means,
-                int m,
-                int *iterations,
-                double *V) {
+bool psMMkmeans(psImage *D,      //< Input data to fit.                                          
+                int dim,         //< dimension of the data                                       
+                long N,          //< number of data points                                       
+                psVector *modes, //< Output vector of which mode each data point belongs to      
+                psImage *means,  //< Matrix containing mode mean values for each dimension       
+                int m,           //< Number of modes to fit                                      
+                int *iterations, //< Number of iterations to perform                             
+                double *V) {     //< Deviation value.  Pseudo chi^2                           
+  PS_ASSERT_IMAGE_NON_NULL(D,false);
+  PS_ASSERT_IMAGE_TYPE(D,PS_TYPE_F32,false);
+  // We need to scan the input for invalid d
+  long Ncensored;
+  psVector *offsets = psVectorAlloc(N,PS_TYPE_S32);
+  psImage  *Dcensored = psMMCensorData(D,&Ncensored,offsets);
+
+  psVector *tempModes = psVectorAlloc(N - Ncensored,PS_TYPE_F32);
+  bool success = psMMkmeansUncensored(Dcensored,
+                                      dim,
+                                      N - Ncensored,
+                                      tempModes,
+                                      means,
+                                      m,
+                                      iterations,
+                                      V);
+  if (Ncensored == 0) { // We can directly return the modes and means calculated
+    if (modes) {
+      psFree(modes);
+    }
+    modes = psMemIncrRefCounter(tempModes);
+  }
+  else {
+    if (!modes) {
+      modes = psVectorAlloc(N,PS_TYPE_F32);
+    }
+
+    psVectorInit(modes,NAN);
+    for (int i = 0; i < N - Ncensored; i++) {
+      modes->data.F32[i - offsets->data.S32[i]] = tempModes->data.F32[i];
+    }
+  }
+  psFree(offsets);
+  psFree(Dcensored);
+  psFree(tempModes);
+  return(success);
+}
+
+
+bool psMMkmeansUncensored(psImage *D,       //< Input data to fit.  
+                          int dim,          //< dimension of the data
+                          long N,           //< number of data points
+                          psVector *modes,  //< Output vector of which mode each data point belongs to
+                          psImage *means,   //< Matrix containing mode mean values for each dimension
+                          int m,            //< Number of modes to fit
+                          int *iterations,  //< Number of iterations to perform
+                          double *V) {      //< Deviation value.  Pseudo chi^2                           
   // Assertions
   PS_ASSERT_IMAGE_NON_NULL(D,false);
@@ -135 +214 @@
   //  PS_ASSERT_IMAGE_SIZE(D,N,dim,false);
 
+  
   // Allocate things we can allocate here if they're not already allocated
   if (!modes) {
-    modes = psVectorAlloc(m,PS_TYPE_F32);
+    modes = psVectorAlloc(N,PS_TYPE_F32);
   }
   if (!means) {
@@ -277 +357 @@
 }
 
-bool psMMGMM(psImage *D,
-             int dim,
-             long N,
-             psVector *modes,
-             psImage *means,
-             psArray *sigma,
-             psVector *pi,
-             psImage *P,
-             int m,
-             int *iterations,
-             double *V) {
+bool psMMGMM(psImage *D,           //< Input data to fit
+             int dim,              //< Dimension of each data vector
+             long N,               //< Number of data vectors
+             psVector *modes,      //< Output vector assigning each data point to a mode
+             psImage *means,       //< Output matrix containing vectors of the means for each mode
+             psArray *sigma,       //< Output array containing covariance matrices for each mode
+             psVector *pi,         //< Output vector contining population fractions for each mode
+             psImage *P,           //< Output matrix containing the probabilities a data point is in each mode
+             int m,                //< Number of modes to fit
+             int *iterations,      //< Number of iterations performed
+             double *V) {          //< Log-likelihood of final solution.
+
+  PS_ASSERT_IMAGE_NON_NULL(D,false);
+  PS_ASSERT_IMAGE_TYPE(D,PS_TYPE_F32,false);
+  // We need to scan the input for invalid d
+  long Ncensored;
+  psVector *offsets = psVectorAlloc(N,PS_TYPE_F32);
+  psImage  *Dcensored = psMMCensorData(D,&Ncensored,offsets);
+
+  psVector *tempModes = psVectorAlloc(N - Ncensored,PS_TYPE_F32);
+  psImage  *tempP     = psImageAlloc(m,N - Ncensored,PS_TYPE_F32);
+  bool success = psMMGMMUncensored(Dcensored,
+                                   dim,
+                                   N - Ncensored,
+                                   tempModes,
+                                   means,
+                                   sigma,
+                                   pi,
+                                   tempP,
+                                   m,
+                                   iterations,
+                                   V);
+  if (Ncensored == 0) { // We can directly return the modes and means calculated
+    if (modes) {
+      psFree(modes);
+    }
+    if (P) {
+      psFree(P);
+    }
+    modes = psMemIncrRefCounter(tempModes);
+    P     = psMemIncrRefCounter(tempP);
+  }
+  else {
+    if (!modes) {
+      modes = psVectorAlloc(N,PS_TYPE_F32);
+    }
+    if (!P) {
+      P = psImageAlloc(N,m,PS_TYPE_F32);
+    }
+    psVectorInit(modes,NAN);
+    psImageInit(P,NAN);
+
+    for (int i = 0; i < N - Ncensored; 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];
+      }
+    }
+  }
+  psFree(offsets);
+  psFree(Dcensored);
+  psFree(tempModes);
+  psFree(tempP);
+  return(success);
+}
+
+
+bool psMMGMMUncensored(psImage *D,           //< Input data to fit
+                       int dim,              //< Dimension of each data vector
+                       long N,               //< Number of data vectors
+                       psVector *modes,      //< Output vector assigning each data point to a mode
+                       psImage *means,       //< Output matrix containing vectors of the means for each mode
+                       psArray *sigma,       //< Output array containing covariance matrices for each mode
+                       psVector *pi,         //< Output vector contining population fractions for each mode
+                       psImage *P,           //< Output matrix containing the probabilities a data point is in each mode
+                       int m,                //< Number of modes to fit
+                       int *iterations,      //< Number of iterations performed
+                       double *V) {          //< Log-likelihood of final solution.
   // Assertions
   PS_ASSERT_IMAGE_NON_NULL(D,false);

trunk/psLib/src/math/psMixtureModels.h

@@ -2 +2 @@
 #define PS_MIXTUREMODELS_H
 
+// These are the recommended functions to call.  They in turn call code
+// that scans the input data matrix/image and censored data points with
+// non-finite values.  The censored data is then passed to the next set
+// of functions.
 
 bool psMMkmeans(psImage *D,
@@ -24 +28 @@
              double *V);
 
+// These functions assume that the input data matrix/image is cleaned of
+// non-finite values.  They should work identically as the ones above on
+// data that is well defined.
+bool psMMkmeansUncensored(psImage *D,
+                          int dim,
+                          long N,
+                          psVector *modes,
+                          psImage *means,
+                          int m,
+                          int *iterations,
+                          double *V);
+
+bool psMMGMMUncensored(psImage *D,
+                       int dim,
+                       long N,
+                       psVector *modes,
+                       psImage *means,
+                       psArray *sigma,
+                       psVector *pi,
+                       psImage *P,
+                       int m,
+                       int *iterations,
+                       double *V);
+
+// This function does the classification in five stages:
+//  1) pass data to K-means code with m modes for initial values.
+//  2) pass data+K-mean values to GMM code with m modes for best values.
+//  3) repeat K-means pass with m-1 modes.
+//  4) repeat GMM pass with m-1 modes.
+//  5) calculate Pvalue that the model with m modes is better than that with m-1
 bool psMMClass(psImage *D,
                int dim,
@@ -34 +68 @@
                int m,
                double *Pval);
-  
+
+
+// These functions are one dimensional equivalent wrappers of the above functions.
+// The psVector data is inserted into the appropriate sized psImage structure to
+// avoid having to do so in each program calling this library.  Since this is the
+// only special case (2-dim is a matrix, n-dim is matrix as well), no other wrappers
+// are needed.
 bool psMM1Dkmeans(psVector *D,
                   long N,