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