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