]> git.jsancho.org Git - lugaru.git/blob - libvorbis-1.0.1/vq/vqgen.c
Simplified web browser launching on Linux.
[lugaru.git] / libvorbis-1.0.1 / vq / vqgen.c
1 /********************************************************************
2  *                                                                  *
3  * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE.   *
4  * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS     *
5  * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE *
6  * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING.       *
7  *                                                                  *
8  * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2001             *
9  * by the XIPHOPHORUS Company http://www.xiph.org/                  *
10  *                                                                  *
11  ********************************************************************
12
13  function: train a VQ codebook 
14  last mod: $Id: vqgen.c,v 1.41 2002/10/11 07:44:28 xiphmont Exp $
15
16  ********************************************************************/
17
18 /* This code is *not* part of libvorbis.  It is used to generate
19    trained codebooks offline and then spit the results into a
20    pregenerated codebook that is compiled into libvorbis.  It is an
21    expensive (but good) algorithm.  Run it on big iron. */
22
23 /* There are so many optimizations to explore in *both* stages that
24    considering the undertaking is almost withering.  For now, we brute
25    force it all */
26
27 #include <stdlib.h>
28 #include <stdio.h>
29 #include <math.h>
30 #include <string.h>
31
32 #include "vqgen.h"
33 #include "bookutil.h"
34
35 /* Codebook generation happens in two steps: 
36
37    1) Train the codebook with data collected from the encoder: We use
38    one of a few error metrics (which represent the distance between a
39    given data point and a candidate point in the training set) to
40    divide the training set up into cells representing roughly equal
41    probability of occurring. 
42
43    2) Generate the codebook and auxiliary data from the trained data set
44 */
45
46 /* Codebook training ****************************************************
47  *
48  * The basic idea here is that a VQ codebook is like an m-dimensional
49  * foam with n bubbles.  The bubbles compete for space/volume and are
50  * 'pressurized' [biased] according to some metric.  The basic alg
51  * iterates through allowing the bubbles to compete for space until
52  * they converge (if the damping is dome properly) on a steady-state
53  * solution. Individual input points, collected from libvorbis, are
54  * used to train the algorithm monte-carlo style.  */
55
56 /* internal helpers *****************************************************/
57 #define vN(data,i) (data+v->elements*i)
58
59 /* default metric; squared 'distance' from desired value. */
60 float _dist(vqgen *v,float *a, float *b){
61   int i;
62   int el=v->elements;
63   float acc=0.f;
64   for(i=0;i<el;i++){
65     float val=(a[i]-b[i]);
66     acc+=val*val;
67   }
68   return sqrt(acc);
69 }
70
71 float *_weight_null(vqgen *v,float *a){
72   return a;
73 }
74
75 /* *must* be beefed up. */
76 void _vqgen_seed(vqgen *v){
77   long i;
78   for(i=0;i<v->entries;i++)
79     memcpy(_now(v,i),_point(v,i),sizeof(float)*v->elements);
80   v->seeded=1;
81 }
82
83 int directdsort(const void *a, const void *b){
84   float av=*((float *)a);
85   float bv=*((float *)b);
86   return (av<bv)-(av>bv);
87 }
88
89 void vqgen_cellmetric(vqgen *v){
90   int j,k;
91   float min=-1.f,max=-1.f,mean=0.f,acc=0.f;
92   long dup=0,unused=0;
93  #ifdef NOISY
94   int i;
95    char buff[80];
96    float spacings[v->entries];
97    int count=0;
98    FILE *cells;
99    sprintf(buff,"cellspace%d.m",v->it);
100    cells=fopen(buff,"w");
101 #endif
102
103   /* minimum, maximum, cell spacing */
104   for(j=0;j<v->entries;j++){
105     float localmin=-1.;
106
107     for(k=0;k<v->entries;k++){
108       if(j!=k){
109         float this=_dist(v,_now(v,j),_now(v,k));
110         if(this>0){
111           if(v->assigned[k] && (localmin==-1 || this<localmin))
112             localmin=this;
113         }else{  
114           if(k<j){
115             dup++;
116             break;
117           }
118         }
119       }
120     }
121     if(k<v->entries)continue;
122
123     if(v->assigned[j]==0){
124       unused++;
125       continue;
126     }
127     
128     localmin=v->max[j]+localmin/2; /* this gives us rough diameter */
129     if(min==-1 || localmin<min)min=localmin;
130     if(max==-1 || localmin>max)max=localmin;
131     mean+=localmin;
132     acc++;
133 #ifdef NOISY
134     spacings[count++]=localmin;
135 #endif
136   }
137
138   fprintf(stderr,"cell diameter: %.03g::%.03g::%.03g (%ld unused/%ld dup)\n",
139           min,mean/acc,max,unused,dup);
140
141 #ifdef NOISY
142   qsort(spacings,count,sizeof(float),directdsort);
143   for(i=0;i<count;i++)
144     fprintf(cells,"%g\n",spacings[i]);
145   fclose(cells);
146 #endif      
147
148 }
149
150 /* External calls *******************************************************/
151
152 /* We have two forms of quantization; in the first, each vector
153    element in the codebook entry is orthogonal.  Residues would use this
154    quantization for example.
155
156    In the second, we have a sequence of monotonically increasing
157    values that we wish to quantize as deltas (to save space).  We
158    still need to quantize so that absolute values are accurate. For
159    example, LSP quantizes all absolute values, but the book encodes
160    distance between values because each successive value is larger
161    than the preceeding value.  Thus the desired quantibits apply to
162    the encoded (delta) values, not abs positions. This requires minor
163    additional encode-side trickery. */
164
165 void vqgen_quantize(vqgen *v,quant_meta *q){
166
167   float maxdel;
168   float mindel;
169
170   float delta;
171   float maxquant=((1<<q->quant)-1);
172
173   int j,k;
174
175   mindel=maxdel=_now(v,0)[0];
176   
177   for(j=0;j<v->entries;j++){
178     float last=0.f;
179     for(k=0;k<v->elements;k++){
180       if(mindel>_now(v,j)[k]-last)mindel=_now(v,j)[k]-last;
181       if(maxdel<_now(v,j)[k]-last)maxdel=_now(v,j)[k]-last;
182       if(q->sequencep)last=_now(v,j)[k];
183     }
184   }
185
186
187   /* first find the basic delta amount from the maximum span to be
188      encoded.  Loosen the delta slightly to allow for additional error
189      during sequence quantization */
190
191   delta=(maxdel-mindel)/((1<<q->quant)-1.5f);
192
193   q->min=_float32_pack(mindel);
194   q->delta=_float32_pack(delta);
195
196   mindel=_float32_unpack(q->min);
197   delta=_float32_unpack(q->delta);
198
199   for(j=0;j<v->entries;j++){
200     float last=0;
201     for(k=0;k<v->elements;k++){
202       float val=_now(v,j)[k];
203       float now=rint((val-last-mindel)/delta);
204       
205       _now(v,j)[k]=now;
206       if(now<0){
207         /* be paranoid; this should be impossible */
208         fprintf(stderr,"fault; quantized value<0\n");
209         exit(1);
210       }
211
212       if(now>maxquant){
213         /* be paranoid; this should be impossible */
214         fprintf(stderr,"fault; quantized value>max\n");
215         exit(1);
216       }
217       if(q->sequencep)last=(now*delta)+mindel+last;
218     }
219   }
220 }
221
222 /* much easier :-).  Unlike in the codebook, we don't un-log log
223    scales; we just make sure they're properly offset. */
224 void vqgen_unquantize(vqgen *v,quant_meta *q){
225   long j,k;
226   float mindel=_float32_unpack(q->min);
227   float delta=_float32_unpack(q->delta);
228
229   for(j=0;j<v->entries;j++){
230     float last=0.f;
231     for(k=0;k<v->elements;k++){
232       float now=_now(v,j)[k];
233       now=fabs(now)*delta+last+mindel;
234       if(q->sequencep)last=now;
235       _now(v,j)[k]=now;
236     }
237   }
238 }
239
240 void vqgen_init(vqgen *v,int elements,int aux,int entries,float mindist,
241                 float  (*metric)(vqgen *,float *, float *),
242                 float *(*weight)(vqgen *,float *),int centroid){
243   memset(v,0,sizeof(vqgen));
244
245   v->centroid=centroid;
246   v->elements=elements;
247   v->aux=aux;
248   v->mindist=mindist;
249   v->allocated=32768;
250   v->pointlist=_ogg_malloc(v->allocated*(v->elements+v->aux)*sizeof(float));
251
252   v->entries=entries;
253   v->entrylist=_ogg_malloc(v->entries*v->elements*sizeof(float));
254   v->assigned=_ogg_malloc(v->entries*sizeof(long));
255   v->bias=_ogg_calloc(v->entries,sizeof(float));
256   v->max=_ogg_calloc(v->entries,sizeof(float));
257   if(metric)
258     v->metric_func=metric;
259   else
260     v->metric_func=_dist;
261   if(weight)
262     v->weight_func=weight;
263   else
264     v->weight_func=_weight_null;
265
266   v->asciipoints=tmpfile();
267
268 }
269
270 void vqgen_addpoint(vqgen *v, float *p,float *a){
271   int k;
272   for(k=0;k<v->elements;k++)
273     fprintf(v->asciipoints,"%.12g\n",p[k]);
274   for(k=0;k<v->aux;k++)
275     fprintf(v->asciipoints,"%.12g\n",a[k]);
276
277   if(v->points>=v->allocated){
278     v->allocated*=2;
279     v->pointlist=_ogg_realloc(v->pointlist,v->allocated*(v->elements+v->aux)*
280                          sizeof(float));
281   }
282
283   memcpy(_point(v,v->points),p,sizeof(float)*v->elements);
284   if(v->aux)memcpy(_point(v,v->points)+v->elements,a,sizeof(float)*v->aux);
285  
286   /* quantize to the density mesh if it's selected */
287   if(v->mindist>0.f){
288     /* quantize to the mesh */
289     for(k=0;k<v->elements+v->aux;k++)
290       _point(v,v->points)[k]=
291         rint(_point(v,v->points)[k]/v->mindist)*v->mindist;
292   }
293   v->points++;
294   if(!(v->points&0xff))spinnit("loading... ",v->points);
295 }
296
297 /* yes, not threadsafe.  These utils aren't */
298 static int sortit=0;
299 static int sortsize=0;
300 static int meshcomp(const void *a,const void *b){
301   if(((sortit++)&0xfff)==0)spinnit("sorting mesh...",sortit);
302   return(memcmp(a,b,sortsize));
303 }
304
305 void vqgen_sortmesh(vqgen *v){
306   sortit=0;
307   if(v->mindist>0.f){
308     long i,march=1;
309
310     /* sort to make uniqueness detection trivial */
311     sortsize=(v->elements+v->aux)*sizeof(float);
312     qsort(v->pointlist,v->points,sortsize,meshcomp);
313
314     /* now march through and eliminate dupes */
315     for(i=1;i<v->points;i++){
316       if(memcmp(_point(v,i),_point(v,i-1),sortsize)){
317         /* a new, unique entry.  march it down */
318         if(i>march)memcpy(_point(v,march),_point(v,i),sortsize);
319         march++;
320       }
321       spinnit("eliminating density... ",v->points-i);
322     }
323
324     /* we're done */
325     fprintf(stderr,"\r%ld training points remining out of %ld"
326             " after density mesh (%ld%%)\n",march,v->points,march*100/v->points);
327     v->points=march;
328
329   }
330   v->sorted=1;
331 }
332
333 float vqgen_iterate(vqgen *v,int biasp){
334   long   i,j,k;
335
336   float fdesired;
337   long  desired;
338   long  desired2;
339
340   float asserror=0.f;
341   float meterror=0.f;
342   float *new;
343   float *new2;
344   long   *nearcount;
345   float *nearbias;
346  #ifdef NOISY
347    char buff[80];
348    FILE *assig;
349    FILE *bias;
350    FILE *cells;
351    sprintf(buff,"cells%d.m",v->it);
352    cells=fopen(buff,"w");
353    sprintf(buff,"assig%d.m",v->it);
354    assig=fopen(buff,"w");
355    sprintf(buff,"bias%d.m",v->it);
356    bias=fopen(buff,"w");
357  #endif
358  
359
360   if(v->entries<2){
361     fprintf(stderr,"generation requires at least two entries\n");
362     exit(1);
363   }
364
365   if(!v->sorted)vqgen_sortmesh(v);
366   if(!v->seeded)_vqgen_seed(v);
367
368   fdesired=(float)v->points/v->entries;
369   desired=fdesired;
370   desired2=desired*2;
371   new=_ogg_malloc(sizeof(float)*v->entries*v->elements);
372   new2=_ogg_malloc(sizeof(float)*v->entries*v->elements);
373   nearcount=_ogg_malloc(v->entries*sizeof(long));
374   nearbias=_ogg_malloc(v->entries*desired2*sizeof(float));
375
376   /* fill in nearest points for entry biasing */
377   /*memset(v->bias,0,sizeof(float)*v->entries);*/
378   memset(nearcount,0,sizeof(long)*v->entries);
379   memset(v->assigned,0,sizeof(long)*v->entries);
380   if(biasp){
381     for(i=0;i<v->points;i++){
382       float *ppt=v->weight_func(v,_point(v,i));
383       float firstmetric=v->metric_func(v,_now(v,0),ppt)+v->bias[0];
384       float secondmetric=v->metric_func(v,_now(v,1),ppt)+v->bias[1];
385       long   firstentry=0;
386       long   secondentry=1;
387       
388       if(!(i&0xff))spinnit("biasing... ",v->points+v->points+v->entries-i);
389       
390       if(firstmetric>secondmetric){
391         float temp=firstmetric;
392         firstmetric=secondmetric;
393         secondmetric=temp;
394         firstentry=1;
395         secondentry=0;
396       }
397       
398       for(j=2;j<v->entries;j++){
399         float thismetric=v->metric_func(v,_now(v,j),ppt)+v->bias[j];
400         if(thismetric<secondmetric){
401           if(thismetric<firstmetric){
402             secondmetric=firstmetric;
403             secondentry=firstentry;
404             firstmetric=thismetric;
405             firstentry=j;
406           }else{
407             secondmetric=thismetric;
408             secondentry=j;
409           }
410         }
411       }
412       
413       j=firstentry;
414       for(j=0;j<v->entries;j++){
415         
416         float thismetric,localmetric;
417         float *nearbiasptr=nearbias+desired2*j;
418         long k=nearcount[j];
419         
420         localmetric=v->metric_func(v,_now(v,j),ppt);
421         /* 'thismetric' is to be the bias value necessary in the current
422            arrangement for entry j to capture point i */
423         if(firstentry==j){
424           /* use the secondary entry as the threshhold */
425           thismetric=secondmetric-localmetric;
426         }else{
427           /* use the primary entry as the threshhold */
428           thismetric=firstmetric-localmetric;
429         }
430         
431         /* support the idea of 'minimum distance'... if we want the
432            cells in a codebook to be roughly some minimum size (as with
433            the low resolution residue books) */
434         
435         /* a cute two-stage delayed sorting hack */
436         if(k<desired){
437           nearbiasptr[k]=thismetric;
438           k++;
439           if(k==desired){
440             spinnit("biasing... ",v->points+v->points+v->entries-i);
441             qsort(nearbiasptr,desired,sizeof(float),directdsort);
442           }
443           
444         }else if(thismetric>nearbiasptr[desired-1]){
445           nearbiasptr[k]=thismetric;
446           k++;
447           if(k==desired2){
448             spinnit("biasing... ",v->points+v->points+v->entries-i);
449             qsort(nearbiasptr,desired2,sizeof(float),directdsort);
450             k=desired;
451           }
452         }
453         nearcount[j]=k;
454       }
455     }
456     
457     /* inflate/deflate */
458     
459     for(i=0;i<v->entries;i++){
460       float *nearbiasptr=nearbias+desired2*i;
461       
462       spinnit("biasing... ",v->points+v->entries-i);
463       
464       /* due to the delayed sorting, we likely need to finish it off....*/
465       if(nearcount[i]>desired)
466         qsort(nearbiasptr,nearcount[i],sizeof(float),directdsort);
467
468       v->bias[i]=nearbiasptr[desired-1];
469
470     }
471   }else{ 
472     memset(v->bias,0,v->entries*sizeof(float));
473   }
474
475   /* Now assign with new bias and find new midpoints */
476   for(i=0;i<v->points;i++){
477     float *ppt=v->weight_func(v,_point(v,i));
478     float firstmetric=v->metric_func(v,_now(v,0),ppt)+v->bias[0];
479     long   firstentry=0;
480
481     if(!(i&0xff))spinnit("centering... ",v->points-i);
482
483     for(j=0;j<v->entries;j++){
484       float thismetric=v->metric_func(v,_now(v,j),ppt)+v->bias[j];
485       if(thismetric<firstmetric){
486         firstmetric=thismetric;
487         firstentry=j;
488       }
489     }
490
491     j=firstentry;
492       
493 #ifdef NOISY
494     fprintf(cells,"%g %g\n%g %g\n\n",
495           _now(v,j)[0],_now(v,j)[1],
496           ppt[0],ppt[1]);
497 #endif
498
499     firstmetric-=v->bias[j];
500     meterror+=firstmetric;
501
502     if(v->centroid==0){
503       /* set up midpoints for next iter */
504       if(v->assigned[j]++){
505         for(k=0;k<v->elements;k++)
506           vN(new,j)[k]+=ppt[k];
507         if(firstmetric>v->max[j])v->max[j]=firstmetric;
508       }else{
509         for(k=0;k<v->elements;k++)
510           vN(new,j)[k]=ppt[k];
511         v->max[j]=firstmetric;
512       }
513     }else{
514       /* centroid */
515       if(v->assigned[j]++){
516         for(k=0;k<v->elements;k++){
517           if(vN(new,j)[k]>ppt[k])vN(new,j)[k]=ppt[k];
518           if(vN(new2,j)[k]<ppt[k])vN(new2,j)[k]=ppt[k];
519         }
520         if(firstmetric>v->max[firstentry])v->max[j]=firstmetric;
521       }else{
522         for(k=0;k<v->elements;k++){
523           vN(new,j)[k]=ppt[k];
524           vN(new2,j)[k]=ppt[k];
525         }
526         v->max[firstentry]=firstmetric;
527       }
528     }
529   }
530
531   /* assign midpoints */
532
533   for(j=0;j<v->entries;j++){
534 #ifdef NOISY
535     fprintf(assig,"%ld\n",v->assigned[j]);
536     fprintf(bias,"%g\n",v->bias[j]);
537 #endif
538     asserror+=fabs(v->assigned[j]-fdesired);
539     if(v->assigned[j]){
540       if(v->centroid==0){
541         for(k=0;k<v->elements;k++)
542           _now(v,j)[k]=vN(new,j)[k]/v->assigned[j];
543       }else{
544         for(k=0;k<v->elements;k++)
545           _now(v,j)[k]=(vN(new,j)[k]+vN(new2,j)[k])/2.f;
546       }
547     }
548   }
549
550   asserror/=(v->entries*fdesired);
551
552   fprintf(stderr,"Pass #%d... ",v->it);
553   fprintf(stderr,": dist %g(%g) metric error=%g \n",
554           asserror,fdesired,meterror/v->points);
555   v->it++;
556   
557   free(new);
558   free(nearcount);
559   free(nearbias);
560 #ifdef NOISY
561   fclose(assig);
562   fclose(bias);
563   fclose(cells);
564 #endif
565   return(asserror);
566 }
567