vfem-engine.c

数据挖掘方面的源码

#include "vfem-engine.h"
#include <math.h>
#include <stdlib.h>

extern float gNeededDelta;
extern float gSigmaSquare;
extern int gNormalApprox;

IterationStatsPtr IterationStatsInitial(VoidAListPtr initialClusters) {
   IterationStatsPtr is;
   int i;

   is = MNewPtr(sizeof(IterationStats));
   is->centroids = VLNew();
   for(i = 0 ; i < VALLength(initialClusters) ; i++) {
      VALAppend(is->centroids, ExampleClone(VALIndex(initialClusters, i)));
   }

   is->cMax = VLNew();
   for(i = 0 ; i < VALLength(initialClusters) ; i++) {
      VALAppend(is->cMax, ExampleClone(VALIndex(initialClusters, i)));
   }

   is->cMaxAssignment = VLNew();
   for(i = 0 ; i < VALLength(initialClusters) ; i++) {
      VALAppend(is->cMaxAssignment,
                 ExampleClone(VALIndex(initialClusters, i)));
   }

   is->cMin = VLNew();
   for(i = 0 ; i < VALLength(initialClusters) ; i++) {
      VALAppend(is->cMin, ExampleClone(VALIndex(initialClusters, i)));
   }

   is->cMinAssignment = VLNew();
   for(i = 0 ; i < VALLength(initialClusters) ; i++) {
      VALAppend(is->cMinAssignment,
                 ExampleClone(VALIndex(initialClusters, i)));
   }
   is->lastBound = MNewPtr(sizeof(float) * VALLength(is->centroids));
   is->lastAssignmentBound = MNewPtr(sizeof(float) * VALLength(is->centroids));

   is->w = MNewPtr(sizeof(double) * VALLength(is->centroids));
   is->wPlus = MNewPtr(sizeof(double) * VALLength(is->centroids));
   is->wPlusSquare = MNewPtr(sizeof(double) * VALLength(is->centroids));
   is->wMinus = MNewPtr(sizeof(double) * VALLength(is->centroids));


   is->wx = MNewPtr(sizeof(double *) * VALLength(is->centroids));
   is->wxPlus = MNewPtr(sizeof(double *) * VALLength(is->centroids));
   is->wxMinus = MNewPtr(sizeof(double *) * VALLength(is->centroids));

   is->deltaPlus = MNewPtr(sizeof(double *) * VALLength(is->centroids));
   is->deltaMinus = MNewPtr(sizeof(double *) * VALLength(is->centroids));

   for(i = 0 ; i < VALLength(initialClusters) ; i++) {
      is->wx[i] = MNewPtr(sizeof(double) * 
          ExampleGetNumAttributes(VALIndex(is->centroids, 0)));
      is->wxPlus[i] = MNewPtr(sizeof(double) * 
          ExampleGetNumAttributes(VALIndex(is->centroids, 0)));
      is->wxMinus[i] = MNewPtr(sizeof(double) * 
          ExampleGetNumAttributes(VALIndex(is->centroids, 0)));
      is->deltaPlus[i] = MNewPtr(sizeof(double) * 
          ExampleGetNumAttributes(VALIndex(is->centroids, 0)));
      is->deltaMinus[i] = MNewPtr(sizeof(double) * 
          ExampleGetNumAttributes(VALIndex(is->centroids, 0)));
   }

   IterationStatsClearCounts(is);

   /* we always have a perfect bound for an initial is */
   is->foundBound = 1;
   is->isFirstIteration = 1;

   return is;
}

void IterationStatsClearCounts(IterationStatsPtr is) {
   int i, j;

   is->n = 0;

   for(i = 0 ; i < VALLength(is->centroids) ; i++) {
      is->lastBound[i] = 0;
      is->lastAssignmentBound[i] = 0;
      is->w[i] = 0;
      is->wPlus[i] = 0;
      is->wPlusSquare[i] = 0;
      is->wMinus[i] = 0;

      for(j = 0 ; j < ExampleGetNumAttributes(VALIndex(is->centroids, 0)) ; j++) {
         is->wx[i][j] = 0;
         is->wxPlus[i][j] = 0;
         is->wxMinus[i][j] = 0;
         is->deltaPlus[i][j] = 0;
         is->deltaMinus[i][j] = 0;
      }
   }

   is->maxEkd = 0;
   is->possibleIDConverge = 0;
   is->guarenteeIDConverge = 0;
   is->wouldEMConverge = 0;
   is->convergeVFEM = 0;
}

void IterationStatsFree(IterationStatsPtr is) {
   int numCentroids = VALLength(is->centroids);
   int i;
   ExamplePtr e;

   for(i = 0 ; i < numCentroids ; i++) {
      e = VALIndex(is->centroids, i);
      ExampleFree(e);
      e = VALIndex(is->cMin, i);
      ExampleFree(e);
      e = VALIndex(is->cMinAssignment, i);
      ExampleFree(e);
      e = VALIndex(is->cMax, i);
      ExampleFree(e);
      e = VALIndex(is->cMaxAssignment, i);
      ExampleFree(e);
   }
   VALFree(is->centroids);
   VALFree(is->cMinAssignment);
   VALFree(is->cMaxAssignment);
   VALFree(is->cMin);
   VALFree(is->cMax);

   MFreePtr(is->lastBound);
   MFreePtr(is->lastAssignmentBound);

   for(i = 0 ; i < VALLength(is->cMin) ; i++) {
      MFreePtr(is->wx[i]);
      MFreePtr(is->wxPlus[i]);
      MFreePtr(is->wxMinus[i]);
      MFreePtr(is->deltaMinus[i]);
      MFreePtr(is->deltaPlus[i]);
   }

   MFreePtr(is->w);
   MFreePtr(is->wx);
   MFreePtr(is->wPlus);
   MFreePtr(is->wPlusSquare);
   MFreePtr(is->wMinus);
   MFreePtr(is->wxPlus);
   MFreePtr(is->wxMinus);
   MFreePtr(is->deltaMinus);
   MFreePtr(is->deltaPlus);

   MFreePtr(is);
}

/* this is so terrible */
static void _MakeNewCentroids(VoidAListPtr centers, VoidAListPtr newCenters,
    VoidAListPtr newAssignment, double **wx, double *w, double *wPlusSquare,
    double *wMinus, long n, ExampleSpecPtr es, float delta, int useSampleError,
       int boundUpper, int pedroBound, double **deltaPlus, 
              double **deltaMinus) {
   double epsilon, *thisWx;
   ExamplePtr newCenter, newCenterAssignment = 0, e;
   int i,j;
   double assignmentError;

   /* free any old centers from the list */
   while(VALLength(newCenters) > 0) {
      i = VALLength(newCenters) - 1;
      e = VALRemove(newCenters, i);
      ExampleFree(e);
   }

   if(useSampleError) {
      while(VALLength(newAssignment) > 0) {
         i = VALLength(newAssignment) - 1;
         e = VALRemove(newAssignment, i);
         ExampleFree(e);
      }
   }

   // This is my bound
   //if(useBound) {
   //   epsilon = ((double)n) * log(2.0 / delta);
   //   epsilon /= 2.0;
   //   epsilon = sqrt(epsilon);
   //} else {
   //   epsilon = 0;
   //}

   if(gMessageLevel > 2) {
      printf("Making new centroids.\n");
      fflush(stdout);
   }

   /* make the list of new centers */
   for(i = 0 ; i < VALLength(centers) ; i++) {
      if(useSampleError) {
         newCenterAssignment = ExampleNew(es);
         ExampleSetClass(newCenterAssignment, 0);
      }

      newCenter = ExampleNew(es);
      ExampleSetClass(newCenter, 0);

      if(useSampleError) {
         epsilon = wPlusSquare[i] / (wMinus[i] * wMinus[i]);
         if(epsilon > 1) {
            epsilon = 1;
         } else if(epsilon < (1.0 / n)) {
            epsilon = 1.0 / n;
         }

         epsilon *= log(2.0 / delta);
         epsilon /= 2.0;
         epsilon = sqrt(epsilon);

         if(gNormalApprox) {
            epsilon /= 50;
         }

         //printf("i %d j %d\tw/w %lf epsilon %lf\n", 
         //i, j, thisWSquare[j] / (thisW[j] * thisW[j]), epsilon);
      } else {
         epsilon = 0;
      }

      thisWx = wx[i];

      for(j = 0 ; j < ExampleSpecGetNumAttributes(es) ; j++) {

         if(pedroBound) {
            assignmentError = max(deltaPlus[i][j], deltaMinus[i][j]) / 
                         wMinus[i];
            if(gMessageLevel > 2) {
               printf("   assignment bound cluster %d dim %d: %.2f\n", i, j,
                        assignmentError);
            }
         } else {
            assignmentError = 0;
         }

         if(ExampleSpecIsAttributeContinuous(es, j)) {
            if(!(w[i] <= 0 || (boundUpper && (w[i] - epsilon <= 0)))) {
               if(boundUpper) {
                  ExampleSetContinuousAttributeValue(newCenter, j,
                      (thisWx[j] / w[i]) + epsilon + assignmentError);
               } else { /* bound lower */
                  /* HERE remove this min at some point */
                  ExampleSetContinuousAttributeValue(newCenter, j,
                      max((thisWx[j] / w[i]) - (epsilon + assignmentError),
                           0));
               }

               if(useSampleError) {
                  /* make the assignment error only center */
                  if(boundUpper) {
                     ExampleSetContinuousAttributeValue(newCenterAssignment, j,
                         (thisWx[j] / w[i]) + assignmentError);
                  } else {
                  /* HERE remove this min at some point */
                     ExampleSetContinuousAttributeValue(newCenterAssignment, j,
                         max((thisWx[j] / w[i]) - assignmentError, 0));
                  }
               }
            } else {
               /* not enough weight so leave the center the same */
               if(gMessageLevel > 1) {
                  printf("not enough weight: %f epsilon: %f, possible div by zero c: %d a: %d!\n",
                                w[i], epsilon, i,j);
               }
               ExampleSetContinuousAttributeValue(newCenter, j,
                 ExampleGetContinuousAttributeValue(VALIndex(centers, i),j));
            }
         } else {
            /* HERE what to do about discrete attributes */
         }
      }

      VALAppend(newCenters, newCenter);
      if(useSampleError) {
         VALAppend(newAssignment, newCenterAssignment);
      }
   }
}

float AssignmentScaledDeltaMax(ExamplePtr e,
                          ExamplePtr centroid, ExamplePtr min, 
                           ExamplePtr max, float epsilon);
float AssignmentScaledDeltaMin(ExamplePtr e,
                          ExamplePtr centroid, ExamplePtr min, 
                           ExamplePtr max, float epsilon);

static void _RecordPedroBoundInfo(FILE *boundData, IterationStatsPtr is, VoidAListPtr newCentroids, ExampleSpecPtr es) {