#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) {