Home | History | Annotate | Download | only in vq 
      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 Xiph.Org Foundation http://www.xiph.org/                  *
     10  *                                                                  *
     11  ********************************************************************
     12 
     13  function: train a VQ codebook
     14  last mod: $Id: vqgen.c 16037 2009-05-26 21:10:58Z xiphmont $
     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 
    568