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: Dec 09 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 /* There are so many optimizations to explore in *both* stages that
  27    considering the undertaking is almost withering.  For now, we brute
  28    force it all */
  29
  30 #include <stdlib.h>
  31 #include <stdio.h>
  32 #include <math.h>
  33 #include <string.h>
  34 #include "minit.h"
  35
  36 /* Codebook generation happens in two steps:
  37
  38    1) Train the codebook with data collected from the encoder: We use
  39    one of a few error metrics (which represent the distance between a
  40    given data point and a candidate point in the training set) to
  41    divide the training set up into cells representing roughly equal
  42    probability of occurring.
  43
  44    2) Generate the codebook and auxiliary data from the trained data set
  45 */
  46
  47 /* Codebook training ****************************************************
  48  *
  49  * The basic idea here is that a VQ codebook is like an m-dimensional
  50  * foam with n bubbles.  The bubbles compete for space/volume and are
  51  * 'pressurized' [biased] according to some metric.  The basic alg
  52  * iterates through allowing the bubbles to compete for space until
  53  * they converge (if the damping is dome properly) on a steady-state
  54  * solution. Individual input points, collected from libvorbis, are
  55  * used to train the algorithm monte-carlo style.  */
  56
  57 typedef struct vqgen{
  58   int it;
  59
  60   int    elements;
  61   double errspread;
  62
  63   /* point cache */
  64   double *pointlist;
  65   long   *first;
  66   long   *second;
  67   long   points;
  68   long   allocated;
  69
  70   /* entries */
  71   double *entrylist;
  72   long   *assigned;
  73   double *bias;
  74   long   entries;
  75
  76   double (*metric_func)   (struct vqgen *v,double *a,double *b);
  77 } vqgen;
  78
  79 /* internal helpers *****************************************************/
  80 #define vN(data,i) (data+v->elements*i)
  81
  82 double *_point(vqgen *v,long ptr){
  83   return v->pointlist+(v->elements*ptr);
  84 }
  85
  86 double *_now(vqgen *v,long ptr){
  87   return v->entrylist+(v->elements*ptr);
  88 }
  89
  90 /* default metric; squared 'distance' from desired value. */
  91 double _dist_sq(vqgen *v,double *a, double *b){
  92   int i;
  93   int el=v->elements;
  94   double acc=0.;
  95   for(i=0;i<el;i++){
  96     double val=(a[i]-b[i]);
  97     acc+=val*val;
  98   }
  99   return acc;
 100 }
 101
 102 /* A metric for LSP codes */
 103                             /* candidate,actual */
 104 double _dist_and_pos(vqgen *v,double *b, double *a){
 105   int i;
 106   int el=v->elements;
 107   double acc=0.;
 108   double lastb=0.;
 109   for(i=0;i<el;i++){
 110     double actualdist=(a[i]-lastb);
 111     double testdist=(b[i]-lastb);
 112     if(actualdist>0 && testdist>0){
 113       double val;
 114       if(actualdist>testdist)
 115         val=actualdist/testdist-1.;
 116       else
 117         val=testdist/actualdist-1.;
 118       acc+=val;
 119     }else{
 120       acc+=999999.;
 121     }
 122     lastb=b[i];
 123   }
 124   return acc;
 125 }
 126
 127 /* *must* be beefed up. */
 128 void _vqgen_seed(vqgen *v){
 129   memcpy(v->entrylist,v->pointlist,sizeof(double)*v->entries*v->elements);
 130 }
 131
 132 /* External calls *******************************************************/
 133
 134 void vqgen_init(vqgen *v,int elements,int entries,
 135                 double (*metric)(vqgen *,double *, double *),
 136                 double spread){
 137   memset(v,0,sizeof(vqgen));
 138
 139   v->elements=elements;
 140   v->errspread=spread;
 141   v->allocated=32768;
 142   v->pointlist=malloc(v->allocated*v->elements*sizeof(double));
 143   v->first=malloc(v->allocated*sizeof(long));
 144   v->second=malloc(v->allocated*sizeof(long));
 145
 146   v->entries=entries;
 147   v->entrylist=malloc(v->entries*v->elements*sizeof(double));
 148   v->assigned=malloc(v->entries*sizeof(long));
 149   v->bias=calloc(v->entries,sizeof(double));
 150   if(metric)
 151     v->metric_func=metric;
 152   else
 153     v->metric_func=_dist_and_pos;
 154 }
 155
 156 void vqgen_addpoint(vqgen *v, double *p){
 157   if(v->points>=v->allocated){
 158     v->allocated*=2;
 159     v->pointlist=realloc(v->pointlist,v->allocated*v->elements*sizeof(double));
 160     v->first=realloc(v->first,v->allocated*sizeof(long));
 161     v->second=realloc(v->second,v->allocated*sizeof(long));
 162   }
 163
 164   memcpy(_point(v,v->points),p,sizeof(double)*v->elements);
 165   v->points++;
 166   if(v->points==v->entries)_vqgen_seed(v);
 167 }
 168
 169 double vqgen_iterate(vqgen *v){
 170   long   i,j,k;
 171   double fdesired=(double)v->points/v->entries;
 172   long  desired=fdesired;
 173   double asserror=0.;
 174   double meterror=0.;
 175   double *new=malloc(sizeof(double)*v->entries*v->elements);
 176   long   *nearcount=malloc(v->entries*sizeof(long));
 177   double *nearbias=malloc(v->entries*desired*sizeof(double));
 178
 179 #ifdef NOISY
 180   char buff[80];
 181   FILE *assig;
 182   FILE *bias;
 183   FILE *cells;
 184   sprintf(buff,"cells%d.m",v->it);
 185   cells=fopen(buff,"w");
 186   sprintf(buff,"assig%d.m",v->it);
 187   assig=fopen(buff,"w");
 188   sprintf(buff,"bias%d.m",v->it);
 189   bias=fopen(buff,"w");
 190 #endif
 191
 192   fprintf(stderr,"Pass #%d... ",v->it);
 193
 194   if(v->entries<2){
 195     fprintf(stderr,"generation requires at least two entries\n");
 196     exit(1);
 197   }
 198
 199   /* fill in nearest points for entries */
 200   /*memset(v->bias,0,sizeof(double)*v->entries);*/
 201   memset(nearcount,0,sizeof(long)*v->entries);
 202   memset(v->assigned,0,sizeof(long)*v->entries);
 203   for(i=0;i<v->points;i++){
 204     double *ppt=_point(v,i);
 205     double firstmetric=v->metric_func(v,_now(v,0),ppt)+v->bias[0];
 206     double secondmetric=v->metric_func(v,_now(v,1),ppt)+v->bias[1];
 207     long   firstentry=0;
 208     long   secondentry=1;
 209     if(firstmetric>secondmetric){
 210       double temp=firstmetric;
 211       firstmetric=secondmetric;
 212       secondmetric=temp;
 213       firstentry=1;
 214       secondentry=0;
 215     }
 216
 217     for(j=2;j<v->entries;j++){
 218       double thismetric=v->metric_func(v,_now(v,j),_point(v,i))+v->bias[j];
 219       if(thismetric<secondmetric){
 220         if(thismetric<firstmetric){
 221           secondmetric=firstmetric;
 222           secondentry=firstentry;
 223           firstmetric=thismetric;
 224           firstentry=j;
 225         }else{
 226           secondmetric=thismetric;
 227           secondentry=j;
 228         }
 229       }
 230     }
 231
 232     j=firstentry;
 233     meterror+=firstmetric-v->bias[firstentry];
 234     /* set up midpoints for next iter */
 235     if(v->assigned[j]++)
 236       for(k=0;k<v->elements;k++)
 237         vN(new,j)[k]+=_point(v,i)[k];
 238     else
 239       for(k=0;k<v->elements;k++)
 240         vN(new,j)[k]=_point(v,i)[k];
 241
 242
 243 #ifdef NOISY
 244     fprintf(cells,"%g %g\n%g %g\n\n",
 245             _now(v,j)[0],_now(v,j)[1],
 246             _point(v,i)[0],_point(v,i)[1]);
 247 #endif
 248
 249     for(j=0;j<v->entries;j++){
 250
 251       double thismetric;
 252       double *nearbiasptr=nearbias+desired*j;
 253       long k=nearcount[j]-1;
 254
 255       /* 'thismetric' is to be the bias value necessary in the current
 256          arrangement for entry j to capture point i */
 257       if(firstentry==j){
 258         /* use the secondary entry as the threshhold */
 259         thismetric=secondmetric-v->metric_func(v,_now(v,j),_point(v,i));
 260       }else{
 261         /* use the primary entry as the threshhold */
 262         thismetric=firstmetric-v->metric_func(v,_now(v,j),_point(v,i));
 263       }
 264
 265       if(k>=0 && thismetric>nearbiasptr[k]){
 266
 267         /* start at the end and search backward for where this entry
 268            belongs */
 269
 270         for(;k>0;k--) if(nearbiasptr[k-1]>=thismetric)break;
 271
 272         /* insert at k.  Shift and inject. */
 273         memmove(nearbiasptr+k+1,nearbiasptr+k,(desired-k-1)*sizeof(double));
 274         nearbiasptr[k]=thismetric;
 275
 276         if(nearcount[j]<desired)nearcount[j]++;
 277
 278       }else{
 279         if(nearcount[j]<desired){
 280           /* we checked the thresh earlier.  We know this is the
 281              last entry */
 282           nearbiasptr[nearcount[j]++]=thismetric;
 283         }
 284       }
 285     }
 286   }
 287
 288   /* inflate/deflate */
 289   for(i=0;i<v->entries;i++)
 290     v->bias[i]=nearbias[(i+1)*desired-1];
 291
 292   /* last, assign midpoints */
 293   for(j=0;j<v->entries;j++){
 294     asserror+=fabs(v->assigned[j]-fdesired);
 295     if(v->assigned[j])
 296       for(k=0;k<v->elements;k++)
 297         _now(v,j)[k]=vN(new,j)[k]/v->assigned[j];
 298 #ifdef NOISY
 299     fprintf(assig,"%ld\n",v->assigned[j]);
 300     fprintf(bias,"%g\n",v->bias[j]);
 301 #endif
 302   }
 303
 304   fprintf(stderr,": dist %g(%g) metric error=%g \n",
 305           asserror/v->entries,fdesired,meterror/v->points);
 306   v->it++;
 307
 308   free(new);
 309   free(nearcount);
 310   free(nearbias);
 311 #ifdef NOISY
 312   fclose(assig);
 313   fclose(bias);
 314   fclose(cells);
 315 #endif
 316   return(asserror);
 317 }
 318
 319 /* Building a codebook from trained set **********************************
 320
 321    The codebook in raw form is technically finished once it's trained.
 322    However, we want to finalize the representative codebook values for
 323    each entry and generate auxiliary information to optimize encoding.
 324    We generate the auxiliary coding tree using collected data,
 325    probably the same data as in the original training */
 326
 327 /* At each recursion, the data set is split in half.  Cells with data
 328    points on side A go into set A, same with set B.  The sets may
 329    overlap.  If the cell overlaps the deviding line only very slightly
 330    (provided parameter), we may choose to ignore the overlap in order
 331    to pare the tree down */
 332
 333 typedef struct vqbook{
 334   long elements;
 335   long entries;
 336   double *valuelist;
 337   long   *codelist;
 338   long   *lengthlist;
 339
 340   /* auxiliary encoding/decoding information */
 341   long   *ptr0;
 342   long   *ptr1;
 343
 344   /* auxiliary encoding information */
 345   double *n;
 346   double *c;
 347   long   aux;
 348   long   alloc;
 349
 350 } vqbook;
 351
 352
 353 double *sortvals;
 354 int els;
 355 int iascsort(const void *a,const void *b){
 356   double av=sortvals[*((long *)a) * els];
 357   double bv=sortvals[*((long *)b) * els];
 358   if(av<bv)return(-1);
 359   return(1);
 360 }
 361
 362 /* goes through the split, but just counts it and returns a metric*/
 363 void lp_count(vqgen *v,long *entryindex,long entries,
 364                 long *pointindex,long points,
 365                 long *entryA,long *entryB,
 366                 double *n, double c, double slack,
 367                 long *entriesA,long *entriesB,long *entriesC){
 368   long i,j,k;
 369   long A=0,B=0,C=0;
 370
 371   memset(entryA,0,sizeof(long)*entries);
 372   memset(entryB,0,sizeof(long)*entries);
 373
 374   for(i=0;i<points;i++){
 375     double *ppt=_point(v,pointindex[i]);
 376     long   firstentry=0;
 377     double firstmetric=_dist_sq(v,_now(v,entryindex[0]),ppt);
 378     double position=-c;
 379
 380     for(j=1;j<entries;j++){
 381       double thismetric=_dist_sq(v,_now(v,entryindex[j]),ppt);
 382       if(thismetric<firstmetric){
 383         firstmetric=thismetric;
 384         firstentry=j;
 385       }
 386     }
 387
 388     /* count point split */
 389     for(k=0;k<v->elements;k++)
 390       position+=ppt[k]*n[k];
 391     if(position>0.){
 392       entryA[firstentry]++;
 393     }else{
 394       entryB[firstentry]++;
 395     }
 396   }
 397
 398   /* look to see if entries are in the slack zone */
 399   /* The entry splitting isn't total, so that storage has to be
 400      allocated for recursion.  Reuse the entryA/entryB vectors */
 401   for(j=0;j<entries;j++){
 402     long total=entryA[j]+entryB[j];
 403     if((double)entryA[j]/total<slack){
 404       entryA[j]=0;
 405     }else if((double)entryB[j]/total<slack){
 406       entryB[j]=0;
 407     }
 408     if(entryA[j] && entryB[j])C++;
 409     if(entryA[j])entryA[A++]=entryindex[j];
 410     if(entryB[j])entryB[B++]=entryindex[j];
 411   }
 412   *entriesA=A;
 413   *entriesB=B;
 414   *entriesC=C;
 415 }
 416
 417 void pq_in_out(vqgen *v,double *n,double *c,double *p,double *q){
 418   int k;
 419   *c=0.;
 420   for(k=0;k<v->elements;k++){
 421     double center=(p[k]+q[k])/2.;
 422     n[k]=(center-q[k])*2.;
 423     *c+=center*n[k];
 424   }
 425 }
 426
 427 void pq_center_out(vqgen *v,double *n,double *c,double *center,double *q){
 428   int k;
 429   *c=0.;
 430   for(k=0;k<v->elements;k++){
 431     n[k]=(center[k]-q[k])*2.;
 432     *c+=center[k]*n[k];
 433   }
 434 }
 435
 436 int lp_split(vqgen *v,vqbook *b,
 437              long *entryindex,long entries,
 438              long *pointindex,long points,
 439              long depth,double slack){
 440
 441   /* The encoder, regardless of book, will be using a straight
 442      euclidian distance-to-point metric to determine closest point.
 443      Thus we split the cells using the same (we've already trained the
 444      codebook set spacing and distribution using special metrics and
 445      even a midpoint division won't disturb the basic properties) */
 446
 447   long ret;
 448   double *p;
 449   double *q;
 450   double *n;
 451   double c;
 452   long *entryA=calloc(entries,sizeof(long));
 453   long *entryB=calloc(entries,sizeof(long));
 454   long entriesA=0;
 455   long entriesB=0;
 456   long entriesC=0;
 457   long pointsA=0;
 458   long i,j,k;
 459
 460   p=alloca(sizeof(double)*v->elements);
 461   q=alloca(sizeof(double)*v->elements);
 462   n=alloca(sizeof(double)*v->elements);
 463   memset(p,0,sizeof(double)*v->elements);
 464
 465   /* We need to find the dividing hyperplane. find the median of each
 466      axis as the centerpoint and the normal facing farthest point */
 467
 468   /* more than one way to do this part.  For small sets, we can brute
 469      force it. */
 470
 471   if(entries<32){
 472     /* try every pair possibility */
 473     double best=0;
 474     long   besti=0;
 475     long   bestj=0;
 476     double this;
 477     for(i=0;i<entries-1;i++){
 478       for(j=i+1;j<entries;j++){
 479         pq_in_out(v,n,&c,_now(v,entryindex[i]),_now(v,entryindex[j]));
 480         lp_count(v,entryindex,entries,
 481                  pointindex,points,
 482                  entryA,entryB,
 483                  n, c, slack,
 484                  &entriesA,&entriesB,&entriesC);
 485         this=(entriesA-entriesC)*(entriesB-entriesC);
 486
 487         if(this>best){
 488           best=this;
 489           besti=i;
 490           bestj=j;
 491         }
 492       }
 493     }
 494     pq_in_out(v,n,&c,_now(v,entryindex[besti]),_now(v,entryindex[bestj]));
 495   }else{
 496     double best=0.;
 497     long   bestj=0;
 498
 499     /* try COG/normal and furthest pairs */
 500     /* medianpoint */
 501     for(k=0;k<v->elements;k++){
 502       /* just sort the index array */
 503       sortvals=v->pointlist+k;
 504       els=v->elements;
 505       qsort(pointindex,points,sizeof(long),iascsort);
 506       if(points&0x1){
 507         p[k]=v->pointlist[(pointindex[points/2])*v->elements+k];
 508       }else{
 509         p[k]=(v->pointlist[(pointindex[points/2])*v->elements+k]+
 510               v->pointlist[(pointindex[points/2-1])*v->elements+k])/2.;
 511       }
 512     }
 513
 514     /* try every normal, but just for distance */
 515     for(j=0;j<entries;j++){
 516       double *ppj=_now(v,entryindex[j]);
 517       double this=_dist_sq(v,p,ppj);
 518       if(this>best){
 519         best=this;
 520         bestj=j;
 521       }
 522     }
 523
 524     pq_center_out(v,n,&c,p,_point(v,pointindex[bestj]));
 525
 526
 527   }
 528
 529   /* find cells enclosing points */
 530   /* count A/B points */
 531
 532   lp_count(v,entryindex,entries,
 533            pointindex,points,
 534            entryA,entryB,
 535            n, c, slack,
 536            &entriesA,&entriesB,&entriesC);
 537
 538   /* the point index is split evenly, so we do an Order n
 539      rearrangement into A first/B last and just pass it on */
 540   {
 541     long Aptr=0;
 542     long Bptr=points-1;
 543     while(Aptr<=Bptr){
 544       while(Aptr<=Bptr){
 545         double position=-c;
 546         for(k=0;k<v->elements;k++)
 547           position+=_point(v,pointindex[Aptr])[k]*n[k];
 548         if(position<=0.)break; /* not in A */
 549         Aptr++;
 550       }
 551       while(Aptr<=Bptr){
 552         double position=-c;
 553         for(k=0;k<v->elements;k++)
 554           position+=_point(v,pointindex[Bptr])[k]*n[k];
 555         if(position>0.)break; /* not in B */
 556         Bptr--;
 557       }
 558       if(Aptr<Bptr){
 559         long temp=pointindex[Aptr];
 560         pointindex[Aptr]=pointindex[Bptr];
 561         pointindex[Bptr]=temp;
 562       }
 563       pointsA=Aptr;
 564     }
 565   }
 566
 567   fprintf(stderr,"split: total=%ld depth=%ld set A=%ld:%ld:%ld=B\n",
 568           entries,depth,entriesA-entriesC,entriesC,entriesB-entriesC);
 569   {
 570     long thisaux=b->aux++;
 571     if(b->aux>=b->alloc){
 572       b->alloc*=2;
 573       b->ptr0=realloc(b->ptr0,sizeof(long)*b->alloc);
 574       b->ptr1=realloc(b->ptr1,sizeof(long)*b->alloc);
 575       b->n=realloc(b->n,sizeof(double)*b->elements*b->alloc);
 576       b->c=realloc(b->c,sizeof(double)*b->alloc);
 577     }
 578
 579     memcpy(b->n+b->elements*thisaux,n,sizeof(double)*v->elements);
 580     b->c[thisaux]=c;
 581
 582     if(entriesA==1){
 583       ret=1;
 584       b->ptr0[thisaux]=entryA[0];
 585     }else{
 586       b->ptr0[thisaux]= -b->aux;
 587       ret=lp_split(v,b,entryA,entriesA,pointindex,pointsA,depth+1,slack);
 588     }
 589     if(entriesB==1){
 590       ret++;
 591       b->ptr1[thisaux]=entryB[0];
 592     }else{
 593       b->ptr1[thisaux]= -b->aux;
 594       ret=lp_split(v,b,entryB,entriesB,pointindex+pointsA,points-pointsA,
 595                    depth+1,slack);
 596     }
 597   }
 598   free(entryA);
 599   free(entryB);
 600   return(ret);
 601 }
 602
 603 int vqenc_entry(vqbook *b,double *val){
 604   int ptr=0,k;
 605   while(1){
 606     double c= -b->c[ptr];
 607     double *nptr=b->n+b->elements*ptr;
 608     for(k=0;k<b->elements;k++)
 609       c+=nptr[k]*val[k];
 610     if(c>0.) /* in A */
 611       ptr= -b->ptr0[ptr];
 612     else     /* in B */
 613       ptr= -b->ptr1[ptr];
 614     if(ptr<=0)break;
 615   }
 616   return(-ptr);
 617 }
 618
 619 void vqgen_book(vqgen *v,vqbook *b){
 620   long i;
 621   long *entryindex=malloc(sizeof(double)*v->entries);
 622   long *pointindex=malloc(sizeof(double)*v->points);
 623
 624   memset(b,0,sizeof(vqbook));
 625   for(i=0;i<v->entries;i++)entryindex[i]=i;
 626   for(i=0;i<v->points;i++)pointindex[i]=i;
 627   b->elements=v->elements;
 628   b->entries=v->entries;
 629   b->alloc=4096;
 630   b->ptr0=malloc(sizeof(long)*b->alloc);
 631   b->ptr1=malloc(sizeof(long)*b->alloc);
 632   b->n=malloc(sizeof(double)*b->elements*b->alloc);
 633   b->c=malloc(sizeof(double)*b->alloc);
 634
 635   b->valuelist=malloc(sizeof(double)*b->elements*b->entries);
 636   b->codelist=malloc(sizeof(long)*b->entries);
 637   b->lengthlist=malloc(b->entries*sizeof(long));
 638
 639   /* first, generate the encoding decision heirarchy */
 640   lp_split(v,b,entryindex,v->entries, pointindex,v->points,0,0);
 641
 642   /* run all training points through the decision tree to get a final
 643      probability count */
 644   {
 645     long *probability=calloc(b->entries*2,sizeof(long));
 646     for(i=0;i<v->points;i++){
 647       int ret=vqenc_entry(b,v->pointlist+i*v->elements);
 648       probability[ret]++;
 649     }
 650
 651     for(i=0;i<b->entries;i++){
 652       fprintf(stderr,"point %ld: %ld\n",i,probability[i]);
 653     }
 654
 655   }
 656
 657   /* generate the codewords (short to long) */
 658
 659   free(entryindex);
 660   free(pointindex);
 661 }
 662
 663 static double testset24[48]={
 664
 665 };
 666
 667 static double testset256[1024]={
 668 0.334427,0.567149,0.749324,0.838460,
 669 0.212558,0.435792,0.661945,0.781195,
 670 0.316791,0.551190,0.700236,0.813437,
 671 0.240661,0.398754,0.570649,0.666909,
 672 0.233216,0.549751,0.715668,0.844393,
 673 0.275589,0.460791,0.687872,0.786855,
 674 0.145195,0.478688,0.678451,0.788536,
 675 0.332491,0.577111,0.770913,0.875295,
 676 0.162683,0.250422,0.399770,0.688120,
 677 0.301578,0.424074,0.595725,0.829533,
 678 0.118704,0.276233,0.492329,0.605532,
 679 0.240679,0.468862,0.670704,0.809156,
 680 0.109799,0.227123,0.533265,0.686982,
 681 0.119657,0.423823,0.597676,0.710405,
 682 0.207444,0.437121,0.671331,0.809054,
 683 0.222787,0.324523,0.466905,0.652758,
 684 0.117279,0.268017,0.410036,0.593785,
 685 0.316531,0.493867,0.642588,0.815791,
 686 0.218069,0.315630,0.554662,0.749348,
 687 0.152218,0.412864,0.664135,0.814275,
 688 0.144242,0.221503,0.457281,0.631403,
 689 0.290454,0.458806,0.603994,0.751802,
 690 0.236719,0.341259,0.586677,0.707653,
 691 0.257617,0.406941,0.602074,0.696609,
 692 0.141833,0.254061,0.523788,0.701270,
 693 0.136685,0.333920,0.515585,0.642137,
 694 0.191560,0.485163,0.665121,0.755073,
 695 0.086062,0.375207,0.623168,0.763659,
 696 0.184017,0.404999,0.609708,0.746118,
 697 0.208200,0.414804,0.618734,0.732746,
 698 0.184643,0.390701,0.666137,0.791445,
 699 0.216594,0.519656,0.713291,0.810097,
 700 0.118095,0.361088,0.577184,0.667553,
 701 0.243858,0.394640,0.550682,0.774866,
 702 0.266410,0.430656,0.629812,0.729035,
 703 0.140187,0.323258,0.534031,0.762617,
 704 0.097902,0.230404,0.382920,0.707036,
 705 0.197244,0.424752,0.578097,0.741857,
 706 0.216830,0.460883,0.622823,0.747350,
 707 0.096586,0.252731,0.471721,0.766199,
 708 0.385920,0.650477,0.846797,0.972222,
 709 0.249895,0.442701,0.706660,0.824760,
 710 0.207781,0.355197,0.538083,0.686972,
 711 0.166064,0.261252,0.493243,0.642141,
 712 0.146441,0.385197,0.595104,0.739789,
 713 0.207555,0.311523,0.533963,0.666395,
 714 0.172735,0.455702,0.624525,0.713503,
 715 0.132351,0.420175,0.528877,0.716175,
 716 0.113027,0.270164,0.462866,0.653668,
 717 0.306970,0.575262,0.785717,0.942889,
 718 0.122649,0.451363,0.652393,0.749624,
 719 0.337459,0.561671,0.713605,0.878091,
 720 0.154890,0.324078,0.469466,0.677888,
 721 0.242478,0.420602,0.651274,0.762209,
 722 0.162266,0.401114,0.639885,0.769255,
 723 0.275493,0.500799,0.753756,0.884563,
 724 0.219840,0.410217,0.634590,0.784979,
 725 0.197822,0.418074,0.514919,0.784821,
 726 0.176855,0.327821,0.624861,0.768299,
 727 0.187133,0.373303,0.614062,0.764326,
 728 0.148866,0.339151,0.571398,0.693017,
 729 0.126516,0.364169,0.643363,0.784478,
 730 0.176790,0.354609,0.553963,0.728466,
 731 0.169151,0.294951,0.559738,0.718315,
 732 0.189233,0.369137,0.615112,0.736337,
 733 0.211442,0.459206,0.587994,0.715496,
 734 0.203544,0.332755,0.496152,0.707495,
 735 0.264950,0.377661,0.601633,0.738910,
 736 0.144156,0.332387,0.533986,0.858065,
 737 0.236705,0.493254,0.646422,0.780372,
 738 0.247066,0.480770,0.725958,0.847099,
 739 0.121273,0.415837,0.575449,0.785059,
 740 0.109050,0.326623,0.585313,0.725416,
 741 0.079853,0.281243,0.579608,0.721251,
 742 0.164131,0.452931,0.608512,0.772505,
 743 0.183958,0.382443,0.595812,0.708979,
 744 0.232813,0.460979,0.612592,0.785098,
 745 0.372632,0.625661,0.794886,0.915078,
 746 0.197199,0.422499,0.630676,0.764479,
 747 0.112259,0.228495,0.423574,0.529402,
 748 0.313775,0.495449,0.684949,0.906438,
 749 0.121923,0.272895,0.590304,0.762871,
 750 0.114215,0.298370,0.523344,0.690963,
 751 0.195814,0.446435,0.567172,0.888808,
 752 0.292057,0.491911,0.707296,0.836253,
 753 0.198244,0.410661,0.555814,0.699832,
 754 0.152291,0.332174,0.609291,0.732500,
 755 0.051002,0.303083,0.609333,0.769208,
 756 0.086086,0.233029,0.475899,0.590504,
 757 0.096449,0.366927,0.502638,0.761407,
 758 0.256632,0.423960,0.650693,0.805566,
 759 0.235739,0.542293,0.776812,0.939851,
 760 0.210050,0.310752,0.439640,0.760752,
 761 0.168588,0.352547,0.620034,0.811813,
 762 0.217251,0.438798,0.704758,0.854804,
 763 0.112327,0.304063,0.583490,0.813809,
 764 0.183413,0.263207,0.503397,0.753610,
 765 0.241265,0.367568,0.648774,0.795899,
 766 0.262796,0.558878,0.784660,0.908486,
 767 0.221425,0.479136,0.683794,0.830981,
 768 0.325039,0.573201,0.837740,0.994924,
 769 0.263655,0.636322,0.878818,1.052555,
 770 0.228884,0.510746,0.742781,0.903443,
 771 0.169799,0.354550,0.667261,0.880659,
 772 0.429209,0.689908,0.942083,1.075538,
 773 0.178843,0.393247,0.698620,0.846058,
 774 0.410735,0.684575,0.900212,1.026610,
 775 0.212082,0.407634,0.677794,0.824167,
 776 0.183638,0.456984,0.701680,0.828091,
 777 0.198668,0.290621,0.586525,0.796566,
 778 0.212584,0.394786,0.624992,0.756258,
 779 0.218545,0.327989,0.671635,0.851695,
 780 0.118076,0.487280,0.758482,0.920421,
 781 0.184608,0.426118,0.642211,0.797123,
 782 0.248725,0.534657,0.737666,0.869414,
 783 0.075312,0.218113,0.574791,0.767607,
 784 0.184571,0.487955,0.742725,0.865892,
 785 0.148818,0.284078,0.481156,0.816008,
 786 0.205061,0.336518,0.506794,0.610053,
 787 0.054274,0.238425,0.414471,0.577301,
 788 0.082819,0.289140,0.632090,0.789039,
 789 0.095066,0.249569,0.666867,0.857491,
 790 0.129159,0.215969,0.524985,0.808056,
 791 0.120949,0.473613,0.723959,0.844274,
 792 0.225445,0.469240,0.781258,0.937501,
 793 0.190737,0.487420,0.664112,0.798043,
 794 0.249655,0.533041,0.808386,0.984218,
 795 0.236544,0.473767,0.723955,0.872327,
 796 0.093351,0.432384,0.672316,0.803956,
 797 0.171295,0.491945,0.662536,0.844976,
 798 0.165080,0.389674,0.533508,0.664579,
 799 0.224994,0.356513,0.535555,0.640525,
 800 0.158366,0.244279,0.437682,0.579242,
 801 0.267249,0.423934,0.542629,0.635151,
 802 0.105840,0.183727,0.395015,0.546705,
 803 0.128701,0.214397,0.323019,0.569708,
 804 0.205691,0.367797,0.491331,0.579727,
 805 0.136376,0.295543,0.461283,0.567165,
 806 0.090912,0.326863,0.449785,0.549052,
 807 0.280052,0.463123,0.590556,0.676314,
 808 0.053220,0.327610,0.496921,0.611028,
 809 0.074906,0.374797,0.563476,0.691693,
 810 0.064942,0.286131,0.414479,0.492806,
 811 0.303210,0.504309,0.652008,0.744313,
 812 0.322904,0.532050,0.682097,0.777012,
 813 0.247437,0.397965,0.504882,0.680784,
 814 0.076782,0.185773,0.355825,0.478750,
 815 0.170565,0.428925,0.587157,0.683893,
 816 0.186692,0.282933,0.468043,0.583405,
 817 0.286847,0.489852,0.626934,0.714018,
 818 0.196398,0.333200,0.447551,0.551070,
 819 0.061757,0.242351,0.503337,0.676114,
 820 0.119199,0.359308,0.514024,0.595332,
 821 0.282369,0.428107,0.549797,0.721344,
 822 0.120772,0.328758,0.465806,0.797635,
 823 0.251676,0.385411,0.692374,0.930173,
 824 0.265590,0.605988,0.828534,0.964088,
 825 0.118910,0.347784,0.690580,0.839342,
 826 0.148664,0.429126,0.564376,0.640949,
 827 0.273239,0.395545,0.492952,0.803872,
 828 0.085855,0.175468,0.449037,0.649337,
 829 0.101341,0.351796,0.455442,0.645061,
 830 0.181516,0.301506,0.381633,0.628656,
 831 0.165340,0.385706,0.549717,0.787060,
 832 0.102530,0.247439,0.362454,0.570843,
 833 0.115626,0.291845,0.411639,0.730918,
 834 0.130168,0.402586,0.531264,0.842151,
 835 0.312057,0.497535,0.608141,0.933018,
 836 0.061317,0.329143,0.457434,0.717368,
 837 0.260419,0.364385,0.522931,0.733634,
 838 0.287346,0.448695,0.703337,0.864549,
 839 0.158190,0.329088,0.444748,0.609073,
 840 0.057217,0.172765,0.500069,0.763990,
 841 0.185582,0.346521,0.426838,0.724792,
 842 0.152203,0.253061,0.580445,0.835384,
 843 0.223681,0.376501,0.463553,0.758668,
 844 0.153621,0.385013,0.482490,0.723052,
 845 0.165862,0.300622,0.515314,0.628741,
 846 0.087543,0.300967,0.529657,0.641654,
 847 0.257693,0.367954,0.453859,0.654989,
 848 0.077763,0.391224,0.586210,0.846828,
 849 0.140653,0.296311,0.431419,0.524173,
 850 0.276106,0.566792,0.723587,0.928366,
 851 0.345768,0.530982,0.711197,1.011075,
 852 0.328643,0.578965,0.775015,1.036448,
 853 0.295489,0.451392,0.541007,0.831967,
 854 0.249971,0.517649,0.662078,0.872210,
 855 0.171802,0.382312,0.492482,0.629808,
 856 0.192277,0.478000,0.611564,0.846436,
 857 0.217769,0.335272,0.508986,0.813008,
 858 0.265875,0.429911,0.625892,0.958544,
 859 0.183329,0.299421,0.565539,0.885282,
 860 0.311307,0.432098,0.756371,0.937483,
 861 0.153908,0.264165,0.391052,0.502304,
 862 0.267528,0.476930,0.862271,1.099421,
 863 0.221527,0.513611,0.684819,0.993058,
 864 0.220330,0.352263,0.575808,0.823017,
 865 0.084526,0.372939,0.705268,0.879314,
 866 0.179752,0.287765,0.507176,0.683204,
 867 0.234005,0.381653,0.537092,0.887657,
 868 0.068805,0.135538,0.363455,0.732055,
 869 0.266267,0.440971,0.719076,1.025917,
 870 0.169728,0.426510,0.670186,0.982939,
 871 0.248443,0.463944,0.736167,0.898651,
 872 0.074789,0.155652,0.601987,0.864470,
 873 0.394203,0.665959,0.950879,1.183016,
 874 0.056927,0.247163,0.620166,0.886961,
 875 0.081110,0.260015,0.725344,0.969987,
 876 0.187566,0.605658,0.907862,1.113592,
 877 0.188174,0.417467,0.746757,0.912718,
 878 0.214944,0.444902,0.813762,1.006999,
 879 0.177305,0.322288,0.623538,0.972940,
 880 0.108933,0.560083,0.892489,1.105198,
 881 0.326010,0.486890,0.808596,1.012135,
 882 0.507702,0.815903,1.085892,1.267201,
 883 0.121027,0.368451,0.715029,1.000456,
 884 0.424229,0.713942,0.971720,1.123569,
 885 0.158148,0.411905,0.812352,1.057000,
 886 0.330896,0.639677,0.814721,1.083593,
 887 0.116143,0.314050,0.643242,0.928001,
 888 0.228088,0.364437,0.747642,1.045211,
 889 0.457076,0.783579,1.021362,1.174221,
 890 0.069613,0.430727,0.845847,1.061835,
 891 0.343281,0.553869,0.865545,1.077270,
 892 0.152973,0.304560,0.642854,0.830441,
 893 0.157602,0.568564,0.815800,1.036183,
 894 0.104003,0.505075,0.796932,0.999392,
 895 0.377220,0.622243,0.901761,1.072988,
 896 0.060039,0.378222,0.710910,0.939836,
 897 0.220206,0.539490,0.849893,1.056103,
 898 0.116884,0.454544,0.627779,0.924877,
 899 0.293087,0.633251,0.927781,1.114120,
 900 0.234659,0.458666,0.693916,0.925308,
 901 0.271802,0.385786,0.661508,0.862002,
 902 0.185212,0.282842,0.431237,0.663853,
 903 0.152441,0.280173,0.711942,0.961418,
 904 0.264191,0.647776,0.976337,1.171963,
 905 0.227340,0.567363,0.753334,1.049249,
 906 0.173306,0.508941,0.761828,0.970236,
 907 0.271942,0.496903,0.759621,0.963074,
 908 0.187050,0.355811,0.752185,0.952287,
 909 0.172003,0.363447,0.486011,0.829142,
 910 0.364755,0.598200,0.868222,1.025048,
 911 0.226504,0.432744,0.649002,0.852690,
 912 0.200857,0.436929,0.651778,0.905448,
 913 0.137860,0.385245,0.609566,0.867573,
 914 0.217881,0.368030,0.601419,0.954770,
 915 0.199237,0.398640,0.593487,0.830535,
 916 0.324623,0.518989,0.782275,0.926822,
 917 0.245185,0.419081,0.717432,0.883462,
 918 0.277741,0.494080,0.765001,0.915227,
 919 0.200418,0.463680,0.705028,0.875441,
 920 0.171638,0.526057,0.697163,0.915452,
 921 0.142351,0.416539,0.707280,0.876499,
 922 0.178309,0.553184,0.777786,0.903243,
 923 0.335361,0.540299,0.812589,0.965039,
 924
 925
 926 };
 927
 928 int main(int argc,char *argv[]){
 929   FILE *in=fopen(argv[1],"r");
 930   vqgen v;
 931   vqbook b;
 932   char buffer[160];
 933   int i,j;
 934
 935   vqgen_init(&v,4,256,_dist_and_pos,0.);
 936
 937   while(fgets(buffer,160,in)){
 938     double a[8];
 939     if(sscanf(buffer,"%lf %lf %lf %lf",
 940               a,a+1,a+2,a+3)==4)
 941       vqgen_addpoint(&v,a);
 942   }
 943   fclose(in);
 944
 945   /*for(i=0;i<200;i++){
 946     vqgen_iterate(&v);
 947   }
 948
 949   for(i=0;i<v.entries;i++){
 950     printf("\n");
 951     for(j=0;j<v.elements;j++)
 952       printf("%f,",_now(&v,i)[j]);
 953   }
 954   printf("\n");
 955
 956   exit(0);*/
 957
 958
 959   memcpy(v.entrylist,testset256,sizeof(testset256));
 960
 961   vqgen_book(&v,&b);
 962
 963   return(0);
 964 }
 965
 966
 967
 968
 969
 970
 971
 972
 973
 974
 975
 976
 977
 978