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-2000             *
   9  * by Monty <monty@xiph.org> and The XIPHOPHORUS Company            *
  10  * http://www.xiph.org/                                             *
  11  *                                                                  *
  12  ********************************************************************
  13
  14  function: train a VQ codebook
  15  last mod: $Id: vqgen.c,v 1.26 2000/01/05 10:14:57 xiphmont Exp $
  16
  17  ********************************************************************/
  18
  19 /* This code is *not* part of libvorbis.  It is used to generate
  20    trained codebooks offline and then spit the results into a
  21    pregenerated codebook that is compiled into libvorbis.  It is an
  22    expensive (but good) algorithm.  Run it on big iron. */
  23
  24 /* There are so many optimizations to explore in *both* stages that
  25    considering the undertaking is almost withering.  For now, we brute
  26    force it all */
  27
  28 #include <stdlib.h>
  29 #include <stdio.h>
  30 #include <math.h>
  31 #include <string.h>
  32
  33 #include "vqgen.h"
  34 #include "bookutil.h"
  35
  36 /* Codebook generation happens in two steps:
  37
  38    1) Train the codebook with data collected from the encoder: We use
  39    one of a few error metrics (which represent the distance between a
  40    given data point and a candidate point in the training set) to
  41    divide the training set up into cells representing roughly equal
  42    probability of occurring.
  43
  44    2) Generate the codebook and auxiliary data from the trained data set
  45 */
  46
  47 /* Codebook training ****************************************************
  48  *
  49  * The basic idea here is that a VQ codebook is like an m-dimensional
  50  * foam with n bubbles.  The bubbles compete for space/volume and are
  51  * 'pressurized' [biased] according to some metric.  The basic alg
  52  * iterates through allowing the bubbles to compete for space until
  53  * they converge (if the damping is dome properly) on a steady-state
  54  * solution. Individual input points, collected from libvorbis, are
  55  * used to train the algorithm monte-carlo style.  */
  56
  57 /* internal helpers *****************************************************/
  58 #define vN(data,i) (data+v->elements*i)
  59
  60 /* default metric; squared 'distance' from desired value. */
  61 double _dist_sq(vqgen *v,double *a, double *b){
  62   int i;
  63   int el=v->elements;
  64   double acc=0.;
  65   for(i=0;i<el;i++){
  66     double val=(a[i]-b[i]);
  67     acc+=val*val;
  68   }
  69   return acc;
  70 }
  71
  72 double *_weight_null(vqgen *v,double *a){
  73   return a;
  74 }
  75
  76 /* *must* be beefed up. */
  77 void _vqgen_seed(vqgen *v){
  78   long i;
  79   for(i=0;i<v->entries;i++)
  80     memcpy(_now(v,i),_point(v,i),sizeof(double)*v->elements);
  81 }
  82
  83 /* External calls *******************************************************/
  84
  85 void spinnit(char *s,int n){
  86   static int p=0;
  87   static long lasttime=0;
  88   long test;
  89   struct timeval thistime;
  90
  91   gettimeofday(&thistime,NULL);
  92   test=thistime.tv_sec*10+thistime.tv_usec/100000;
  93   if(lasttime!=test){
  94     lasttime=test;
  95
  96     fprintf(stderr,"%s%d ",s,n);
  97
  98     p++;if(p>3)p=0;
  99     switch(p){
 100     case 0:
 101       fprintf(stderr,"|    \r");
 102       break;
 103     case 1:
 104       fprintf(stderr,"/    \r");
 105       break;
 106     case 2:
 107       fprintf(stderr,"-    \r");
 108       break;
 109     case 3:
 110       fprintf(stderr,"\\    \r");
 111       break;
 112     }
 113     fflush(stderr);
 114   }
 115 }
 116
 117 /* We have two forms of quantization; in the first, each vector
 118    element in the codebook entry is orthogonal.  Residues would use this
 119    quantization for example.
 120
 121    In the second, we have a sequence of monotonically increasing
 122    values that we wish to quantize as deltas (to save space).  We
 123    still need to quantize so that absolute values are accurate. For
 124    example, LSP quantizes all absolute values, but the book encodes
 125    distance between values because each successive value is larger
 126    than the preceeding value.  Thus the desired quantibits apply to
 127    the encoded (delta) values, not abs positions. This requires minor
 128    additional encode-side trickery. */
 129
 130 void vqgen_quantize(vqgen *v,quant_meta *q){
 131
 132   double maxdel;
 133   double mindel;
 134
 135   double delta;
 136   double maxquant=((1<<q->quant)-1);
 137
 138   int j,k;
 139
 140   mindel=maxdel=_now(v,0)[0];
 141
 142   for(j=0;j<v->entries;j++){
 143     double last=0.;
 144     for(k=0;k<v->elements;k++){
 145       if(mindel>_now(v,j)[k]-last)mindel=_now(v,j)[k]-last;
 146       if(maxdel<_now(v,j)[k]-last)maxdel=_now(v,j)[k]-last;
 147       if(q->sequencep)last=_now(v,j)[k];
 148     }
 149   }
 150
 151
 152   /* first find the basic delta amount from the maximum span to be
 153      encoded.  Loosen the delta slightly to allow for additional error
 154      during sequence quantization */
 155
 156   delta=(maxdel-mindel)/((1<<q->quant)-1.5);
 157
 158   q->min=float24_pack(mindel);
 159   q->delta=float24_pack(delta);
 160
 161   for(j=0;j<v->entries;j++){
 162     double last=0;
 163     for(k=0;k<v->elements;k++){
 164       double val=_now(v,j)[k];
 165       double now=rint((val-last-mindel)/delta);
 166
 167       _now(v,j)[k]=now;
 168       if(now<0){
 169         /* be paranoid; this should be impossible */
 170         fprintf(stderr,"fault; quantized value<0\n");
 171         exit(1);
 172       }
 173
 174       if(now>maxquant){
 175         /* be paranoid; this should be impossible */
 176         fprintf(stderr,"fault; quantized value>max\n");
 177         exit(1);
 178       }
 179       if(q->sequencep)last=(now*delta)+mindel+last;
 180     }
 181   }
 182 }
 183
 184 /* much easier :-) */
 185 void vqgen_unquantize(vqgen *v,quant_meta *q){
 186   long j,k;
 187   double mindel=float24_unpack(q->min);
 188   double delta=float24_unpack(q->delta);
 189
 190   for(j=0;j<v->entries;j++){
 191     double last=0.;
 192     for(k=0;k<v->elements;k++){
 193       double now=_now(v,j)[k]*delta+last+mindel;
 194       _now(v,j)[k]=now;
 195       if(q->sequencep)last=now;
 196
 197     }
 198   }
 199 }
 200
 201 void vqgen_init(vqgen *v,int elements,int aux,int entries,
 202                 double  (*metric)(vqgen *,double *, double *),
 203                 double *(*weight)(vqgen *,double *)){
 204   memset(v,0,sizeof(vqgen));
 205
 206   v->elements=elements;
 207   v->aux=aux;
 208   v->allocated=32768;
 209   v->pointlist=malloc(v->allocated*(v->elements+v->aux)*sizeof(double));
 210
 211   v->entries=entries;
 212   v->entrylist=malloc(v->entries*v->elements*sizeof(double));
 213   v->assigned=malloc(v->entries*sizeof(long));
 214   v->bias=calloc(v->entries,sizeof(double));
 215   if(metric)
 216     v->metric_func=metric;
 217   else
 218     v->metric_func=_dist_sq;
 219   if(weight)
 220     v->weight_func=weight;
 221   else
 222     v->weight_func=_weight_null;
 223
 224 }
 225
 226 void vqgen_addpoint(vqgen *v, double *p,double *a){
 227   if(v->points>=v->allocated){
 228     v->allocated*=2;
 229     v->pointlist=realloc(v->pointlist,v->allocated*(v->elements+v->aux)*
 230                          sizeof(double));
 231   }
 232
 233   memcpy(_point(v,v->points),p,sizeof(double)*v->elements);
 234   if(v->aux)memcpy(_point(v,v->points)+v->elements,a,sizeof(double)*v->aux);
 235   v->points++;
 236   if(v->points==v->entries)_vqgen_seed(v);
 237 }
 238
 239 int directdsort(const void *a, const void *b){
 240   double av=*((double *)a);
 241   double bv=*((double *)b);
 242   if(av>bv)return(-1);
 243   return(1);
 244 }
 245
 246 double vqgen_iterate(vqgen *v){
 247   long   i,j,k;
 248   double fdesired=(double)v->points/v->entries;
 249   long  desired=fdesired;
 250   long  desired2=desired*2;
 251   double asserror=0.;
 252   double meterror=0.;
 253   double *new=malloc(sizeof(double)*v->entries*v->elements);
 254   long   *nearcount=malloc(v->entries*sizeof(long));
 255   double *nearbias=malloc(v->entries*desired2*sizeof(double));
 256
 257   if(v->entries<2){
 258     fprintf(stderr,"generation requires at least two entries\n");
 259     exit(1);
 260   }
 261
 262   /* fill in nearest points for entries */
 263   /*memset(v->bias,0,sizeof(double)*v->entries);*/
 264   memset(nearcount,0,sizeof(long)*v->entries);
 265   memset(v->assigned,0,sizeof(long)*v->entries);
 266   for(i=0;i<v->points;i++){
 267     double *ppt=v->weight_func(v,_point(v,i));
 268     double firstmetric=v->metric_func(v,_now(v,0),ppt)+v->bias[0];
 269     double secondmetric=v->metric_func(v,_now(v,1),ppt)+v->bias[1];
 270     long   firstentry=0;
 271     long   secondentry=1;
 272
 273     spinnit("working... ",v->points-i);
 274
 275     if(firstmetric>secondmetric){
 276       double temp=firstmetric;
 277       firstmetric=secondmetric;
 278       secondmetric=temp;
 279       firstentry=1;
 280       secondentry=0;
 281     }
 282
 283     for(j=2;j<v->entries;j++){
 284       double thismetric=v->metric_func(v,_now(v,j),ppt)+v->bias[j];
 285       if(thismetric<secondmetric){
 286         if(thismetric<firstmetric){
 287           secondmetric=firstmetric;
 288           secondentry=firstentry;
 289           firstmetric=thismetric;
 290           firstentry=j;
 291         }else{
 292           secondmetric=thismetric;
 293           secondentry=j;
 294         }
 295       }
 296     }
 297
 298     j=firstentry;
 299     meterror+=sqrt(_dist_sq(v,_now(v,j),ppt)/v->elements);
 300     /* set up midpoints for next iter */
 301     if(v->assigned[j]++)
 302       for(k=0;k<v->elements;k++)
 303         vN(new,j)[k]+=ppt[k];
 304     else
 305       for(k=0;k<v->elements;k++)
 306         vN(new,j)[k]=ppt[k];
 307
 308     for(j=0;j<v->entries;j++){
 309
 310       double thismetric;
 311       double *nearbiasptr=nearbias+desired2*j;
 312       long k=nearcount[j];
 313
 314       /* 'thismetric' is to be the bias value necessary in the current
 315          arrangement for entry j to capture point i */
 316       if(firstentry==j){
 317         /* use the secondary entry as the threshhold */
 318         thismetric=secondmetric-v->metric_func(v,_now(v,j),ppt);
 319       }else{
 320         /* use the primary entry as the threshhold */
 321         thismetric=firstmetric-v->metric_func(v,_now(v,j),ppt);
 322       }
 323
 324       /* a cute two-stage delayed sorting hack */
 325       if(k<desired){
 326         nearbiasptr[k]=thismetric;
 327         k++;
 328         if(k==desired)
 329           qsort(nearbiasptr,desired,sizeof(double),directdsort);
 330
 331       }else if(thismetric>nearbiasptr[desired-1]){
 332         nearbiasptr[k]=thismetric;
 333         k++;
 334         if(k==desired2){
 335           qsort(nearbiasptr,desired2,sizeof(double),directdsort);
 336           k=desired;
 337         }
 338       }
 339       nearcount[j]=k;
 340     }
 341   }
 342
 343   /* inflate/deflate */
 344   for(i=0;i<v->entries;i++){
 345     double *nearbiasptr=nearbias+desired2*i;
 346
 347     spinnit("working... ",v->entries-i);
 348
 349     /* due to the delayed sorting, we likely need to finish it off....*/
 350     if(nearcount[i]>desired)
 351       qsort(nearbiasptr,nearcount[i],sizeof(double),directdsort);
 352
 353     v->bias[i]=nearbiasptr[desired-1];
 354   }
 355
 356   /* assign midpoints */
 357
 358   for(j=0;j<v->entries;j++){
 359     asserror+=fabs(v->assigned[j]-fdesired);
 360     if(v->assigned[j])
 361       for(k=0;k<v->elements;k++)
 362         _now(v,j)[k]=vN(new,j)[k]/v->assigned[j];
 363   }
 364
 365   asserror/=(v->entries*fdesired);
 366
 367   fprintf(stderr,"Pass #%d... ",v->it);
 368   fprintf(stderr,": dist %g(%g) metric error=%g \n",
 369           asserror,fdesired,meterror/v->points);
 370   v->it++;
 371
 372   free(new);
 373   free(nearcount);
 374   free(nearbias);
 375 #ifdef NOISY
 376   fclose(assig);
 377   fclose(bias);
 378   fclose(cells);
 379 #endif
 380   return(asserror);
 381 }
 382