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.27 2000/01/05 15:05:01 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(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 sqrt(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;
 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  #ifdef NOISY
 257    char buff[80];
 258    FILE *assig;
 259    FILE *bias;
 260    FILE *cells;
 261    sprintf(buff,"cells%d.m",v->it);
 262    cells=fopen(buff,"w");
 263    sprintf(buff,"assig%d.m",v->it);
 264    assig=fopen(buff,"w");
 265    sprintf(buff,"bias%d.m",v->it);
 266    bias=fopen(buff,"w");
 267  #endif
 268
 269
 270   if(v->entries<2){
 271     fprintf(stderr,"generation requires at least two entries\n");
 272     exit(1);
 273   }
 274
 275   /* fill in nearest points for entries */
 276   /*memset(v->bias,0,sizeof(double)*v->entries);*/
 277   memset(nearcount,0,sizeof(long)*v->entries);
 278   memset(v->assigned,0,sizeof(long)*v->entries);
 279   for(i=0;i<v->points;i++){
 280     double *ppt=v->weight_func(v,_point(v,i));
 281     double firstmetric=v->metric_func(v,_now(v,0),ppt)+v->bias[0];
 282     double secondmetric=v->metric_func(v,_now(v,1),ppt)+v->bias[1];
 283     long   firstentry=0;
 284     long   secondentry=1;
 285
 286     spinnit("working... ",v->points-i);
 287
 288     if(firstmetric>secondmetric){
 289       double temp=firstmetric;
 290       firstmetric=secondmetric;
 291       secondmetric=temp;
 292       firstentry=1;
 293       secondentry=0;
 294     }
 295
 296     for(j=2;j<v->entries;j++){
 297       double thismetric=v->metric_func(v,_now(v,j),ppt)+v->bias[j];
 298       if(thismetric<secondmetric){
 299         if(thismetric<firstmetric){
 300           secondmetric=firstmetric;
 301           secondentry=firstentry;
 302           firstmetric=thismetric;
 303           firstentry=j;
 304         }else{
 305           secondmetric=thismetric;
 306           secondentry=j;
 307         }
 308       }
 309     }
 310
 311     j=firstentry;
 312
 313 #ifdef NOISY
 314     fprintf(cells,"%g %g\n%g %g\n\n",
 315           _now(v,j)[0],_now(v,j)[1],
 316           ppt[0],ppt[1]);
 317 #endif
 318
 319     j=firstentry;
 320     meterror+=firstmetric-v->bias[firstentry];
 321     /* set up midpoints for next iter */
 322     if(v->assigned[j]++)
 323       for(k=0;k<v->elements;k++)
 324         vN(new,j)[k]+=ppt[k];
 325     else
 326       for(k=0;k<v->elements;k++)
 327         vN(new,j)[k]=ppt[k];
 328
 329     for(j=0;j<v->entries;j++){
 330
 331       double thismetric;
 332       double *nearbiasptr=nearbias+desired2*j;
 333       long k=nearcount[j];
 334
 335       /* 'thismetric' is to be the bias value necessary in the current
 336          arrangement for entry j to capture point i */
 337       if(firstentry==j){
 338         /* use the secondary entry as the threshhold */
 339         thismetric=secondmetric-v->metric_func(v,_now(v,j),ppt);
 340       }else{
 341         /* use the primary entry as the threshhold */
 342         thismetric=firstmetric-v->metric_func(v,_now(v,j),ppt);
 343       }
 344
 345       /* a cute two-stage delayed sorting hack */
 346       if(k<desired){
 347         nearbiasptr[k]=thismetric;
 348         k++;
 349         if(k==desired)
 350           qsort(nearbiasptr,desired,sizeof(double),directdsort);
 351
 352       }else if(thismetric>nearbiasptr[desired-1]){
 353         nearbiasptr[k]=thismetric;
 354         k++;
 355         if(k==desired2){
 356           qsort(nearbiasptr,desired2,sizeof(double),directdsort);
 357           k=desired;
 358         }
 359       }
 360       nearcount[j]=k;
 361     }
 362   }
 363
 364   /* inflate/deflate */
 365   for(i=0;i<v->entries;i++){
 366     double *nearbiasptr=nearbias+desired2*i;
 367
 368     spinnit("working... ",v->entries-i);
 369
 370     /* due to the delayed sorting, we likely need to finish it off....*/
 371     if(nearcount[i]>desired)
 372       qsort(nearbiasptr,nearcount[i],sizeof(double),directdsort);
 373
 374     v->bias[i]=nearbiasptr[desired-1];
 375   }
 376
 377   /* assign midpoints */
 378
 379   for(j=0;j<v->entries;j++){
 380 #ifdef NOISY
 381     fprintf(assig,"%ld\n",v->assigned[j]);
 382     fprintf(bias,"%g\n",v->bias[j]);
 383 #endif
 384     asserror+=fabs(v->assigned[j]-fdesired);
 385     if(v->assigned[j])
 386       for(k=0;k<v->elements;k++)
 387         _now(v,j)[k]=vN(new,j)[k]/v->assigned[j];
 388   }
 389
 390   asserror/=(v->entries*fdesired);
 391
 392   fprintf(stderr,"Pass #%d... ",v->it);
 393   fprintf(stderr,": dist %g(%g) metric error=%g \n",
 394           asserror,fdesired,meterror/v->points);
 395   v->it++;
 396
 397   free(new);
 398   free(nearcount);
 399   free(nearbias);
 400 #ifdef NOISY
 401   fclose(assig);
 402   fclose(bias);
 403   fclose(cells);
 404 #endif
 405   return(asserror);
 406 }
 407