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 04 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]);
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);
95 void _limited_sort(double *vals,long *index,long total,int desired){
98 if(bisect>1)_limited_sort(vals,index,bisect,desired);
99 if(total-bisect>1)_limited_sort(vals,index+bisect,total-bisect,
101 if(desired>total)desired=total;
107 long *temp=alloca(desired*sizeof(long));
109 for(i=0;i<desired;i++){
110 if(vals[index[ptra]]>vals[index[ptrb]]){
111 temp[i]=index[ptra++];
113 temp[i]=index[ptrb++];
116 for(i++;i<desired;i++)temp[i]=index[ptrb++];
119 for(i++;i<desired;i++)temp[i]=index[ptra++];
123 memcpy(index,temp,desired*sizeof(long));
128 /* External calls *******************************************************/
130 void vqgen_init(vqgen *v,int elements,int entries,
131 double (*metric)(vqgen *,double *, double *),
133 memset(v,0,sizeof(vqgen));
135 v->elements=elements;
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));
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));
147 v->metric_func=metric;
149 v->metric_func=_dist_sq;
152 void vqgen_addpoint(vqgen *v, double *p){
153 if(v->points>=v->allocated){
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));
160 memcpy(_point(v,v->points),p,sizeof(double)*v->elements);
162 if(v->points==v->entries)_vqgen_seed(v);
165 void vqgen_recenter(vqgen *v){
167 double fdesired=(double)v->points/v->entries;
171 if(v->entries<2)exit(1);
173 /* first round: new midpoints */
175 double *newlo=malloc(sizeof(double)*v->entries*v->elements);
176 double *newhi=malloc(sizeof(double)*v->entries*v->elements);
178 memset(v->assigned,0,sizeof(long)*v->entries);
180 for(i=0;i<v->points;i++){
181 double firstmetric=v->metric_func(v,_now(v,0),_point(v,i))+v->bias[0];
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;
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];
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];
207 for(j=0;j<v->entries;j++){
209 for(k=0;k<v->elements;k++)
210 _now(v,j)[k]=(vN(newlo,j)[k]+vN(newhi,j)[k])/2.;
217 for(i=0;i<v->entries;i++){
218 asserror+=fabs(fdesired-v->assigned[i]);
221 fprintf(stderr,"recenter: dist error=%g(%g) metric error=%g\n",
222 asserror/v->entries,fdesired,meterror/v->points);
224 memset(v->bias,0,sizeof(double)*v->entries);
227 void vqgen_rebias(vqgen *v){
229 double fdesired=(float)v->points/v->entries;
230 long desired=(v->points+v->entries-1)/v->entries;
234 long *nearcount=calloc(v->entries,sizeof(long));
235 double *nearbias=malloc(v->entries*desired*sizeof(double));
237 if(v->entries<2)exit(1);
239 /* second round: fill in nearest points for entries */
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));
248 if(firstmetric>secondmetric){
249 double temp=firstmetric;
250 firstmetric=secondmetric;
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;
265 secondmetric=thismetric;
270 v->assigned[firstentry]++;
271 meterror+=firstmetric;
273 for(j=0;j<v->entries;j++){
276 double *nearbiasptr=nearbias+desired*j;
277 long k=nearcount[j]-1;
279 /* 'thismetric' is to be the bias value necessary in the current
280 arrangement for entry j to capture point i */
282 /* use the secondary entry as the threshhold */
283 thismetric=secondmetric-v->metric_func(v,_now(v,j),_point(v,i));
285 /* use the primary entry as the threshhold */
286 thismetric=firstmetric-v->metric_func(v,_now(v,j),_point(v,i));
289 if(k>=0 && thismetric>nearbiasptr[k]){
291 /* start at the end and search backward for where this entry
294 for(;k>0;k--) if(nearbiasptr[k-1]>=thismetric)break;
296 /* insert at k. Shift and inject. */
297 memmove(nearbiasptr+k+1,nearbiasptr+k,(desired-k-1)*sizeof(double));
298 nearbiasptr[k]=thismetric;
300 if(nearcount[j]<desired)nearcount[j]++;
303 if(nearcount[j]<desired){
304 /* we checked the thresh earlier. We know this is the
306 nearbiasptr[nearcount[j]++]=thismetric;
312 /* third round: inflate/deflate */
314 for(i=0;i<v->entries;i++){
315 v->bias[i]=nearbias[(i+1)*desired-1];
319 for(i=0;i<v->entries;i++){
320 asserror+=fabs(fdesired-v->assigned[i]);
323 fprintf(stderr,"rebias: dist error=%g(%g) metric error=%g\n",
324 asserror/v->entries,fdesired,meterror/v->points);
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 */
334 int lp_split(vqgen *v,long *index,long points,
335 double *additional_a, double *additional_b, long depth){
336 static int frameno=0;
350 printf("leaf: points %ld, depth %ld\n",points,depth);
355 printf("leaf: point %ld, depth %ld\n",index[0],depth);
360 /* The result must be an even split */
363 printf("even split: depth %ld\n",depth);
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){
372 best=_dist_sq(v,_now(v,index[i]),_now(v,index[j]));
377 /* We have our endpoints. initial divvy */
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));
387 minit_init(&l,v->elements,points+depth-1,0);
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];
398 for(i=0;i<points+depth;i++)
399 a[i]=aa+i*v->elements;
401 /* check each bubble to see if it intersects the divisor hyperplane */
402 for(j=0;j<points;j++){
413 for(i=0;i<points;i++){
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;
425 /* additional bounding hyperplanes from previous splits */
426 for(i=0;i<depth;i++){
428 for(k=0;k<v->elements;k++)
429 a[count][k]=additional_a[m++];
430 b[count++]=additional_b[i];
433 /* on what side of the dividing hyperplane is the current test
435 for(i=0;i<v->elements;i++)
436 d1+=_now(v,index[j])[i]*cA[i];
441 ret=minit_solve(&l,a,b,cA,1e-6,NULL,NULL,&d2);
444 ret=minit_solve(&l,a,b,cB,1e-6,NULL,NULL,&d2);
450 /* bounded solution */
453 if(d2>dA){cp++;B[bp++]=index[j];}
456 if(d2<dA){cp++;A[ap++]=index[j];}
460 /* unbounded solution or no solution; we're in both sets */
474 sprintf(buffer,"set%d.m",frameno);
476 for(i=0;i<points;i++)
477 fprintf(o,"%g %g\n\n",_now(v,index[i])[0],_now(v,index[i])[1]);
480 sprintf(buffer,"setA%d.m",frameno);
483 fprintf(o,"%g %g\n\n",_now(v,A[i])[0],_now(v,A[i])[1]);
486 sprintf(buffer,"setB%d.m",frameno);
489 fprintf(o,"%g %g\n\n",_now(v,B[i])[0],_now(v,B[i])[1]);
492 sprintf(buffer,"div%d.m",frameno);
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]);
501 sprintf(buffer,"bound%d.m",frameno);
504 fprintf(o,"%g %g\n%g %g\n\n",
505 _now(v,index[besti])[0],
507 additional_a[i*2]*_now(v,index[besti])[0])/
509 _now(v,index[bestj])[0],
511 additional_a[i*2]*_now(v,index[bestj])[0])/
512 additional_a[i*2+1]);
523 fprintf(stderr,"split: total=%ld depth=%ld set A=%ld:%ld:%ld=B\n",
524 points,depth,ap-cp,cp,bp-cp);
526 /* handle adding this divisor hyperplane to the additional
527 constraints. The 'inside' side is different for A and B */
530 /* marginally wasteful */
531 double *newadd_a=alloca(sizeof(double)*v->elements*(depth+1));
532 double *newadd_b=alloca(sizeof(double)*(depth+1));
534 memcpy(newadd_a,additional_a,sizeof(double)*v->elements*depth);
535 memcpy(newadd_b,additional_b,sizeof(double)*depth);
538 memcpy(newadd_a+v->elements*depth,cA,sizeof(double)*v->elements);
540 leaves=lp_split(v,A,ap,newadd_a,newadd_b,depth+1);
542 memcpy(newadd_a+v->elements*depth,cB,sizeof(double)*v->elements);
544 leaves+=lp_split(v,B,bp,newadd_a,newadd_b,depth+1);
554 void vqgen_book(vqgen *v){
561 static double testset24[48]={
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,
848 int main(int argc,char *argv[]){
849 FILE *in=fopen(argv[1],"r");
854 vqgen_init(&v,4,256,_dist_sq,0.);
856 while(fgets(buffer,160,in)){
858 if(sscanf(buffer,"%lf %lf %lf %lf",
860 vqgen_addpoint(&v,a);
870 fprintf(stderr,"%d ",i);
877 memcpy(v.entrylist,testset256,sizeof(testset256));
879 for(i=0;i<v.entries;i++){
881 for(j=0;j<v.elements;j++)
882 printf("%f,",_now(&v,i)[j]);
887 /* try recursively splitting the space using LP
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));