Incremental so I can continue working on a faster box ;-)
authorMonty <xiphmont@xiph.org>
Mon, 6 Dec 1999 02:55:34 +0000 (02:55 +0000)
committerMonty <xiphmont@xiph.org>
Mon, 6 Dec 1999 02:55:34 +0000 (02:55 +0000)
Monty

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

vq/vqgen.c

index 09a4646..ba745ba 100644 (file)
@@ -14,7 +14,7 @@
  function: build a VQ codebook 
  author: Monty <xiphmont@mit.edu>
  modifications by: Monty
- last modification date: Nov 18 1999
+ last modification date: Oct 04 1999
 
  ********************************************************************/
 
@@ -27,6 +27,7 @@
 #include <stdio.h>
 #include <math.h>
 #include <string.h>
+#include "minit.h"
 
 /************************************************************************
  * The basic idea here is that a VQ codebook is like an m-dimensional
  * 'pressurized' [biased] according to some metric.  The basic alg
  * iterates through allowing the bubbles to compete for space until
  * they converge (if the damping is dome properly) on a steady-state
- * solution.
- *
- * We use the ratio of local to average 'error' as the metric to bias a
- * variable-length word codebook, and probability of occurrence within
- * that bubble as the metric to bias fixed length word
- * codebooks. Individual input points, collected from libvorbis, are
+ * solution. Individual input points, collected from libvorbis, are
  * used to train the algorithm monte-carlo style.  */
 
 typedef struct vqgen{
@@ -54,29 +50,32 @@ typedef struct vqgen{
   long   allocated;
 
   /* entries */
-  double *entry_now;
-  double *entry_next;
+  double *entrylist;
   long   *assigned;
-  double *metric;
   double *bias;
   long   entries;
 
   double (*metric_func)   (struct vqgen *v,double *a,double *b);
 } vqgen;
 
+typedef struct vqbook{
+
+
+
+} vqbook;
+
 /* internal helpers *****************************************************/
+#define vN(data,i) (data+v->elements*i)
+
 double *_point(vqgen *v,long ptr){
   return v->pointlist+(v->elements*ptr);
 }
 
 double *_now(vqgen *v,long ptr){
-  return v->entry_now+(v->elements*ptr);
-}
-
-double *_next(vqgen *v,long ptr){
-  return v->entry_next+(v->elements*ptr);
+  return v->entrylist+(v->elements*ptr);
 }
 
+/* default metric; squared 'distance' from desired value. */
 double _dist_sq(vqgen *v,double *a, double *b){
   int i;
   int el=v->elements;
@@ -88,6 +87,46 @@ double _dist_sq(vqgen *v,double *a, double *b){
   return acc;
 }
 
+/* *must* be beefed up.  Perhaps a Floyd-Steinberg like scattering? */
+void _vqgen_seed(vqgen *v){
+  memcpy(v->entrylist,v->pointlist,sizeof(double)*v->entries*v->elements);
+}
+
+void _limited_sort(double *vals,long *index,long total,int desired){
+  if(total>1){
+    long bisect=total>>1;
+    if(bisect>1)_limited_sort(vals,index,bisect,desired);
+    if(total-bisect>1)_limited_sort(vals,index+bisect,total-bisect,
+                                   desired);
+    if(desired>total)desired=total;
+
+    {
+      int i;
+      long ptra=0;
+      long ptrb=bisect;
+      long *temp=alloca(desired*sizeof(long));
+
+      for(i=0;i<desired;i++){
+       if(vals[index[ptra]]>vals[index[ptrb]]){
+         temp[i]=index[ptra++];
+       }else{
+         temp[i]=index[ptrb++];
+       }
+       if(ptra==bisect)
+         for(i++;i<desired;i++)temp[i]=index[ptrb++];
+       else
+         if(ptrb==total)
+           for(i++;i<desired;i++)temp[i]=index[ptra++];
+          
+      }
+
+      memcpy(index,temp,desired*sizeof(long));
+    }
+  }
+}
+
+/* External calls *******************************************************/
+
 void vqgen_init(vqgen *v,int elements,int entries,
                double (*metric)(vqgen *,double *, double *),
                double spread){
@@ -101,10 +140,8 @@ void vqgen_init(vqgen *v,int elements,int entries,
   v->second=malloc(v->allocated*sizeof(long));
 
   v->entries=entries;
-  v->entry_now=malloc(v->entries*v->elements*sizeof(double));
-  v->entry_next=malloc(v->entries*v->elements*sizeof(double));
+  v->entrylist=malloc(v->entries*v->elements*sizeof(double));
   v->assigned=malloc(v->entries*sizeof(long));
-  v->metric=malloc(v->entries*sizeof(double));
   v->bias=calloc(v->entries,sizeof(double));
   if(metric)
     v->metric_func=metric;
@@ -112,11 +149,6 @@ void vqgen_init(vqgen *v,int elements,int entries,
     v->metric_func=_dist_sq;
 }
 
-/* *must* be beefed up.  Perhaps a Floyd-Steinberg like scattering? */
-void _vqgen_seed(vqgen *v){
-  memcpy(v->entry_now,v->pointlist,sizeof(double)*v->entries*v->elements);
-}
-
 void vqgen_addpoint(vqgen *v, double *p){
   if(v->points>=v->allocated){
     v->allocated*=2;
@@ -130,65 +162,99 @@ void vqgen_addpoint(vqgen *v, double *p){
   if(v->points==v->entries)_vqgen_seed(v);
 }
 
-double *sort_t_vals;
-int sort_t_compare(const void *ap,const void *bp){
-  long a=*(long *)ap;
-  long b=*(long *)bp;
-  double val=sort_t_vals[a]-sort_t_vals[b];
-  if(val<0)return(1);
-  if(val>0)return(-1);
-  return(0);
-}
+void vqgen_recenter(vqgen *v){
+  long   i,j,k;
+  double fdesired=(double)v->points/v->entries;
+  double asserror=0.;
+  double meterror=0.;
+
+  if(v->entries<2)exit(1);
 
-void vqgen_iterate(vqgen *v,int biasp){
-  static int iteration=0;
-  long i,j;
-  double avmetric=0.;
+  /* first round: new midpoints */
+  {
+    double *newlo=malloc(sizeof(double)*v->entries*v->elements);
+    double *newhi=malloc(sizeof(double)*v->entries*v->elements);
 
-  FILE *as;
-  FILE *graph;
-  FILE *err;
-  FILE *bias;
-  char name[80];
-  
+    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];
+      long   firstentry=0;
+    
+      for(j=1;j<v->entries;j++){
+       double thismetric=v->metric_func(v,_now(v,j),_point(v,i))+v->bias[j];
+       if(thismetric<firstmetric){
+         firstmetric=thismetric;
+         firstentry=j;
+       }
+      }
 
-  sprintf(name,"space%d.m",iteration);
-  graph=fopen(name,"w");
+      j=firstentry;
+      meterror+=firstmetric-v->bias[firstentry];
+      if(v->assigned[j]++){
+       for(k=0;k<v->elements;k++){
+         if(_point(v,i)[k]<vN(newlo,j)[k])vN(newlo,j)[k]=_point(v,i)[k];
+         if(_point(v,i)[k]>vN(newhi,j)[k])vN(newhi,j)[k]=_point(v,i)[k];
+       }
+      }else{
+       for(k=0;k<v->elements;k++){
+         vN(newlo,j)[k]=_point(v,i)[k];
+         vN(newhi,j)[k]=_point(v,i)[k];
+       }       
+      }
+    }
 
-  sprintf(name,"err%d.m",iteration);
-  err=fopen(name,"w");
+    for(j=0;j<v->entries;j++){
+      if(v->assigned[j]){
+       for(k=0;k<v->elements;k++)
+         _now(v,j)[k]=(vN(newlo,j)[k]+vN(newhi,j)[k])/2.;
+      }
+    }
+    free(newlo);
+    free(newhi);
+  }
 
-  sprintf(name,"bias%d.m",iteration);
-  bias=fopen(name,"w");
+  for(i=0;i<v->entries;i++){
+    asserror+=fabs(fdesired-v->assigned[i]);
+  }
 
-  sprintf(name,"as%d.m",iteration);
-  as=fopen(name,"w");
+  fprintf(stderr,"recenter: dist error=%g(%g) metric error=%g\n",
+         asserror/v->entries,fdesired,meterror/v->points);
+  
+  memset(v->bias,0,sizeof(double)*v->entries);
+}
 
+void vqgen_rebias(vqgen *v){
+  long   i,j,k;
+  double fdesired=(float)v->points/v->entries;
+  long   desired=(v->points+v->entries-1)/v->entries;
 
-  /* init */
-  memset(v->entry_next,0,sizeof(double)*v->elements*v->entries);
-  memset(v->assigned,0,sizeof(long)*v->entries);
-  memset(v->metric,0,sizeof(double)*v->entries);
+  double asserror=0.;
+  double meterror=0.;
+  long   *nearcount=calloc(v->entries,sizeof(long));
+  double *nearbias=malloc(v->entries*desired*sizeof(double));
 
   if(v->entries<2)exit(1);
 
-  /* assign all the points, accumulate metric */
+  /* second round: fill in nearest points for 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];
+    /* not clear this should be biased */
+    double firstmetric=v->metric_func(v,_now(v,0),_point(v,i));
+    double secondmetric=v->metric_func(v,_now(v,1),_point(v,i));
     long   firstentry=0;
-    long   secondentry=0;
-    
+    long   secondentry=1;
     if(firstmetric>secondmetric){
-      double tempmetric=firstmetric;
+      double temp=firstmetric;
       firstmetric=secondmetric;
-      secondmetric=tempmetric;
+      secondmetric=temp;
       firstentry=1;
-      secondentry=1;
+      secondentry=0;
     }
-
+      
     for(j=2;j<v->entries;j++){
-      double thismetric=v->metric_func(v,_now(v,j),_point(v,i))+v->bias[j];
+      double thismetric=v->metric_func(v,_now(v,j),_point(v,i));
       if(thismetric<secondmetric){
        if(thismetric<firstmetric){
          secondmetric=firstmetric;
@@ -201,122 +267,592 @@ void vqgen_iterate(vqgen *v,int biasp){
        }
       }
     }
-    
-    v->first[i]=firstentry;
-    v->second[i]=secondentry;
-    v->metric[firstentry]+=firstmetric-v->bias[firstentry];
     v->assigned[firstentry]++;
+    meterror+=firstmetric;
+    
+    for(j=0;j<v->entries;j++){
+      
+      double thismetric;
+      double *nearbiasptr=nearbias+desired*j;
+      long k=nearcount[j]-1;
+      
+      /* 'thismetric' is to be the bias value necessary in the current
+         arrangement for entry j to capture point i */
+      if(firstentry==j){
+       /* use the secondary entry as the threshhold */
+       thismetric=secondmetric-v->metric_func(v,_now(v,j),_point(v,i));
+      }else{
+       /* use the primary entry as the threshhold */
+       thismetric=firstmetric-v->metric_func(v,_now(v,j),_point(v,i));
+      }
 
-    {
-      double *no=_now(v,firstentry);
-      double *ne=_next(v,firstentry);
-      double *p=_point(v,i);
-      for(j=0;j<v->elements;j++)
-       ne[j]+=p[j];
-      fprintf(graph,"%g %g\n%g %g\n\n",p[0],p[1],no[0],no[1]);
+      if(k>=0 && thismetric>nearbiasptr[k]){
+       
+       /* start at the end and search backward for where this entry
+          belongs */
+       
+       for(;k>0;k--) if(nearbiasptr[k-1]>=thismetric)break;
+       
+       /* insert at k.  Shift and inject. */
+       memmove(nearbiasptr+k+1,nearbiasptr+k,(desired-k-1)*sizeof(double));
+       nearbiasptr[k]=thismetric;
+       
+       if(nearcount[j]<desired)nearcount[j]++;
+       
+      }else{
+       if(nearcount[j]<desired){
+         /* we checked the thresh earlier.  We know this is the
+            last entry */
+         nearbiasptr[nearcount[j]++]=thismetric;
+       }
+      }
+    }
+  }
+  
+  /* third round: inflate/deflate */
+  {
+    for(i=0;i<v->entries;i++){
+      v->bias[i]=nearbias[(i+1)*desired-1];
     }
   }
 
-  /* new midpoints */
   for(i=0;i<v->entries;i++){
-    double *next=_next(v,i);
-    double *now=_now(v,i);
-    if(v->assigned[i]){
-      for(j=0;j<v->elements;j++)
-       next[j]/=v->assigned[i];
-    }else{
-      for(j=0;j<v->elements;j++)
-       next[j]=now[j];
-    }
+    asserror+=fabs(fdesired-v->assigned[i]);
   }
+
+  fprintf(stderr,"rebias: dist error=%g(%g) metric error=%g\n",
+         asserror/v->entries,fdesired,meterror/v->points);
   
-  /* average global metric */
-  for(i=0;i<v->entries;i++){
-    avmetric+=v->metric[i];
+  free(nearcount);
+  free(nearbias);
+}
+
+/* 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;
+  long i,j,k;
+  long leaves=0;
+
+  double best=0.;
+  long besti=-1;
+  long bestj=-1;
+
+  if(depth>7){
+    printf("leaf: points %ld, depth %ld\n",points,depth);
+    return(1);
   }
-  avmetric/=v->entries;
-  fprintf(stderr,"%ld... ",iteration);
-
-  /* positive/negative 'pressure' */  
-  if(biasp){
-    /* Another round of n*m. Above we had to index by point.  Now, we
-       have to index by entry */
-    long   *index=alloca(v->points*sizeof(long));
-    double *vals =alloca(v->points*sizeof(double));
-    double *met  =alloca(v->points*sizeof(double));
-    for(i=0;i<v->entries;i++){
 
-      /* sort all n-thousand points by effective metric */
+  if(points==1){
+    printf("leaf: point %ld, depth %ld\n",index[0],depth);
+    return(1);
+  }
 
-      for(j=0;j<v->points;j++){
-       int threshentry;
-       index[j]=j;
-       if(v->first[j]==i){
-         /* use the secondary entry as the threshhold */
-         threshentry=v->second[j];
-       }else{
-         /* use the primary entry as the threshhold */
-         threshentry=v->first[j];
+  if(points==2){
+    /* The result must be an even split */
+    B[bp++]=index[0];
+    A[ap++]=index[1];
+    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){
+         besti=i;
+         bestj=j;
+         best=_dist_sq(v,_now(v,index[i]),_now(v,index[j]));
        }
-       
-       /* bias value right at the threshhold to grab this point */
-       met[j]=v->metric_func(v,_now(v,i),_point(v,j));
-       vals[j]=
-         v->metric_func(v,_now(v,threshentry),_point(v,j))+
-         v->bias[threshentry]-
-         v->metric_func(v,_now(v,i),_point(v,j));
       }
-      
-      sort_t_vals=vals;
-      qsort(index,v->points,sizeof(long),sort_t_compare);
+    }
 
-      {
-       /* find the desired number of points and level of bias.  We
-          nominally want all entries to represent the same number of
-          points, but we may also have a metric range restriction */
-       int desired=(v->points+v->entries-1)/v->entries;
+    /* 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;
+      }
+      
+      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(v->errspread>=1.){
-         double acc=0.;
-         for(j=0;j<desired;j++){
-           acc+=met[index[j]];
-           if(acc*v->errspread>avmetric)break;
+       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++;
+           }
          }
-       }else
-         j=desired;
-       v->bias[i]=vals[index[j-1]];
+         
+         /* 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;
+         }
+       }
+      }
+
+      /*{
+       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);
+       }
+
+       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);
+      
       }
+      free(cA);
+      free(cB);    
+      return(leaves);
     }
-  }else{
-    memset(v->bias,0,sizeof(double)*v->entries);
   }
+}
+
+void vqgen_book(vqgen *v){
+
 
-  /* dump state, report metric */
-  for(i=0;i<v->entries;i++){
-    for(j=0;j<v->assigned[i];j++)
-      fprintf(err,"%g\n",v->metric[i]/v->assigned[i]);
-    fprintf(bias,"%g\n",v->bias[i]);
-    fprintf(as,"%ld\n",v->assigned[i]);
-  }
 
-  fprintf(stderr,"average metric: %g\n",avmetric/v->elements/v->points*v->entries);
 
-  fclose(as);
-  fclose(err);
-  fclose(bias);
-  fclose(graph);
-  memcpy(v->entry_now,v->entry_next,sizeof(double)*v->entries*v->elements);
-  iteration++;
 }
 
+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,
+
+};
+
 int main(int argc,char *argv[]){
   FILE *in=fopen(argv[1],"r");
   vqgen v;
   char buffer[160];
-  int i;
-
-  vqgen_init(&v,4,24,_dist_sq,1.);
+  int i,j;
 
+  vqgen_init(&v,4,256,_dist_sq,0.);
+  
   while(fgets(buffer,160,in)){
     double a[8];
     if(sscanf(buffer,"%lf %lf %lf %lf",
@@ -324,17 +860,51 @@ int main(int argc,char *argv[]){
       vqgen_addpoint(&v,a);
   }
   fclose(in);
+  
+  for(i=0;i<10;i++)
+    vqgen_recenter(&v);
+  
+  for(i=0;i<100;i++){
+    vqgen_rebias(&v);
+    vqgen_recenter(&v);
+    fprintf(stderr,"%d ",i);
+  }
 
-  for(i=0;i<20;i++)
-    vqgen_iterate(&v,0);
+  vqgen_recenter(&v);
 
-  for(i=0;i<100;i++)
-    vqgen_iterate(&v,1);
+  exit(0);
 
-  for(i=0;i<5;i++)
-    vqgen_iterate(&v,0);
+  memcpy(v.entrylist,testset256,sizeof(testset256));
 
+  for(i=0;i<v.entries;i++){
+    printf("\n");
+    for(j=0;j<v.elements;j++)
+      printf("%f,",_now(&v,i)[j]);
+  }
+  printf("\n");
+    
+
+  /* try recursively splitting the space using LP 
+  {
+    long index[v.entries];
+    for(i=0;i<v.entries;i++)index[i]=i;
+    fprintf(stderr,"\n\nleaves=%d\n",
+           lp_split(&v,index,v.entries,NULL,NULL,0));
+           }*/
+  
   return(0);
 }
 
 
+
+
+
+
+
+
+
+
+
+
+
+