New split code.
authorMonty <xiphmont@xiph.org>
Thu, 9 Dec 1999 10:16:13 +0000 (10:16 +0000)
committerMonty <xiphmont@xiph.org>
Thu, 9 Dec 1999 10:16:13 +0000 (10:16 +0000)
Monty

svn path=/trunk/vorbis/; revision=183

vq/vqgen.c

index 06f4a71..21c3b00 100644 (file)
    pregenerated codebook that is compiled into libvorbis.  It is an
    expensive (but good) algorithm.  Run it on big iron. */
 
+/* There are so many optimizations to explore in *both* stages that
+   considering the undertaking is almost withering.  For now, we brute
+   force it all */
+
 #include <stdlib.h>
 #include <stdio.h>
 #include <math.h>
 #include <string.h>
 #include "minit.h"
 
-/************************************************************************
+/* Codebook generation happens in two steps: 
+
+   1) Train the codebook with data collected from the encoder: We use
+   one of a few error metrics (which represent the distance between a
+   given data point and a candidate point in the training set) to
+   divide the training set up into cells representing roughly equal
+   probability of occurring. 
+
+   2) Generate the codebook and auxiliary data from the trained data set
+*/
+
+/* Codebook training ****************************************************
+ *
  * The basic idea here is that a VQ codebook is like an m-dimensional
  * foam with n bubbles.  The bubbles compete for space/volume and are
  * 'pressurized' [biased] according to some metric.  The basic alg
@@ -88,6 +104,8 @@ double _dist_sq(vqgen *v,double *a, double *b){
   }
   return acc;
 }
+
+/* A metric for LSP codes */
                             /* candidate,actual */
 double _dist_and_pos(vqgen *v,double *b, double *a){
   int i;
@@ -163,7 +181,6 @@ double vqgen_iterate(vqgen *v){
   double *new=malloc(sizeof(double)*v->entries*v->elements);
   long   *nearcount=malloc(v->entries*sizeof(long));
   double *nearbias=malloc(v->entries*desired*sizeof(double));
-  long biascount;
 
 #ifdef NOISY
   char buff[80];
@@ -178,7 +195,7 @@ double vqgen_iterate(vqgen *v){
   bias=fopen(buff,"w");
 #endif
 
-  fprintf(stderr,"Pass #%ld... ",v->it);
+  fprintf(stderr,"Pass #%d... ",v->it);
 
   if(v->entries<2){
     fprintf(stderr,"generation requires at least two entries\n");
@@ -190,8 +207,9 @@ double vqgen_iterate(vqgen *v){
   memset(nearcount,0,sizeof(long)*v->entries);
   memset(v->assigned,0,sizeof(long)*v->entries);
   for(i=0;i<v->points;i++){
-    double firstmetric=v->metric_func(v,_now(v,0),_point(v,i))+v->bias[0];
-    double secondmetric=v->metric_func(v,_now(v,1),_point(v,i))+v->bias[1];
+    double *ppt=_point(v,i);
+    double firstmetric=v->metric_func(v,_now(v,0),ppt)+v->bias[0];
+    double secondmetric=v->metric_func(v,_now(v,1),ppt)+v->bias[1];
     long   firstentry=0;
     long   secondentry=1;
     if(firstmetric>secondmetric){
@@ -304,228 +322,260 @@ double vqgen_iterate(vqgen *v){
   return(asserror);
 }
 
-/* the additional fields are the hyperplanes we've already used to
-   split the set.  The additional hyperplanes are necessary to bound
-   later splits.  One per depth */
-
-int lp_split(vqgen *v,long *index,long points,
-             double *additional_a, double *additional_b, long depth){
-  static int frameno=0;
-  long A[points];
-  long B[points];
-  long ap=0;
-  long bp=0;
-  long cp=0;
+/* Building a codebook from trained set **********************************
+
+   The codebook in raw form is technically finished once it's trained.
+   However, we want to finalize the representative codebook values for
+   each entry and generate auxiliary information to optimize encoding.
+   We generate the auxiliary coding tree using collected data,
+   probably the same data as in the original training */
+
+/* At each recursion, the data set is split in half.  Cells with data
+   points on side A go into set A, same with set B.  The sets may
+   overlap.  If the cell overlaps the deviding line only very slightly
+   (provided parameter), we may choose to ignore the overlap in order
+   to pare the tree down */
+
+double *sortvals;
+int els;
+int iascsort(const void *a,const void *b){
+  double av=sortvals[*((long *)a) * els];
+  double bv=sortvals[*((long *)b) * els];
+  if(av<bv)return(-1);
+  return(1);
+}
+
+/* goes through the split, but just counts it and returns a metric*/
+void lp_count(vqgen *v,long *entryindex,long entries, 
+               long *pointindex,long points,
+               long *entryA,long *entryB,
+               double *n, double c, double slack,
+               long *entriesA,long *entriesB,long *entriesC){
   long i,j,k;
-  long leaves=0;
+  long A=0,B=0,C=0;
 
-  double best=0.;
-  long besti=-1;
-  long bestj=-1;
+  memset(entryA,0,sizeof(long)*entries);
+  memset(entryB,0,sizeof(long)*entries);
 
-  if(depth>7){
-    printf("leaf: points %ld, depth %ld\n",points,depth);
-    return(1);
+  for(i=0;i<points;i++){
+    double *ppt=_point(v,pointindex[i]);
+    long   firstentry=0;
+    double firstmetric=_dist_sq(v,_now(v,entryindex[0]),ppt);
+    double position=-c;
+    
+    for(j=1;j<entries;j++){
+      double thismetric=_dist_sq(v,_now(v,entryindex[j]),ppt);
+      if(thismetric<firstmetric){
+       firstmetric=thismetric;
+       firstentry=j;
+      }
+    }
+
+    /* count point split */
+    for(k=0;k<v->elements;k++)
+      position+=ppt[k]*n[k];
+    if(position>0.){
+      entryA[firstentry]++;
+    }else{
+      entryB[firstentry]++;
+    }
+  }
+
+  /* look to see if entries are in the slack zone */
+  /* The entry splitting isn't total, so that storage has to be
+     allocated for recursion.  Reuse the entryA/entryB vectors */
+  for(j=0;j<entries;j++){
+    long total=entryA[j]+entryB[j];
+    if((double)entryA[j]/total<slack){
+      entryA[j]=0;
+    }else if((double)entryB[j]/total<slack){
+      entryB[j]=0;
+    }
+    if(entryA[j] && entryB[j])C++;
+    if(entryA[j])entryA[A++]=entryindex[j];
+    if(entryB[j])entryB[B++]=entryindex[j];
   }
+  *entriesA=A;
+  *entriesB=B;
+  *entriesC=C;
+}
+
+void pq_in_out(vqgen *v,double *n,double *c,double *p,double *q){
+  int k;
+  *c=0.;
+  for(k=0;k<v->elements;k++){
+    double center=(p[k]+q[k])/2.;
+    n[k]=(center-q[k])*2.;
+    *c+=center*n[k];
+  }
+}
 
-  if(points==1){
-    printf("leaf: point %ld, depth %ld\n",index[0],depth);
+void pq_center_out(vqgen *v,double *n,double *c,double *center,double *q){
+  int k;
+  *c=0.;
+  for(k=0;k<v->elements;k++){
+    n[k]=(center[k]-q[k])*2.;
+    *c+=center[k]*n[k];
+  }
+}
+
+int lp_split(vqgen *v,long *entryindex,long entries, 
+            long *pointindex,long points,long depth,
+            double slack){
+
+  /* The encoder, regardless of book, will be using a straight
+     euclidian distance-to-point metric to determine closest point.
+     Thus we split the cells using the same (we've already trained the
+     codebook set spacing and distribution using special metrics and
+     even a midpoint division won't disturb the basic properties) */
+
+  long ret;
+  double *p;
+  double *q;
+  double *n;
+  double c;
+  long *entryA=calloc(entries,sizeof(long));
+  long *entryB=calloc(entries,sizeof(long));
+  long entriesA=0;
+  long entriesB=0;
+  long entriesC=0;
+  long pointsA=0;
+  long i,j,k;
+
+  p=alloca(sizeof(double)*v->elements);
+  q=alloca(sizeof(double)*v->elements);
+  n=alloca(sizeof(double)*v->elements);
+  memset(p,0,sizeof(double)*v->elements);
+
+  /* depth limit */
+  if(depth>20){
+    printf("leaf: entries %ld, depth %ld\n",entries,depth);
+    return(entries);
+  }
+
+  /* nothing to do */
+  if(entries==1){
+    printf("leaf: entry %ld, depth %ld\n",entryindex[0],depth);
     return(1);
   }
 
-  if(points==2){
-    /* The result must be an even split */
-    B[bp++]=index[0];
-    A[ap++]=index[1];
+  /* The result must be an even split */
+  if(entries==2){
     printf("even split: depth %ld\n",depth);
     return(2);
-  }else{
-    /* We need to find the best split */
-    for(i=0;i<points-1;i++){
-      for(j=i+1;j<points;j++){
-       if(_dist_sq(v,_now(v,index[i]),_now(v,index[j]))>best){
+  }
+
+  /* We need to find the dividing hyperplane. find the median of each
+     axis as the centerpoint and the normal facing farthest point */
+
+  /* more than one way to do this part.  For small sets, we can brute
+     force it. */
+
+  if(entries<64){
+    /* try every pair possibility */
+    double best=0;
+    long   besti=0;
+    long   bestj=0;
+    double this;
+    for(i=0;i<entries-1;i++){
+      for(j=i+1;j<entries;j++){
+       pq_in_out(v,n,&c,_now(v,entryindex[i]),_now(v,entryindex[j]));
+       lp_count(v,entryindex,entries, 
+                pointindex,points,
+                entryA,entryB,
+                n, c, slack,
+                &entriesA,&entriesB,&entriesC);
+       this=(entriesA-entriesC)*(entriesB-entriesC);
+
+       if(this>best){
+         best=this;
          besti=i;
          bestj=j;
-         best=_dist_sq(v,_now(v,index[i]),_now(v,index[j]));
        }
       }
     }
-
-    /* We have our endpoints. initial divvy */
-    {
-      double *cA=malloc(v->elements*sizeof(double));
-      double *cB=malloc(v->elements*sizeof(double));
-      double **a=malloc((points+depth)*sizeof(double *));
-      double *b=malloc((points+depth)*sizeof(double));
-      double *aa=malloc((points+depth)*v->elements*sizeof(double));
-      double dA=0.,dB=0.;
-      minit l;
-      
-      minit_init(&l,v->elements,points+depth-1,0);
-      
-      /* divisor hyperplane */
-      for(i=0;i<v->elements;i++){
-       double m=(_now(v,index[besti])[i]+_now(v,index[bestj])[i])/2;
-       cA[i]=_now(v,index[besti])[i]-_now(v,index[bestj])[i];
-       cB[i]=-cA[i];
-       dA+=cA[i]*m;
-       dB-=cA[i]*m;
+    pq_in_out(v,n,&c,_now(v,entryindex[besti]),_now(v,entryindex[bestj]));
+  }else{
+    double best=0.;
+    long   bestj=0;
+
+    /* try COG/normal and furthest pairs */
+    /* medianpoint */
+    for(k=0;k<v->elements;k++){
+      /* just sort the index array */
+      sortvals=v->pointlist+k;
+      els=v->elements;
+      qsort(pointindex,points,sizeof(long),iascsort);
+      if(points&0x1){
+       p[k]=v->pointlist[(pointindex[points/2])*v->elements+k];
+      }else{
+       p[k]=(v->pointlist[(pointindex[points/2])*v->elements+k]+
+             v->pointlist[(pointindex[points/2-1])*v->elements+k])/2.;
       }
-      
-      for(i=0;i<points+depth;i++)
-       a[i]=aa+i*v->elements;
-      
-      /* check each bubble to see if it intersects the divisor hyperplane */
-      for(j=0;j<points;j++){
-       long count=0,m=0;
-       double d1=0.,d2=0.;
-       int ret;
-       
-       if(j==besti){
-         B[bp++]=index[j];
-       }else if (j==bestj){
-         A[ap++]=index[j];
-       }else{
-        
-         for(i=0;i<points;i++){
-           if(i!=j){
-             b[count]=0.;
-             for(k=0;k<v->elements;k++){
-               double m=(_now(v,index[i])[k]+_now(v,index[j])[k])/2;
-               a[count][k]=_now(v,index[i])[k]-_now(v,index[j])[k];
-               b[count]+=a[count][k]*m;
-             }
-             count++;
-           }
-         }
-         
-         /* additional bounding hyperplanes from previous splits */
-         for(i=0;i<depth;i++){
-           b[count]=0.;
-           for(k=0;k<v->elements;k++)
-             a[count][k]=additional_a[m++];
-           b[count++]=additional_b[i];
-         }
-         
-         /* on what side of the dividing hyperplane is the current test
-            point? */
-         for(i=0;i<v->elements;i++)
-           d1+=_now(v,index[j])[i]*cA[i];
-         
-         
-         if(d1<dA){
-           fprintf(stderr,"+");
-           ret=minit_solve(&l,a,b,cA,1e-6,NULL,NULL,&d2);
-         }else{
-           fprintf(stderr,"-");
-           ret=minit_solve(&l,a,b,cB,1e-6,NULL,NULL,&d2);
-           d2= -d2;
-         }
-         
-         switch(ret){
-         case 0:
-           /* bounded solution */
-           if(d1<dA){
-             A[ap++]=index[j];
-             if(d2>dA){cp++;B[bp++]=index[j];}
-           }else{
-             B[bp++]=index[j];
-             if(d2<dA){cp++;A[ap++]=index[j];}
-           }
-           break;
-         default:
-           /* unbounded solution or no solution; we're in both sets */
-           B[bp++]=index[j];
-           A[ap++]=index[j];
-           cp++;
-           break;
-         }
-       }
+    }
+    
+    /* try every normal, but just for distance */
+    for(j=0;j<entries;j++){
+      double *ppj=_now(v,entryindex[j]);
+      double this=_dist_sq(v,p,ppj);
+      if(this>best){
+       best=this;
+       bestj=j;
       }
+    }
 
-      /*{
-       char buffer[80];
-       FILE *o;
-       int i;
-
-       sprintf(buffer,"set%d.m",frameno);
-       o=fopen(buffer,"w");
-       for(i=0;i<points;i++)
-         fprintf(o,"%g %g\n\n",_now(v,index[i])[0],_now(v,index[i])[1]);
-       fclose(o);
-
-       sprintf(buffer,"setA%d.m",frameno);
-       o=fopen(buffer,"w");
-       for(i=0;i<ap;i++)
-         fprintf(o,"%g %g\n\n",_now(v,A[i])[0],_now(v,A[i])[1]);
-       fclose(o);
-
-       sprintf(buffer,"setB%d.m",frameno);
-       o=fopen(buffer,"w");
-       for(i=0;i<bp;i++)
-         fprintf(o,"%g %g\n\n",_now(v,B[i])[0],_now(v,B[i])[1]);
-       fclose(o);
-
-       sprintf(buffer,"div%d.m",frameno);
-       o=fopen(buffer,"w");
-       fprintf(o,"%g %g\n%g %g\n\n",
-               _now(v,index[besti])[0],
-               (dA-cA[0]*_now(v,index[besti])[0])/cA[1],
-               _now(v,index[bestj])[0],
-               (dA-cA[0]*_now(v,index[bestj])[0])/cA[1]);
-       fclose(o);
-
-       sprintf(buffer,"bound%d.m",frameno);
-       o=fopen(buffer,"w");
-       for(i=0;i<depth;i++)
-         fprintf(o,"%g %g\n%g %g\n\n",
-                 _now(v,index[besti])[0],
-                 (additional_b[i]-
-                  additional_a[i*2]*_now(v,index[besti])[0])/
-                 additional_a[i*2+1],
-                 _now(v,index[bestj])[0],
-                 (additional_b[i]-
-                  additional_a[i*2]*_now(v,index[bestj])[0])/
-                 additional_a[i*2+1]);
-       fclose(o);
-       frameno++;
-       }*/
-
-
-      minit_clear(&l);
-      free(a);
-      free(b);
-      free(aa);
-
-      fprintf(stderr,"split: total=%ld depth=%ld set A=%ld:%ld:%ld=B\n",
-             points,depth,ap-cp,cp,bp-cp);
-
-      /* handle adding this divisor hyperplane to the additional
-         constraints.  The 'inside' side is different for A and B */
-
-      {
-       /* marginally wasteful */
-       double *newadd_a=alloca(sizeof(double)*v->elements*(depth+1));
-       double *newadd_b=alloca(sizeof(double)*(depth+1));
-       if(additional_a){
-         memcpy(newadd_a,additional_a,sizeof(double)*v->elements*depth);
-         memcpy(newadd_b,additional_b,sizeof(double)*depth);
-       }
+    pq_center_out(v,n,&c,p,_point(v,pointindex[bestj]));
 
-       memcpy(newadd_a+v->elements*depth,cA,sizeof(double)*v->elements);
-       newadd_b[depth]=dA;     
-       leaves=lp_split(v,A,ap,newadd_a,newadd_b,depth+1);
 
-       memcpy(newadd_a+v->elements*depth,cB,sizeof(double)*v->elements);
-       newadd_b[depth]=dB;     
-       leaves+=lp_split(v,B,bp,newadd_a,newadd_b,depth+1);
-      
+  }
+
+  /* find cells enclosing points */
+  /* count A/B points */
+
+  lp_count(v,entryindex,entries, 
+          pointindex,points,
+          entryA,entryB,
+          n, c, slack,
+          &entriesA,&entriesB,&entriesC);
+
+  /* the point index is split evenly, so we do an Order n
+     rearrangement into A first/B last and just pass it on */
+  {
+    long Aptr=0;
+    long Bptr=points-1;
+    while(Aptr<Bptr){
+      while(Aptr<Bptr){
+       double position=-c;
+       for(k=0;k<v->elements;k++)
+         position+=_point(v,pointindex[Aptr])[k]*n[k];
+       if(position<0.)break; /* not in A */
+       Aptr++;
+      }
+      while(Aptr<Bptr){
+       double position=-c;
+       for(k=0;k<v->elements;k++)
+         position+=_point(v,pointindex[Bptr])[k]*n[k];
+       if(position>0.)break; /* not in B */
+       Bptr--;
+      }
+      if(Aptr<Bptr){
+       long temp=pointindex[Aptr];
+       pointindex[Aptr]=pointindex[Bptr];
+       pointindex[Bptr]=temp;
       }
-      free(cA);
-      free(cB);    
-      return(leaves);
+      pointsA=Aptr;
+      Aptr++;
+      Bptr--;
     }
   }
+
+  fprintf(stderr,"split: total=%ld depth=%ld set A=%ld:%ld:%ld=B\n",
+         entries,depth,entriesA-entriesC,entriesC,entriesB-entriesC);
+
+  ret=lp_split(v,entryA,entriesA,pointindex,pointsA,depth+1,slack); 
+  ret+=lp_split(v,entryB,entriesB,pointindex+pointsA,points-pointsA,
+               depth+1,slack); 
+  return(ret);
 }
 
 void vqgen_book(vqgen *v){
@@ -536,289 +586,267 @@ void vqgen_book(vqgen *v){
 }
 
 static double testset24[48]={
-1.047758,1.245406,
-1.007972,1.150078,
-1.064390,1.146266,
-0.761842,0.826055,
-0.862073,1.243707,
-0.746351,1.055368,
-0.956844,1.048223,
-0.877782,1.111871,
-0.964564,1.309219,
-0.920062,1.168658,
-0.787654,0.895163,
-0.974003,1.200115,
-0.788360,1.142671,
-0.978414,1.122249,
-0.869938,0.954436,
-1.131220,1.348747,
-0.934630,1.108564,
-0.896666,1.045772,
-0.707360,0.920954,
-0.788914,0.970626,
-0.828682,1.043673,
-1.042016,1.190406,
-1.031643,1.068194,
-1.120824,1.220326,
+
 };
 
 static double testset256[1024]={
-1.047758,1.245406,1.007972,1.150078,
-1.064390,1.146266,0.761842,0.826055,
-0.862073,1.243707,0.746351,1.055368,
-0.956844,1.048223,0.877782,1.111871,
-0.964564,1.309219,0.920062,1.168658,
-0.787654,0.895163,0.974003,1.200115,
-0.788360,1.142671,0.978414,1.122249,
-0.869938,0.954436,1.131220,1.348747,
-0.934630,1.108564,0.896666,1.045772,
-0.707360,0.920954,0.788914,0.970626,
-0.828682,1.043673,1.042016,1.190406,
-1.031643,1.068194,1.120824,1.220326,
-0.747084,0.888176,1.051535,1.221335,
-0.980283,1.087708,1.311154,1.445530,
-0.927548,1.107237,1.296902,1.480610,
-0.804415,0.921948,1.029952,1.231078,
-0.923413,1.064919,1.143109,1.242755,
-0.826403,1.010624,1.353139,1.455870,
-0.837930,1.073739,1.251770,1.405370,
-0.704437,0.890366,1.016676,1.136689,
-0.840924,1.151858,1.369346,1.558572,
-0.835229,1.095713,1.165585,1.311155,
-0.763157,1.000718,1.168780,1.345318,
-0.912852,1.040279,1.260220,1.468606,
-0.680839,0.920854,1.056199,1.214924,
-0.993572,1.124795,1.220880,1.416346,
-0.766878,0.994942,1.170035,1.432164,
-1.006875,1.069566,1.283594,1.531742,
-1.029033,1.188357,1.330626,1.519446,
-0.915581,1.121514,1.193260,1.431981,
-0.950939,1.117674,1.280092,1.487413,
-0.783569,0.863166,1.198348,1.309199,
-0.845209,0.949087,1.133299,1.253263,
-0.836299,1.050413,1.263712,1.482634,
-1.049458,1.194466,1.332964,1.485894,
-0.921516,1.086016,1.244810,1.349821,
-0.888983,1.111494,1.313032,1.351122,
-1.015785,1.142849,1.292746,1.463028,
-0.860002,1.249262,1.354167,1.534761,
-0.942136,1.053546,1.120840,1.373279,
-0.846661,1.059871,1.294814,1.504904,
-0.796592,1.061316,1.349680,1.494821,
-0.880504,1.033570,1.314379,1.484966,
-0.770438,1.039993,1.238670,1.412317,
-1.036503,1.068649,1.202522,1.487553,
-0.942615,1.056630,1.344355,1.562161,
-0.963759,1.096427,1.324085,1.497266,
-0.899956,1.206083,1.349655,1.572202,
-0.977023,1.132722,1.338437,1.488086,
-1.096643,1.242143,1.357079,1.513244,
-0.795193,1.011915,1.312761,1.419554,
-1.021869,1.075269,1.265230,1.475986,
-1.019546,1.186088,1.347798,1.567607,
-1.051191,1.172904,1.381952,1.511843,
-0.924332,1.192953,1.385378,1.420727,
-0.869539,0.962051,1.400803,1.539414,
-0.949875,1.053939,1.305695,1.422764,
-0.962223,1.085902,1.247702,1.462327,
-0.991777,1.182062,1.287421,1.508453,
-0.918274,1.107269,1.325467,1.454377,
-0.993589,1.184665,1.352294,1.475413,
-0.804960,1.011490,1.259031,1.462335,
-0.845949,1.173858,1.389681,1.441148,
-0.858168,1.019930,1.232463,1.440218,
-0.884411,1.206226,1.318968,1.457990,
-0.893469,1.074271,1.274410,1.433475,
-1.039654,1.061160,1.122426,1.189133,
-1.062243,1.177357,1.344765,1.542071,
-0.900964,1.116803,1.380397,1.419108,
-0.969535,1.216066,1.383582,1.435924,
-1.017913,1.330986,1.403452,1.450709,
-0.771534,0.973324,1.171604,1.236873,
-0.955683,1.029609,1.309284,1.498449,
-0.892165,1.192347,1.322444,1.356025,
-0.890209,1.023145,1.326649,1.431669,
-1.038788,1.188891,1.330844,1.436180,
-1.050355,1.177363,1.283447,1.479040,
-0.825367,0.960892,1.179047,1.381923,
-0.978931,1.122163,1.304581,1.532161,
-0.970750,1.149650,1.390225,1.547422,
-0.805416,0.945116,1.118815,1.343870,
-1.004106,1.185996,1.304086,1.464108,
-0.771591,1.055865,1.150908,1.228683,
-0.950175,1.104237,1.203437,1.387264,
-0.923504,1.122814,1.173989,1.318650,
-0.960635,1.042342,1.196995,1.395394,
-1.035111,1.064804,1.134387,1.341424,
-0.815182,1.083612,1.159518,1.255777,
-0.784079,0.830594,1.056005,1.289902,
-0.737774,0.994134,1.115470,1.264879,
-0.820162,1.105657,1.182947,1.423629,
-1.067552,1.142869,1.191626,1.387014,
-0.752416,0.918026,1.124637,1.414500,
-0.767963,0.840069,0.997625,1.325653,
-0.787595,0.865440,1.090071,1.227348,
-0.798877,1.105239,1.331379,1.440643,
-0.829079,1.133632,1.280774,1.470690,
-1.009195,1.132959,1.284044,1.500589,
-0.698569,0.824611,1.236682,1.462088,
-0.817460,0.985767,1.081910,1.257751,
-0.784033,0.882552,1.149678,1.326920,
-1.039403,1.085310,1.211033,1.558635,
-0.966795,1.168857,1.309960,1.497788,
-0.906244,1.027050,1.198846,1.323671,
-0.776495,1.185285,1.309167,1.380683,
-0.799628,0.927976,1.048238,1.299045,
-0.779808,0.881990,1.167898,1.220314,
-0.770446,0.918281,1.049189,1.179393,
-0.768416,1.037558,1.165760,1.302606,
-0.743121,0.814671,0.990501,1.224236,
-1.037050,1.068240,1.159690,1.426280,
-0.978810,1.214329,1.253336,1.324395,
-0.984003,1.121443,1.376382,1.510519,
-1.146510,1.229726,1.417616,1.781032,
-0.897163,1.147910,1.221186,1.371815,
-1.068525,1.211553,1.343551,1.506743,
-0.762313,1.091082,1.253251,1.472381,
-0.960562,1.041965,1.247053,1.399214,
-0.864482,1.123473,1.163412,1.238620,
-0.963484,1.132803,1.164992,1.250389,
-1.009456,1.139510,1.251339,1.449078,
-0.851837,1.113642,1.170290,1.362806,
-0.857073,0.962039,1.127381,1.471682,
-1.047754,1.213381,1.388899,1.492383,
-0.938921,1.267308,1.337076,1.478427,
-0.790388,0.912816,1.159450,1.273259,
-0.832690,0.997776,1.156639,1.302621,
-0.783009,0.975374,1.080630,1.311112,
-0.819784,1.145093,1.326949,1.525480,
-0.806394,1.089564,1.329564,1.550362,
-0.958608,1.036364,1.379118,1.443043,
-0.753680,0.941781,1.147749,1.297211,
-0.883974,1.091394,1.315093,1.480307,
-1.239591,1.417601,1.495843,1.641941,
-0.881068,0.973780,1.278918,1.429384,
-0.869052,0.977661,1.280744,1.551295,
-1.022685,1.052986,1.105046,1.168670,
-0.981698,1.131448,1.197781,1.467704,
-0.945034,1.163410,1.250872,1.428793,
-1.092055,1.139380,1.187113,1.264586,
-0.788897,0.953764,1.232844,1.506461,
-0.885206,0.978419,1.209467,1.387569,
-0.762835,0.964279,1.064844,1.243153,
-0.974906,1.134763,1.283544,1.460386,
-1.081114,1.169260,1.290987,1.535452,
-0.880531,1.075660,1.378980,1.522978,
-1.039431,1.083792,1.350481,1.588522,
-0.922449,1.060622,1.152309,1.467930,
-0.918935,1.073732,1.283243,1.479567,
-0.992722,1.145398,1.239500,1.509004,
-0.903405,1.243062,1.421386,1.572148,
-0.912531,1.100239,1.310219,1.521424,
-0.875168,0.912787,1.123618,1.363635,
-0.732241,1.108317,1.323119,1.515675,
-0.951518,1.141874,1.341623,1.533409,
-0.992099,1.215801,1.477413,1.734691,
-1.005267,1.122655,1.356975,1.436445,
-0.740659,0.807749,1.167728,1.380618,
-1.078727,1.214706,1.328173,1.556699,
-1.026027,1.168906,1.313249,1.486078,
-0.914877,1.147150,1.389489,1.523984,
-1.062339,1.120684,1.160024,1.234794,
-1.129141,1.197655,1.374217,1.547755,
-1.070011,1.255271,1.360225,1.467869,
-0.779980,0.871696,1.098031,1.284490,
-1.062409,1.177542,1.314581,1.696662,
-0.702935,1.229873,1.370813,1.479500,
-1.029357,1.225167,1.341607,1.478163,
-1.025666,1.141749,1.185959,1.332892,
-0.799462,0.951470,1.214070,1.305787,
-0.740521,0.805457,1.107504,1.317258,
-0.784194,0.838683,1.055934,1.242692,
-0.839416,1.048060,1.391801,1.623786,
-0.692627,0.907677,1.060843,1.341002,
-0.823625,1.244497,1.396901,1.586243,
-0.942859,1.232978,1.348170,1.536735,
-0.894882,1.131376,1.292892,1.462724,
-0.776974,0.904449,1.325557,1.451968,
-1.066188,1.218328,1.376282,1.545381,
-0.990053,1.124279,1.340534,1.559666,
-0.776321,0.935169,1.191081,1.372326,
-0.949935,1.115929,1.397704,1.442549,
-0.936357,1.126885,1.356247,1.579931,
-1.045433,1.274605,1.366947,1.590215,
-1.063890,1.138062,1.220645,1.460005,
-0.751448,0.811338,1.078027,1.403146,
-0.935678,1.102858,1.356557,1.515948,
-1.026417,1.143843,1.299309,1.413976,
-0.821475,1.000237,1.105073,1.379882,
-0.960249,1.031602,1.250804,1.462620,
-0.661745,1.106452,1.291188,1.439529,
-0.691661,0.947241,1.261183,1.457391,
-0.839024,0.914750,1.040695,1.375853,
-1.029401,1.065349,1.121370,1.270670,
-0.888316,1.003349,1.103749,1.341290,
-0.766328,0.879083,1.217779,1.419868,
-0.811672,1.049673,1.186460,1.375742,
-0.969585,1.170126,1.259338,1.595070,
-1.016617,1.145706,1.335413,1.503064,
-0.980227,1.295316,1.417964,1.478684,
-0.885701,1.248099,1.416821,1.465338,
-0.953583,1.266596,1.325829,1.372230,
-1.038619,1.225860,1.351664,1.527048,
-1.104724,1.215839,1.250491,1.335237,
-1.124449,1.269485,1.420756,1.677260,
-0.881337,1.352259,1.433218,1.492756,
-1.023097,1.059205,1.110249,1.212976,
-1.135632,1.282243,1.415540,1.566039,
-1.063524,1.252342,1.399082,1.518433,
-0.790658,0.856337,1.154909,1.274963,
-1.119059,1.382625,1.423651,1.492741,
-1.030526,1.296610,1.330511,1.396550,
-0.947952,1.065235,1.225274,1.554455,
-0.960669,1.282313,1.370362,1.572736,
-0.996042,1.208193,1.385036,1.534877,
-1.003206,1.377432,1.431110,1.497189,
-1.088166,1.227944,1.508129,1.740407,
-0.996566,1.162407,1.347665,1.524235,
-0.944606,1.287026,1.400822,1.437156,
-1.066144,1.238314,1.454451,1.616016,
-1.121065,1.200875,1.316542,1.459583,
-1.158944,1.271448,1.356823,1.450510,
-0.811670,1.278011,1.361550,1.440647,
-0.875620,1.103051,1.230854,1.483429,
-0.959882,1.180685,1.381224,1.572807,
-1.049222,1.186513,1.387883,1.567423,
-1.049920,1.293154,1.507194,1.678371,
-0.872138,1.193727,1.455817,1.665837,
-0.738018,0.946583,1.335543,1.552061,
-1.072856,1.151838,1.321877,1.486028,
-1.026909,1.153575,1.306700,1.550408,
-0.940629,1.121943,1.303228,1.409285,
-1.023544,1.269480,1.316695,1.508732,
-0.924872,1.119626,1.479602,1.698360,
-0.793358,1.160361,1.381276,1.636124,
-0.892504,1.248046,1.320256,1.367114,
-1.165222,1.230265,1.507234,1.820748,
-0.711084,0.876165,1.145850,1.493303,
-0.776482,0.869670,1.124462,1.370831,
-0.785085,0.912534,1.099328,1.281523,
-0.947788,1.168239,1.407409,1.473592,
-0.913992,1.316132,1.619366,1.865030,
-1.024592,1.151907,1.383557,1.545112,
-0.987748,1.168194,1.415323,1.639563,
-0.725008,1.037088,1.306659,1.474607,
-0.965243,1.263523,1.454521,1.638929,
-1.035169,1.219404,1.427315,1.578414,
-0.787215,0.960639,1.245517,1.423549,
-0.902895,1.153039,1.412772,1.613360,
-0.983205,1.338886,1.493960,1.619623,
-0.822576,1.010524,1.232397,1.359007,
-0.773450,1.090005,1.170248,1.369534,
-1.117781,1.354485,1.662228,1.829565,
-1.104316,1.324945,1.525085,1.741040,
-1.135275,1.371208,1.537082,1.673949,
-0.889899,1.195451,1.242648,1.291881,
+0.334427,0.567149,0.749324,0.838460,
+0.212558,0.435792,0.661945,0.781195,
+0.316791,0.551190,0.700236,0.813437,
+0.240661,0.398754,0.570649,0.666909,
+0.233216,0.549751,0.715668,0.844393,
+0.275589,0.460791,0.687872,0.786855,
+0.145195,0.478688,0.678451,0.788536,
+0.332491,0.577111,0.770913,0.875295,
+0.162683,0.250422,0.399770,0.688120,
+0.301578,0.424074,0.595725,0.829533,
+0.118704,0.276233,0.492329,0.605532,
+0.240679,0.468862,0.670704,0.809156,
+0.109799,0.227123,0.533265,0.686982,
+0.119657,0.423823,0.597676,0.710405,
+0.207444,0.437121,0.671331,0.809054,
+0.222787,0.324523,0.466905,0.652758,
+0.117279,0.268017,0.410036,0.593785,
+0.316531,0.493867,0.642588,0.815791,
+0.218069,0.315630,0.554662,0.749348,
+0.152218,0.412864,0.664135,0.814275,
+0.144242,0.221503,0.457281,0.631403,
+0.290454,0.458806,0.603994,0.751802,
+0.236719,0.341259,0.586677,0.707653,
+0.257617,0.406941,0.602074,0.696609,
+0.141833,0.254061,0.523788,0.701270,
+0.136685,0.333920,0.515585,0.642137,
+0.191560,0.485163,0.665121,0.755073,
+0.086062,0.375207,0.623168,0.763659,
+0.184017,0.404999,0.609708,0.746118,
+0.208200,0.414804,0.618734,0.732746,
+0.184643,0.390701,0.666137,0.791445,
+0.216594,0.519656,0.713291,0.810097,
+0.118095,0.361088,0.577184,0.667553,
+0.243858,0.394640,0.550682,0.774866,
+0.266410,0.430656,0.629812,0.729035,
+0.140187,0.323258,0.534031,0.762617,
+0.097902,0.230404,0.382920,0.707036,
+0.197244,0.424752,0.578097,0.741857,
+0.216830,0.460883,0.622823,0.747350,
+0.096586,0.252731,0.471721,0.766199,
+0.385920,0.650477,0.846797,0.972222,
+0.249895,0.442701,0.706660,0.824760,
+0.207781,0.355197,0.538083,0.686972,
+0.166064,0.261252,0.493243,0.642141,
+0.146441,0.385197,0.595104,0.739789,
+0.207555,0.311523,0.533963,0.666395,
+0.172735,0.455702,0.624525,0.713503,
+0.132351,0.420175,0.528877,0.716175,
+0.113027,0.270164,0.462866,0.653668,
+0.306970,0.575262,0.785717,0.942889,
+0.122649,0.451363,0.652393,0.749624,
+0.337459,0.561671,0.713605,0.878091,
+0.154890,0.324078,0.469466,0.677888,
+0.242478,0.420602,0.651274,0.762209,
+0.162266,0.401114,0.639885,0.769255,
+0.275493,0.500799,0.753756,0.884563,
+0.219840,0.410217,0.634590,0.784979,
+0.197822,0.418074,0.514919,0.784821,
+0.176855,0.327821,0.624861,0.768299,
+0.187133,0.373303,0.614062,0.764326,
+0.148866,0.339151,0.571398,0.693017,
+0.126516,0.364169,0.643363,0.784478,
+0.176790,0.354609,0.553963,0.728466,
+0.169151,0.294951,0.559738,0.718315,
+0.189233,0.369137,0.615112,0.736337,
+0.211442,0.459206,0.587994,0.715496,
+0.203544,0.332755,0.496152,0.707495,
+0.264950,0.377661,0.601633,0.738910,
+0.144156,0.332387,0.533986,0.858065,
+0.236705,0.493254,0.646422,0.780372,
+0.247066,0.480770,0.725958,0.847099,
+0.121273,0.415837,0.575449,0.785059,
+0.109050,0.326623,0.585313,0.725416,
+0.079853,0.281243,0.579608,0.721251,
+0.164131,0.452931,0.608512,0.772505,
+0.183958,0.382443,0.595812,0.708979,
+0.232813,0.460979,0.612592,0.785098,
+0.372632,0.625661,0.794886,0.915078,
+0.197199,0.422499,0.630676,0.764479,
+0.112259,0.228495,0.423574,0.529402,
+0.313775,0.495449,0.684949,0.906438,
+0.121923,0.272895,0.590304,0.762871,
+0.114215,0.298370,0.523344,0.690963,
+0.195814,0.446435,0.567172,0.888808,
+0.292057,0.491911,0.707296,0.836253,
+0.198244,0.410661,0.555814,0.699832,
+0.152291,0.332174,0.609291,0.732500,
+0.051002,0.303083,0.609333,0.769208,
+0.086086,0.233029,0.475899,0.590504,
+0.096449,0.366927,0.502638,0.761407,
+0.256632,0.423960,0.650693,0.805566,
+0.235739,0.542293,0.776812,0.939851,
+0.210050,0.310752,0.439640,0.760752,
+0.168588,0.352547,0.620034,0.811813,
+0.217251,0.438798,0.704758,0.854804,
+0.112327,0.304063,0.583490,0.813809,
+0.183413,0.263207,0.503397,0.753610,
+0.241265,0.367568,0.648774,0.795899,
+0.262796,0.558878,0.784660,0.908486,
+0.221425,0.479136,0.683794,0.830981,
+0.325039,0.573201,0.837740,0.994924,
+0.263655,0.636322,0.878818,1.052555,
+0.228884,0.510746,0.742781,0.903443,
+0.169799,0.354550,0.667261,0.880659,
+0.429209,0.689908,0.942083,1.075538,
+0.178843,0.393247,0.698620,0.846058,
+0.410735,0.684575,0.900212,1.026610,
+0.212082,0.407634,0.677794,0.824167,
+0.183638,0.456984,0.701680,0.828091,
+0.198668,0.290621,0.586525,0.796566,
+0.212584,0.394786,0.624992,0.756258,
+0.218545,0.327989,0.671635,0.851695,
+0.118076,0.487280,0.758482,0.920421,
+0.184608,0.426118,0.642211,0.797123,
+0.248725,0.534657,0.737666,0.869414,
+0.075312,0.218113,0.574791,0.767607,
+0.184571,0.487955,0.742725,0.865892,
+0.148818,0.284078,0.481156,0.816008,
+0.205061,0.336518,0.506794,0.610053,
+0.054274,0.238425,0.414471,0.577301,
+0.082819,0.289140,0.632090,0.789039,
+0.095066,0.249569,0.666867,0.857491,
+0.129159,0.215969,0.524985,0.808056,
+0.120949,0.473613,0.723959,0.844274,
+0.225445,0.469240,0.781258,0.937501,
+0.190737,0.487420,0.664112,0.798043,
+0.249655,0.533041,0.808386,0.984218,
+0.236544,0.473767,0.723955,0.872327,
+0.093351,0.432384,0.672316,0.803956,
+0.171295,0.491945,0.662536,0.844976,
+0.165080,0.389674,0.533508,0.664579,
+0.224994,0.356513,0.535555,0.640525,
+0.158366,0.244279,0.437682,0.579242,
+0.267249,0.423934,0.542629,0.635151,
+0.105840,0.183727,0.395015,0.546705,
+0.128701,0.214397,0.323019,0.569708,
+0.205691,0.367797,0.491331,0.579727,
+0.136376,0.295543,0.461283,0.567165,
+0.090912,0.326863,0.449785,0.549052,
+0.280052,0.463123,0.590556,0.676314,
+0.053220,0.327610,0.496921,0.611028,
+0.074906,0.374797,0.563476,0.691693,
+0.064942,0.286131,0.414479,0.492806,
+0.303210,0.504309,0.652008,0.744313,
+0.322904,0.532050,0.682097,0.777012,
+0.247437,0.397965,0.504882,0.680784,
+0.076782,0.185773,0.355825,0.478750,
+0.170565,0.428925,0.587157,0.683893,
+0.186692,0.282933,0.468043,0.583405,
+0.286847,0.489852,0.626934,0.714018,
+0.196398,0.333200,0.447551,0.551070,
+0.061757,0.242351,0.503337,0.676114,
+0.119199,0.359308,0.514024,0.595332,
+0.282369,0.428107,0.549797,0.721344,
+0.120772,0.328758,0.465806,0.797635,
+0.251676,0.385411,0.692374,0.930173,
+0.265590,0.605988,0.828534,0.964088,
+0.118910,0.347784,0.690580,0.839342,
+0.148664,0.429126,0.564376,0.640949,
+0.273239,0.395545,0.492952,0.803872,
+0.085855,0.175468,0.449037,0.649337,
+0.101341,0.351796,0.455442,0.645061,
+0.181516,0.301506,0.381633,0.628656,
+0.165340,0.385706,0.549717,0.787060,
+0.102530,0.247439,0.362454,0.570843,
+0.115626,0.291845,0.411639,0.730918,
+0.130168,0.402586,0.531264,0.842151,
+0.312057,0.497535,0.608141,0.933018,
+0.061317,0.329143,0.457434,0.717368,
+0.260419,0.364385,0.522931,0.733634,
+0.287346,0.448695,0.703337,0.864549,
+0.158190,0.329088,0.444748,0.609073,
+0.057217,0.172765,0.500069,0.763990,
+0.185582,0.346521,0.426838,0.724792,
+0.152203,0.253061,0.580445,0.835384,
+0.223681,0.376501,0.463553,0.758668,
+0.153621,0.385013,0.482490,0.723052,
+0.165862,0.300622,0.515314,0.628741,
+0.087543,0.300967,0.529657,0.641654,
+0.257693,0.367954,0.453859,0.654989,
+0.077763,0.391224,0.586210,0.846828,
+0.140653,0.296311,0.431419,0.524173,
+0.276106,0.566792,0.723587,0.928366,
+0.345768,0.530982,0.711197,1.011075,
+0.328643,0.578965,0.775015,1.036448,
+0.295489,0.451392,0.541007,0.831967,
+0.249971,0.517649,0.662078,0.872210,
+0.171802,0.382312,0.492482,0.629808,
+0.192277,0.478000,0.611564,0.846436,
+0.217769,0.335272,0.508986,0.813008,
+0.265875,0.429911,0.625892,0.958544,
+0.183329,0.299421,0.565539,0.885282,
+0.311307,0.432098,0.756371,0.937483,
+0.153908,0.264165,0.391052,0.502304,
+0.267528,0.476930,0.862271,1.099421,
+0.221527,0.513611,0.684819,0.993058,
+0.220330,0.352263,0.575808,0.823017,
+0.084526,0.372939,0.705268,0.879314,
+0.179752,0.287765,0.507176,0.683204,
+0.234005,0.381653,0.537092,0.887657,
+0.068805,0.135538,0.363455,0.732055,
+0.266267,0.440971,0.719076,1.025917,
+0.169728,0.426510,0.670186,0.982939,
+0.248443,0.463944,0.736167,0.898651,
+0.074789,0.155652,0.601987,0.864470,
+0.394203,0.665959,0.950879,1.183016,
+0.056927,0.247163,0.620166,0.886961,
+0.081110,0.260015,0.725344,0.969987,
+0.187566,0.605658,0.907862,1.113592,
+0.188174,0.417467,0.746757,0.912718,
+0.214944,0.444902,0.813762,1.006999,
+0.177305,0.322288,0.623538,0.972940,
+0.108933,0.560083,0.892489,1.105198,
+0.326010,0.486890,0.808596,1.012135,
+0.507702,0.815903,1.085892,1.267201,
+0.121027,0.368451,0.715029,1.000456,
+0.424229,0.713942,0.971720,1.123569,
+0.158148,0.411905,0.812352,1.057000,
+0.330896,0.639677,0.814721,1.083593,
+0.116143,0.314050,0.643242,0.928001,
+0.228088,0.364437,0.747642,1.045211,
+0.457076,0.783579,1.021362,1.174221,
+0.069613,0.430727,0.845847,1.061835,
+0.343281,0.553869,0.865545,1.077270,
+0.152973,0.304560,0.642854,0.830441,
+0.157602,0.568564,0.815800,1.036183,
+0.104003,0.505075,0.796932,0.999392,
+0.377220,0.622243,0.901761,1.072988,
+0.060039,0.378222,0.710910,0.939836,
+0.220206,0.539490,0.849893,1.056103,
+0.116884,0.454544,0.627779,0.924877,
+0.293087,0.633251,0.927781,1.114120,
+0.234659,0.458666,0.693916,0.925308,
+0.271802,0.385786,0.661508,0.862002,
+0.185212,0.282842,0.431237,0.663853,
+0.152441,0.280173,0.711942,0.961418,
+0.264191,0.647776,0.976337,1.171963,
+0.227340,0.567363,0.753334,1.049249,
+0.173306,0.508941,0.761828,0.970236,
+0.271942,0.496903,0.759621,0.963074,
+0.187050,0.355811,0.752185,0.952287,
+0.172003,0.363447,0.486011,0.829142,
+0.364755,0.598200,0.868222,1.025048,
+0.226504,0.432744,0.649002,0.852690,
+0.200857,0.436929,0.651778,0.905448,
+0.137860,0.385245,0.609566,0.867573,
+0.217881,0.368030,0.601419,0.954770,
+0.199237,0.398640,0.593487,0.830535,
+0.324623,0.518989,0.782275,0.926822,
+0.245185,0.419081,0.717432,0.883462,
+0.277741,0.494080,0.765001,0.915227,
+0.200418,0.463680,0.705028,0.875441,
+0.171638,0.526057,0.697163,0.915452,
+0.142351,0.416539,0.707280,0.876499,
+0.178309,0.553184,0.777786,0.903243,
+0.335361,0.540299,0.812589,0.965039,
+
 
 };
 
@@ -838,14 +866,10 @@ int main(int argc,char *argv[]){
   }
   fclose(in);
 
-  for(i=0;i<200;i++){
+  /*for(i=0;i<200;i++){
     vqgen_iterate(&v);
   }
 
-  exit(0);
-
-  memcpy(v.entrylist,testset256,sizeof(testset256));
-
   for(i=0;i<v.entries;i++){
     printf("\n");
     for(j=0;j<v.elements;j++)
@@ -853,14 +877,20 @@ int main(int argc,char *argv[]){
   }
   printf("\n");
     
+  exit(0);*/
+
+
+  memcpy(v.entrylist,testset256,sizeof(testset256));
 
-  /* try recursively splitting the space using LP 
   {
-    long index[v.entries];
-    for(i=0;i<v.entries;i++)index[i]=i;
+    long entryindex[v.entries];
+    long pointindex[v.points];
+    for(i=0;i<v.entries;i++)entryindex[i]=i;
+    for(i=0;i<v.points;i++)pointindex[i]=i;
+
     fprintf(stderr,"\n\nleaves=%d\n",
-           lp_split(&v,index,v.entries,NULL,NULL,0));
-           }*/
+           lp_split(&v,entryindex,v.entries, pointindex,v.points,0,0));
+  }
   
   return(0);
 }