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