vq/vqgen.c

   1 /********************************************************************
   2  *                                                                  *
   3  * THIS FILE IS PART OF THE Ogg Vorbis SOFTWARE CODEC SOURCE CODE.  *
   4  * USE, DISTRIBUTION AND REPRODUCTION OF THIS SOURCE IS GOVERNED BY *
   5  * THE GNU PUBLIC LICENSE 2, WHICH IS INCLUDED WITH THIS SOURCE.    *
   6  * PLEASE READ THESE TERMS DISTRIBUTING.                            *
   7  *                                                                  *
   8  * THE OggSQUISH SOURCE CODE IS (C) COPYRIGHT 1994-1999             *
   9  * by 1999 Monty <monty@xiph.org> and The XIPHOPHORUS Company       *
  10  * http://www.xiph.org/                                             *
  11  *                                                                  *
  12  ********************************************************************
  13
  14  function: build a VQ codebook
  15  author: Monty <xiphmont@mit.edu>
  16  modifications by: Monty
  17  last modification date: Oct 08 1999
  18
  19  ********************************************************************/
  20
  21 /* This code is *not* part of libvorbis.  It is used to generate
  22    trained codebooks offline and then spit the results into a
  23    pregenerated codebook that is compiled into libvorbis.  It is an
  24    expensive (but good) algorithm.  Run it on big iron. */
  25
  26 #include <stdlib.h>
  27 #include <stdio.h>
  28 #include <math.h>
  29 #include <string.h>
  30 #include "minit.h"
  31
  32 /************************************************************************
  33  * The basic idea here is that a VQ codebook is like an m-dimensional
  34  * foam with n bubbles.  The bubbles compete for space/volume and are
  35  * 'pressurized' [biased] according to some metric.  The basic alg
  36  * iterates through allowing the bubbles to compete for space until
  37  * they converge (if the damping is dome properly) on a steady-state
  38  * solution. Individual input points, collected from libvorbis, are
  39  * used to train the algorithm monte-carlo style.  */
  40
  41 typedef struct vqgen{
  42   int    elements;
  43   double errspread;
  44
  45   /* point cache */
  46   double *pointlist;
  47   long   *first;
  48   long   *second;
  49   long   points;
  50   long   allocated;
  51
  52   /* entries */
  53   double *entrylist;
  54   long   *assigned;
  55   double *bias;
  56   long   entries;
  57
  58   double (*metric_func)   (struct vqgen *v,double *a,double *b);
  59 } vqgen;
  60
  61 typedef struct vqbook{
  62
  63
  64
  65 } vqbook;
  66
  67 /* internal helpers *****************************************************/
  68 #define vN(data,i) (data+v->elements*i)
  69
  70 double *_point(vqgen *v,long ptr){
  71   return v->pointlist+(v->elements*ptr);
  72 }
  73
  74 double *_now(vqgen *v,long ptr){
  75   return v->entrylist+(v->elements*ptr);
  76 }
  77
  78 /* default metric; squared 'distance' from desired value. */
  79 double _dist_sq(vqgen *v,double *a, double *b){
  80   int i;
  81   int el=v->elements;
  82   double acc=0.;
  83   for(i=0;i<el;i++){
  84     double val=(a[i]-b[i]);
  85     acc+=val*val;
  86   }
  87   return acc;
  88 }
  89                             /* candidate,actual */
  90 double _dist_and_pos(vqgen *v,double *b, double *a){
  91   int i;
  92   int el=v->elements;
  93   double acc=0.;
  94   double lastb=0.;
  95   for(i=0;i<el;i++){
  96     double actualdist=(a[i]-lastb);
  97     double testdist=(b[i]-lastb);
  98     if(actualdist>0 && testdist>0){
  99       double val;
 100       if(actualdist>testdist)
 101         val=actualdist/testdist;
 102       else
 103         val=testdist/actualdist;
 104       acc+=val*val-1.;
 105     }else{
 106       fprintf(stderr,"\nA zero (shouldn't happen in our data) \n");
 107       acc+=999999.;
 108     }
 109     lastb=b[i];
 110   }
 111   return acc;
 112 }
 113
 114 /* *must* be beefed up. */
 115 void _vqgen_seed(vqgen *v){
 116   memcpy(v->entrylist,v->pointlist,sizeof(double)*v->entries*v->elements);
 117 }
 118
 119 /* External calls *******************************************************/
 120
 121 void vqgen_init(vqgen *v,int elements,int entries,
 122                 double (*metric)(vqgen *,double *, double *),
 123                 double spread){
 124   memset(v,0,sizeof(vqgen));
 125
 126   v->elements=elements;
 127   v->errspread=spread;
 128   v->allocated=32768;
 129   v->pointlist=malloc(v->allocated*v->elements*sizeof(double));
 130   v->first=malloc(v->allocated*sizeof(long));
 131   v->second=malloc(v->allocated*sizeof(long));
 132
 133   v->entries=entries;
 134   v->entrylist=malloc(v->entries*v->elements*sizeof(double));
 135   v->assigned=malloc(v->entries*sizeof(long));
 136   v->bias=calloc(v->entries,sizeof(double));
 137   if(metric)
 138     v->metric_func=metric;
 139   else
 140     v->metric_func=_dist_and_pos;
 141 }
 142
 143 void vqgen_addpoint(vqgen *v, double *p){
 144   if(v->points>=v->allocated){
 145     v->allocated*=2;
 146     v->pointlist=realloc(v->pointlist,v->allocated*v->elements*sizeof(double));
 147     v->first=realloc(v->first,v->allocated*sizeof(long));
 148     v->second=realloc(v->second,v->allocated*sizeof(long));
 149   }
 150
 151   memcpy(_point(v,v->points),p,sizeof(double)*v->elements);
 152   v->points++;
 153   if(v->points==v->entries)_vqgen_seed(v);
 154 }
 155
 156 static long it=0;
 157 void vqgen_recenter(vqgen *v){
 158   long   i,j,k;
 159   double fdesired=(double)v->points/v->entries;
 160   double asserror=0.;
 161   double meterror=0.;
 162
 163 #ifdef NOISY
 164   char buff[80];
 165   FILE *cells;
 166   FILE *assig;
 167   FILE *bias;
 168   sprintf(buff,"cells%d.m",it);
 169   cells=fopen(buff,"w");
 170   sprintf(buff,"assig%d.m",it);
 171   assig=fopen(buff,"w");
 172   sprintf(buff,"bias%d.m",it);
 173   bias=fopen(buff,"w");
 174 #endif
 175
 176   if(v->entries<2)exit(1);
 177
 178   /* first round: new midpoints */
 179   {
 180     double *newlo=malloc(sizeof(double)*v->entries*v->elements);
 181     double *newhi=malloc(sizeof(double)*v->entries*v->elements);
 182
 183     memset(v->assigned,0,sizeof(long)*v->entries);
 184
 185     for(i=0;i<v->points;i++){
 186       double firstmetric=v->metric_func(v,_now(v,0),_point(v,i))+v->bias[0];
 187       long   firstentry=0;
 188
 189       for(j=1;j<v->entries;j++){
 190         double thismetric=v->metric_func(v,_now(v,j),_point(v,i))+v->bias[j];
 191         if(thismetric<firstmetric){
 192           firstmetric=thismetric;
 193           firstentry=j;
 194         }
 195       }
 196
 197       j=firstentry;
 198 #ifdef NOISY
 199       fprintf(cells,"%g %g\n%g %g\n\n",_now(v,j)[0],_now(v,j)[1],
 200               _point(v,i)[0],_point(v,i)[1]);
 201 #endif
 202       meterror+=firstmetric-v->bias[firstentry];
 203       if(v->assigned[j]++){
 204         for(k=0;k<v->elements;k++){
 205           vN(newlo,j)[k]+=_point(v,i)[k];
 206           /*if(_point(v,i)[k]<vN(newlo,j)[k])vN(newlo,j)[k]=_point(v,i)[k];
 207             if(_point(v,i)[k]>vN(newhi,j)[k])vN(newhi,j)[k]=_point(v,i)[k];*/
 208         }
 209       }else{
 210         for(k=0;k<v->elements;k++){
 211           vN(newlo,j)[k]=_point(v,i)[k];
 212           /*vN(newlo,j)[k]=_point(v,i)[k];
 213             vN(newhi,j)[k]=_point(v,i)[k];*/
 214         }
 215       }
 216     }
 217
 218     for(j=0;j<v->entries;j++){
 219       if(v->assigned[j]){
 220         for(k=0;k<v->elements;k++){
 221           _now(v,j)[k]=vN(newlo,j)[k]/v->assigned[j];
 222           /*_now(v,j)[k]=(vN(newlo,j)[k]+vN(newhi,j)[k])/2.;*/
 223         }
 224       }
 225     }
 226     free(newlo);
 227     free(newhi);
 228   }
 229
 230   for(i=0;i<v->entries;i++){
 231     asserror+=fabs(fdesired-v->assigned[i]);
 232 #ifdef NOISY
 233     fprintf(assig,"%ld\n",v->assigned[i]);
 234     fprintf(bias,"%g\n",v->bias[i]);
 235 #endif
 236   }
 237
 238   fprintf(stderr,"recenter: dist error=%g(%g) metric error=%g\n",
 239           asserror/v->entries,fdesired,meterror/v->points);
 240
 241   it++;
 242
 243 #ifdef NOISY
 244   fclose(cells);
 245   fclose(assig);
 246   fclose(bias);
 247 #endif
 248 }
 249
 250 void vqgen_rebias(vqgen *v){
 251   long   i,j,k;
 252   double fdesired=(float)v->points/v->entries;
 253   long   desired=fdesired;
 254
 255   double asserror=0.;
 256   double meterror=0.;
 257   long   *nearcount=calloc(v->entries,sizeof(long));
 258   double *nearbias=malloc(v->entries*desired*sizeof(double));
 259
 260   if(v->entries<2)exit(1);
 261
 262   /* second round: fill in nearest points for entries */
 263
 264   memset(v->assigned,0,sizeof(long)*v->entries);
 265   for(i=0;i<v->points;i++){
 266     /* not clear this should be biased */
 267     double firstmetric=v->metric_func(v,_now(v,0),_point(v,i))+v->bias[0];
 268     double secondmetric=v->metric_func(v,_now(v,1),_point(v,i))+v->bias[1];
 269     long   firstentry=0;
 270     long   secondentry=1;
 271     if(firstmetric>secondmetric){
 272       double temp=firstmetric;
 273       firstmetric=secondmetric;
 274       secondmetric=temp;
 275       firstentry=1;
 276       secondentry=0;
 277     }
 278
 279     for(j=2;j<v->entries;j++){
 280       double thismetric=v->metric_func(v,_now(v,j),_point(v,i))+v->bias[j];
 281       if(thismetric<secondmetric){
 282         if(thismetric<firstmetric){
 283           secondmetric=firstmetric;
 284           secondentry=firstentry;
 285           firstmetric=thismetric;
 286           firstentry=j;
 287         }else{
 288           secondmetric=thismetric;
 289           secondentry=j;
 290         }
 291       }
 292     }
 293     v->assigned[firstentry]++;
 294     meterror+=firstmetric-v->bias[firstentry];
 295
 296     for(j=0;j<v->entries;j++){
 297
 298       double thismetric;
 299       double *nearbiasptr=nearbias+desired*j;
 300       long k=nearcount[j]-1;
 301
 302       /* 'thismetric' is to be the bias value necessary in the current
 303          arrangement for entry j to capture point i */
 304       if(firstentry==j){
 305         /* use the secondary entry as the threshhold */
 306         thismetric=secondmetric-v->metric_func(v,_now(v,j),_point(v,i));
 307       }else{
 308         /* use the primary entry as the threshhold */
 309         thismetric=firstmetric-v->metric_func(v,_now(v,j),_point(v,i));
 310       }
 311
 312       if(k>=0 && thismetric>nearbiasptr[k]){
 313
 314         /* start at the end and search backward for where this entry
 315            belongs */
 316
 317         for(;k>0;k--) if(nearbiasptr[k-1]>=thismetric)break;
 318
 319         /* insert at k.  Shift and inject. */
 320         memmove(nearbiasptr+k+1,nearbiasptr+k,(desired-k-1)*sizeof(double));
 321         nearbiasptr[k]=thismetric;
 322
 323         if(nearcount[j]<desired)nearcount[j]++;
 324
 325       }else{
 326         if(nearcount[j]<desired){
 327           /* we checked the thresh earlier.  We know this is the
 328              last entry */
 329           nearbiasptr[nearcount[j]++]=thismetric;
 330         }
 331       }
 332     }
 333   }
 334
 335   /* third round: inflate/deflate */
 336   {
 337     for(i=0;i<v->entries;i++){
 338       v->bias[i]=nearbias[(i+1)*desired-1];
 339     }
 340   }
 341
 342   for(i=0;i<v->entries;i++){
 343     asserror+=fabs(fdesired-v->assigned[i]);
 344   }
 345
 346   fprintf(stderr,"rebias: dist error=%g(%g) metric error=%g\n",
 347           asserror/v->entries,fdesired,meterror/v->points);
 348
 349   free(nearcount);
 350   free(nearbias);
 351 }
 352
 353 /* the additional fields are the hyperplanes we've already used to
 354    split the set.  The additional hyperplanes are necessary to bound
 355    later splits.  One per depth */
 356
 357 int lp_split(vqgen *v,long *index,long points,
 358               double *additional_a, double *additional_b, long depth){
 359   static int frameno=0;
 360   long A[points];
 361   long B[points];
 362   long ap=0;
 363   long bp=0;
 364   long cp=0;
 365   long i,j,k;
 366   long leaves=0;
 367
 368   double best=0.;
 369   long besti=-1;
 370   long bestj=-1;
 371
 372   if(depth>7){
 373     printf("leaf: points %ld, depth %ld\n",points,depth);
 374     return(1);
 375   }
 376
 377   if(points==1){
 378     printf("leaf: point %ld, depth %ld\n",index[0],depth);
 379     return(1);
 380   }
 381
 382   if(points==2){
 383     /* The result must be an even split */
 384     B[bp++]=index[0];
 385     A[ap++]=index[1];
 386     printf("even split: depth %ld\n",depth);
 387     return(2);
 388   }else{
 389     /* We need to find the best split */
 390     for(i=0;i<points-1;i++){
 391       for(j=i+1;j<points;j++){
 392         if(_dist_sq(v,_now(v,index[i]),_now(v,index[j]))>best){
 393           besti=i;
 394           bestj=j;
 395           best=_dist_sq(v,_now(v,index[i]),_now(v,index[j]));
 396         }
 397       }
 398     }
 399
 400     /* We have our endpoints. initial divvy */
 401     {
 402       double *cA=malloc(v->elements*sizeof(double));
 403       double *cB=malloc(v->elements*sizeof(double));
 404       double **a=malloc((points+depth)*sizeof(double *));
 405       double *b=malloc((points+depth)*sizeof(double));
 406       double *aa=malloc((points+depth)*v->elements*sizeof(double));
 407       double dA=0.,dB=0.;
 408       minit l;
 409
 410       minit_init(&l,v->elements,points+depth-1,0);
 411
 412       /* divisor hyperplane */
 413       for(i=0;i<v->elements;i++){
 414         double m=(_now(v,index[besti])[i]+_now(v,index[bestj])[i])/2;
 415         cA[i]=_now(v,index[besti])[i]-_now(v,index[bestj])[i];
 416         cB[i]=-cA[i];
 417         dA+=cA[i]*m;
 418         dB-=cA[i]*m;
 419       }
 420
 421       for(i=0;i<points+depth;i++)
 422         a[i]=aa+i*v->elements;
 423
 424       /* check each bubble to see if it intersects the divisor hyperplane */
 425       for(j=0;j<points;j++){
 426         long count=0,m=0;
 427         double d1=0.,d2=0.;
 428         int ret;
 429
 430         if(j==besti){
 431           B[bp++]=index[j];
 432         }else if (j==bestj){
 433           A[ap++]=index[j];
 434         }else{
 435
 436           for(i=0;i<points;i++){
 437             if(i!=j){
 438               b[count]=0.;
 439               for(k=0;k<v->elements;k++){
 440                 double m=(_now(v,index[i])[k]+_now(v,index[j])[k])/2;
 441                 a[count][k]=_now(v,index[i])[k]-_now(v,index[j])[k];
 442                 b[count]+=a[count][k]*m;
 443               }
 444               count++;
 445             }
 446           }
 447
 448           /* additional bounding hyperplanes from previous splits */
 449           for(i=0;i<depth;i++){
 450             b[count]=0.;
 451             for(k=0;k<v->elements;k++)
 452               a[count][k]=additional_a[m++];
 453             b[count++]=additional_b[i];
 454           }
 455
 456           /* on what side of the dividing hyperplane is the current test
 457              point? */
 458           for(i=0;i<v->elements;i++)
 459             d1+=_now(v,index[j])[i]*cA[i];
 460
 461
 462           if(d1<dA){
 463             fprintf(stderr,"+");
 464             ret=minit_solve(&l,a,b,cA,1e-6,NULL,NULL,&d2);
 465           }else{
 466             fprintf(stderr,"-");
 467             ret=minit_solve(&l,a,b,cB,1e-6,NULL,NULL,&d2);
 468             d2= -d2;
 469           }
 470
 471           switch(ret){
 472           case 0:
 473             /* bounded solution */
 474             if(d1<dA){
 475               A[ap++]=index[j];
 476               if(d2>dA){cp++;B[bp++]=index[j];}
 477             }else{
 478               B[bp++]=index[j];
 479               if(d2<dA){cp++;A[ap++]=index[j];}
 480             }
 481             break;
 482           default:
 483             /* unbounded solution or no solution; we're in both sets */
 484             B[bp++]=index[j];
 485             A[ap++]=index[j];
 486             cp++;
 487             break;
 488           }
 489         }
 490       }
 491
 492       /*{
 493         char buffer[80];
 494         FILE *o;
 495         int i;
 496
 497         sprintf(buffer,"set%d.m",frameno);
 498         o=fopen(buffer,"w");
 499         for(i=0;i<points;i++)
 500           fprintf(o,"%g %g\n\n",_now(v,index[i])[0],_now(v,index[i])[1]);
 501         fclose(o);
 502
 503         sprintf(buffer,"setA%d.m",frameno);
 504         o=fopen(buffer,"w");
 505         for(i=0;i<ap;i++)
 506           fprintf(o,"%g %g\n\n",_now(v,A[i])[0],_now(v,A[i])[1]);
 507         fclose(o);
 508
 509         sprintf(buffer,"setB%d.m",frameno);
 510         o=fopen(buffer,"w");
 511         for(i=0;i<bp;i++)
 512           fprintf(o,"%g %g\n\n",_now(v,B[i])[0],_now(v,B[i])[1]);
 513         fclose(o);
 514
 515         sprintf(buffer,"div%d.m",frameno);
 516         o=fopen(buffer,"w");
 517         fprintf(o,"%g %g\n%g %g\n\n",
 518                 _now(v,index[besti])[0],
 519                 (dA-cA[0]*_now(v,index[besti])[0])/cA[1],
 520                 _now(v,index[bestj])[0],
 521                 (dA-cA[0]*_now(v,index[bestj])[0])/cA[1]);
 522         fclose(o);
 523
 524         sprintf(buffer,"bound%d.m",frameno);
 525         o=fopen(buffer,"w");
 526         for(i=0;i<depth;i++)
 527           fprintf(o,"%g %g\n%g %g\n\n",
 528                   _now(v,index[besti])[0],
 529                   (additional_b[i]-
 530                    additional_a[i*2]*_now(v,index[besti])[0])/
 531                   additional_a[i*2+1],
 532                   _now(v,index[bestj])[0],
 533                   (additional_b[i]-
 534                    additional_a[i*2]*_now(v,index[bestj])[0])/
 535                   additional_a[i*2+1]);
 536         fclose(o);
 537         frameno++;
 538         }*/
 539
 540
 541       minit_clear(&l);
 542       free(a);
 543       free(b);
 544       free(aa);
 545
 546       fprintf(stderr,"split: total=%ld depth=%ld set A=%ld:%ld:%ld=B\n",
 547               points,depth,ap-cp,cp,bp-cp);
 548
 549       /* handle adding this divisor hyperplane to the additional
 550          constraints.  The 'inside' side is different for A and B */
 551
 552       {
 553         /* marginally wasteful */
 554         double *newadd_a=alloca(sizeof(double)*v->elements*(depth+1));
 555         double *newadd_b=alloca(sizeof(double)*(depth+1));
 556         if(additional_a){
 557           memcpy(newadd_a,additional_a,sizeof(double)*v->elements*depth);
 558           memcpy(newadd_b,additional_b,sizeof(double)*depth);
 559         }
 560
 561         memcpy(newadd_a+v->elements*depth,cA,sizeof(double)*v->elements);
 562         newadd_b[depth]=dA;
 563         leaves=lp_split(v,A,ap,newadd_a,newadd_b,depth+1);
 564
 565         memcpy(newadd_a+v->elements*depth,cB,sizeof(double)*v->elements);
 566         newadd_b[depth]=dB;
 567         leaves+=lp_split(v,B,bp,newadd_a,newadd_b,depth+1);
 568
 569       }
 570       free(cA);
 571       free(cB);
 572       return(leaves);
 573     }
 574   }
 575 }
 576
 577 void vqgen_book(vqgen *v){
 578
 579
 580
 581
 582 }
 583
 584 static double testset24[48]={
 585 1.047758,1.245406,
 586 1.007972,1.150078,
 587 1.064390,1.146266,
 588 0.761842,0.826055,
 589 0.862073,1.243707,
 590 0.746351,1.055368,
 591 0.956844,1.048223,
 592 0.877782,1.111871,
 593 0.964564,1.309219,
 594 0.920062,1.168658,
 595 0.787654,0.895163,
 596 0.974003,1.200115,
 597 0.788360,1.142671,
 598 0.978414,1.122249,
 599 0.869938,0.954436,
 600 1.131220,1.348747,
 601 0.934630,1.108564,
 602 0.896666,1.045772,
 603 0.707360,0.920954,
 604 0.788914,0.970626,
 605 0.828682,1.043673,
 606 1.042016,1.190406,
 607 1.031643,1.068194,
 608 1.120824,1.220326,
 609 };
 610
 611 static double testset256[1024]={
 612 1.047758,1.245406,1.007972,1.150078,
 613 1.064390,1.146266,0.761842,0.826055,
 614 0.862073,1.243707,0.746351,1.055368,
 615 0.956844,1.048223,0.877782,1.111871,
 616 0.964564,1.309219,0.920062,1.168658,
 617 0.787654,0.895163,0.974003,1.200115,
 618 0.788360,1.142671,0.978414,1.122249,
 619 0.869938,0.954436,1.131220,1.348747,
 620 0.934630,1.108564,0.896666,1.045772,
 621 0.707360,0.920954,0.788914,0.970626,
 622 0.828682,1.043673,1.042016,1.190406,
 623 1.031643,1.068194,1.120824,1.220326,
 624 0.747084,0.888176,1.051535,1.221335,
 625 0.980283,1.087708,1.311154,1.445530,
 626 0.927548,1.107237,1.296902,1.480610,
 627 0.804415,0.921948,1.029952,1.231078,
 628 0.923413,1.064919,1.143109,1.242755,
 629 0.826403,1.010624,1.353139,1.455870,
 630 0.837930,1.073739,1.251770,1.405370,
 631 0.704437,0.890366,1.016676,1.136689,
 632 0.840924,1.151858,1.369346,1.558572,
 633 0.835229,1.095713,1.165585,1.311155,
 634 0.763157,1.000718,1.168780,1.345318,
 635 0.912852,1.040279,1.260220,1.468606,
 636 0.680839,0.920854,1.056199,1.214924,
 637 0.993572,1.124795,1.220880,1.416346,
 638 0.766878,0.994942,1.170035,1.432164,
 639 1.006875,1.069566,1.283594,1.531742,
 640 1.029033,1.188357,1.330626,1.519446,
 641 0.915581,1.121514,1.193260,1.431981,
 642 0.950939,1.117674,1.280092,1.487413,
 643 0.783569,0.863166,1.198348,1.309199,
 644 0.845209,0.949087,1.133299,1.253263,
 645 0.836299,1.050413,1.263712,1.482634,
 646 1.049458,1.194466,1.332964,1.485894,
 647 0.921516,1.086016,1.244810,1.349821,
 648 0.888983,1.111494,1.313032,1.351122,
 649 1.015785,1.142849,1.292746,1.463028,
 650 0.860002,1.249262,1.354167,1.534761,
 651 0.942136,1.053546,1.120840,1.373279,
 652 0.846661,1.059871,1.294814,1.504904,
 653 0.796592,1.061316,1.349680,1.494821,
 654 0.880504,1.033570,1.314379,1.484966,
 655 0.770438,1.039993,1.238670,1.412317,
 656 1.036503,1.068649,1.202522,1.487553,
 657 0.942615,1.056630,1.344355,1.562161,
 658 0.963759,1.096427,1.324085,1.497266,
 659 0.899956,1.206083,1.349655,1.572202,
 660 0.977023,1.132722,1.338437,1.488086,
 661 1.096643,1.242143,1.357079,1.513244,
 662 0.795193,1.011915,1.312761,1.419554,
 663 1.021869,1.075269,1.265230,1.475986,
 664 1.019546,1.186088,1.347798,1.567607,
 665 1.051191,1.172904,1.381952,1.511843,
 666 0.924332,1.192953,1.385378,1.420727,
 667 0.869539,0.962051,1.400803,1.539414,
 668 0.949875,1.053939,1.305695,1.422764,
 669 0.962223,1.085902,1.247702,1.462327,
 670 0.991777,1.182062,1.287421,1.508453,
 671 0.918274,1.107269,1.325467,1.454377,
 672 0.993589,1.184665,1.352294,1.475413,
 673 0.804960,1.011490,1.259031,1.462335,
 674 0.845949,1.173858,1.389681,1.441148,
 675 0.858168,1.019930,1.232463,1.440218,
 676 0.884411,1.206226,1.318968,1.457990,
 677 0.893469,1.074271,1.274410,1.433475,
 678 1.039654,1.061160,1.122426,1.189133,
 679 1.062243,1.177357,1.344765,1.542071,
 680 0.900964,1.116803,1.380397,1.419108,
 681 0.969535,1.216066,1.383582,1.435924,
 682 1.017913,1.330986,1.403452,1.450709,
 683 0.771534,0.973324,1.171604,1.236873,
 684 0.955683,1.029609,1.309284,1.498449,
 685 0.892165,1.192347,1.322444,1.356025,
 686 0.890209,1.023145,1.326649,1.431669,
 687 1.038788,1.188891,1.330844,1.436180,
 688 1.050355,1.177363,1.283447,1.479040,
 689 0.825367,0.960892,1.179047,1.381923,
 690 0.978931,1.122163,1.304581,1.532161,
 691 0.970750,1.149650,1.390225,1.547422,
 692 0.805416,0.945116,1.118815,1.343870,
 693 1.004106,1.185996,1.304086,1.464108,
 694 0.771591,1.055865,1.150908,1.228683,
 695 0.950175,1.104237,1.203437,1.387264,
 696 0.923504,1.122814,1.173989,1.318650,
 697 0.960635,1.042342,1.196995,1.395394,
 698 1.035111,1.064804,1.134387,1.341424,
 699 0.815182,1.083612,1.159518,1.255777,
 700 0.784079,0.830594,1.056005,1.289902,
 701 0.737774,0.994134,1.115470,1.264879,
 702 0.820162,1.105657,1.182947,1.423629,
 703 1.067552,1.142869,1.191626,1.387014,
 704 0.752416,0.918026,1.124637,1.414500,
 705 0.767963,0.840069,0.997625,1.325653,
 706 0.787595,0.865440,1.090071,1.227348,
 707 0.798877,1.105239,1.331379,1.440643,
 708 0.829079,1.133632,1.280774,1.470690,
 709 1.009195,1.132959,1.284044,1.500589,
 710 0.698569,0.824611,1.236682,1.462088,
 711 0.817460,0.985767,1.081910,1.257751,
 712 0.784033,0.882552,1.149678,1.326920,
 713 1.039403,1.085310,1.211033,1.558635,
 714 0.966795,1.168857,1.309960,1.497788,
 715 0.906244,1.027050,1.198846,1.323671,
 716 0.776495,1.185285,1.309167,1.380683,
 717 0.799628,0.927976,1.048238,1.299045,
 718 0.779808,0.881990,1.167898,1.220314,
 719 0.770446,0.918281,1.049189,1.179393,
 720 0.768416,1.037558,1.165760,1.302606,
 721 0.743121,0.814671,0.990501,1.224236,
 722 1.037050,1.068240,1.159690,1.426280,
 723 0.978810,1.214329,1.253336,1.324395,
 724 0.984003,1.121443,1.376382,1.510519,
 725 1.146510,1.229726,1.417616,1.781032,
 726 0.897163,1.147910,1.221186,1.371815,
 727 1.068525,1.211553,1.343551,1.506743,
 728 0.762313,1.091082,1.253251,1.472381,
 729 0.960562,1.041965,1.247053,1.399214,
 730 0.864482,1.123473,1.163412,1.238620,
 731 0.963484,1.132803,1.164992,1.250389,
 732 1.009456,1.139510,1.251339,1.449078,
 733 0.851837,1.113642,1.170290,1.362806,
 734 0.857073,0.962039,1.127381,1.471682,
 735 1.047754,1.213381,1.388899,1.492383,
 736 0.938921,1.267308,1.337076,1.478427,
 737 0.790388,0.912816,1.159450,1.273259,
 738 0.832690,0.997776,1.156639,1.302621,
 739 0.783009,0.975374,1.080630,1.311112,
 740 0.819784,1.145093,1.326949,1.525480,
 741 0.806394,1.089564,1.329564,1.550362,
 742 0.958608,1.036364,1.379118,1.443043,
 743 0.753680,0.941781,1.147749,1.297211,
 744 0.883974,1.091394,1.315093,1.480307,
 745 1.239591,1.417601,1.495843,1.641941,
 746 0.881068,0.973780,1.278918,1.429384,
 747 0.869052,0.977661,1.280744,1.551295,
 748 1.022685,1.052986,1.105046,1.168670,
 749 0.981698,1.131448,1.197781,1.467704,
 750 0.945034,1.163410,1.250872,1.428793,
 751 1.092055,1.139380,1.187113,1.264586,
 752 0.788897,0.953764,1.232844,1.506461,
 753 0.885206,0.978419,1.209467,1.387569,
 754 0.762835,0.964279,1.064844,1.243153,
 755 0.974906,1.134763,1.283544,1.460386,
 756 1.081114,1.169260,1.290987,1.535452,
 757 0.880531,1.075660,1.378980,1.522978,
 758 1.039431,1.083792,1.350481,1.588522,
 759 0.922449,1.060622,1.152309,1.467930,
 760 0.918935,1.073732,1.283243,1.479567,
 761 0.992722,1.145398,1.239500,1.509004,
 762 0.903405,1.243062,1.421386,1.572148,
 763 0.912531,1.100239,1.310219,1.521424,
 764 0.875168,0.912787,1.123618,1.363635,
 765 0.732241,1.108317,1.323119,1.515675,
 766 0.951518,1.141874,1.341623,1.533409,
 767 0.992099,1.215801,1.477413,1.734691,
 768 1.005267,1.122655,1.356975,1.436445,
 769 0.740659,0.807749,1.167728,1.380618,
 770 1.078727,1.214706,1.328173,1.556699,
 771 1.026027,1.168906,1.313249,1.486078,
 772 0.914877,1.147150,1.389489,1.523984,
 773 1.062339,1.120684,1.160024,1.234794,
 774 1.129141,1.197655,1.374217,1.547755,
 775 1.070011,1.255271,1.360225,1.467869,
 776 0.779980,0.871696,1.098031,1.284490,
 777 1.062409,1.177542,1.314581,1.696662,
 778 0.702935,1.229873,1.370813,1.479500,
 779 1.029357,1.225167,1.341607,1.478163,
 780 1.025666,1.141749,1.185959,1.332892,
 781 0.799462,0.951470,1.214070,1.305787,
 782 0.740521,0.805457,1.107504,1.317258,
 783 0.784194,0.838683,1.055934,1.242692,
 784 0.839416,1.048060,1.391801,1.623786,
 785 0.692627,0.907677,1.060843,1.341002,
 786 0.823625,1.244497,1.396901,1.586243,
 787 0.942859,1.232978,1.348170,1.536735,
 788 0.894882,1.131376,1.292892,1.462724,
 789 0.776974,0.904449,1.325557,1.451968,
 790 1.066188,1.218328,1.376282,1.545381,
 791 0.990053,1.124279,1.340534,1.559666,
 792 0.776321,0.935169,1.191081,1.372326,
 793 0.949935,1.115929,1.397704,1.442549,
 794 0.936357,1.126885,1.356247,1.579931,
 795 1.045433,1.274605,1.366947,1.590215,
 796 1.063890,1.138062,1.220645,1.460005,
 797 0.751448,0.811338,1.078027,1.403146,
 798 0.935678,1.102858,1.356557,1.515948,
 799 1.026417,1.143843,1.299309,1.413976,
 800 0.821475,1.000237,1.105073,1.379882,
 801 0.960249,1.031602,1.250804,1.462620,
 802 0.661745,1.106452,1.291188,1.439529,
 803 0.691661,0.947241,1.261183,1.457391,
 804 0.839024,0.914750,1.040695,1.375853,
 805 1.029401,1.065349,1.121370,1.270670,
 806 0.888316,1.003349,1.103749,1.341290,
 807 0.766328,0.879083,1.217779,1.419868,
 808 0.811672,1.049673,1.186460,1.375742,
 809 0.969585,1.170126,1.259338,1.595070,
 810 1.016617,1.145706,1.335413,1.503064,
 811 0.980227,1.295316,1.417964,1.478684,
 812 0.885701,1.248099,1.416821,1.465338,
 813 0.953583,1.266596,1.325829,1.372230,
 814 1.038619,1.225860,1.351664,1.527048,
 815 1.104724,1.215839,1.250491,1.335237,
 816 1.124449,1.269485,1.420756,1.677260,
 817 0.881337,1.352259,1.433218,1.492756,
 818 1.023097,1.059205,1.110249,1.212976,
 819 1.135632,1.282243,1.415540,1.566039,
 820 1.063524,1.252342,1.399082,1.518433,
 821 0.790658,0.856337,1.154909,1.274963,
 822 1.119059,1.382625,1.423651,1.492741,
 823 1.030526,1.296610,1.330511,1.396550,
 824 0.947952,1.065235,1.225274,1.554455,
 825 0.960669,1.282313,1.370362,1.572736,
 826 0.996042,1.208193,1.385036,1.534877,
 827 1.003206,1.377432,1.431110,1.497189,
 828 1.088166,1.227944,1.508129,1.740407,
 829 0.996566,1.162407,1.347665,1.524235,
 830 0.944606,1.287026,1.400822,1.437156,
 831 1.066144,1.238314,1.454451,1.616016,
 832 1.121065,1.200875,1.316542,1.459583,
 833 1.158944,1.271448,1.356823,1.450510,
 834 0.811670,1.278011,1.361550,1.440647,
 835 0.875620,1.103051,1.230854,1.483429,
 836 0.959882,1.180685,1.381224,1.572807,
 837 1.049222,1.186513,1.387883,1.567423,
 838 1.049920,1.293154,1.507194,1.678371,
 839 0.872138,1.193727,1.455817,1.665837,
 840 0.738018,0.946583,1.335543,1.552061,
 841 1.072856,1.151838,1.321877,1.486028,
 842 1.026909,1.153575,1.306700,1.550408,
 843 0.940629,1.121943,1.303228,1.409285,
 844 1.023544,1.269480,1.316695,1.508732,
 845 0.924872,1.119626,1.479602,1.698360,
 846 0.793358,1.160361,1.381276,1.636124,
 847 0.892504,1.248046,1.320256,1.367114,
 848 1.165222,1.230265,1.507234,1.820748,
 849 0.711084,0.876165,1.145850,1.493303,
 850 0.776482,0.869670,1.124462,1.370831,
 851 0.785085,0.912534,1.099328,1.281523,
 852 0.947788,1.168239,1.407409,1.473592,
 853 0.913992,1.316132,1.619366,1.865030,
 854 1.024592,1.151907,1.383557,1.545112,
 855 0.987748,1.168194,1.415323,1.639563,
 856 0.725008,1.037088,1.306659,1.474607,
 857 0.965243,1.263523,1.454521,1.638929,
 858 1.035169,1.219404,1.427315,1.578414,
 859 0.787215,0.960639,1.245517,1.423549,
 860 0.902895,1.153039,1.412772,1.613360,
 861 0.983205,1.338886,1.493960,1.619623,
 862 0.822576,1.010524,1.232397,1.359007,
 863 0.773450,1.090005,1.170248,1.369534,
 864 1.117781,1.354485,1.662228,1.829565,
 865 1.104316,1.324945,1.525085,1.741040,
 866 1.135275,1.371208,1.537082,1.673949,
 867 0.889899,1.195451,1.242648,1.291881,
 868
 869 };
 870
 871 int main(int argc,char *argv[]){
 872   FILE *in=fopen(argv[1],"r");
 873   vqgen v;
 874   char buffer[160];
 875   int i,j;
 876
 877   vqgen_init(&v,4,256,_dist_sq,0.);
 878
 879   while(fgets(buffer,160,in)){
 880     double a[8];
 881     if(sscanf(buffer,"%lf %lf %lf %lf",
 882               a,a+1,a+2,a+3)==4)
 883       vqgen_addpoint(&v,a);
 884   }
 885   fclose(in);
 886
 887   for(i=0;i<10;i++)
 888     vqgen_recenter(&v);
 889
 890   for(i=0;i<200;i++){
 891     vqgen_rebias(&v);
 892     vqgen_rebias(&v);
 893     vqgen_rebias(&v);
 894     vqgen_rebias(&v);
 895     vqgen_rebias(&v);
 896     vqgen_recenter(&v);
 897     fprintf(stderr,"%d ",i);
 898   }
 899
 900   vqgen_recenter(&v);
 901
 902   exit(0);
 903
 904   memcpy(v.entrylist,testset256,sizeof(testset256));
 905
 906   for(i=0;i<v.entries;i++){
 907     printf("\n");
 908     for(j=0;j<v.elements;j++)
 909       printf("%f,",_now(&v,i)[j]);
 910   }
 911   printf("\n");
 912
 913
 914   /* try recursively splitting the space using LP
 915   {
 916     long index[v.entries];
 917     for(i=0;i<v.entries;i++)index[i]=i;
 918     fprintf(stderr,"\n\nleaves=%d\n",
 919             lp_split(&v,index,v.entries,NULL,NULL,0));
 920             }*/
 921
 922   return(0);
 923 }
 924
 925
 926
 927
 928
 929
 930
 931
 932
 933
 934
 935
 936
 937