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