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.28 2000/01/06 13:57:13 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 /* We have two forms of quantization; in the first, each vector
  86    element in the codebook entry is orthogonal.  Residues would use this
  87    quantization for example.
  88
  89    In the second, we have a sequence of monotonically increasing
  90    values that we wish to quantize as deltas (to save space).  We
  91    still need to quantize so that absolute values are accurate. For
  92    example, LSP quantizes all absolute values, but the book encodes
  93    distance between values because each successive value is larger
  94    than the preceeding value.  Thus the desired quantibits apply to
  95    the encoded (delta) values, not abs positions. This requires minor
  96    additional encode-side trickery. */
  97
  98 void vqgen_quantize(vqgen *v,quant_meta *q){
  99
 100   double maxdel;
 101   double mindel;
 102
 103   double delta;
 104   double maxquant=((1<<q->quant)-1);
 105
 106   int j,k;
 107
 108   mindel=maxdel=_now(v,0)[0];
 109
 110   for(j=0;j<v->entries;j++){
 111     double last=0.;
 112     for(k=0;k<v->elements;k++){
 113       if(mindel>_now(v,j)[k]-last)mindel=_now(v,j)[k]-last;
 114       if(maxdel<_now(v,j)[k]-last)maxdel=_now(v,j)[k]-last;
 115       if(q->sequencep)last=_now(v,j)[k];
 116     }
 117   }
 118
 119
 120   /* first find the basic delta amount from the maximum span to be
 121      encoded.  Loosen the delta slightly to allow for additional error
 122      during sequence quantization */
 123
 124   delta=(maxdel-mindel)/((1<<q->quant)-1.5);
 125
 126   q->min=float24_pack(mindel);
 127   q->delta=float24_pack(delta);
 128
 129   for(j=0;j<v->entries;j++){
 130     double last=0;
 131     for(k=0;k<v->elements;k++){
 132       double val=_now(v,j)[k];
 133       double now=rint((val-last-mindel)/delta);
 134
 135       _now(v,j)[k]=now;
 136       if(now<0){
 137         /* be paranoid; this should be impossible */
 138         fprintf(stderr,"fault; quantized value<0\n");
 139         exit(1);
 140       }
 141
 142       if(now>maxquant){
 143         /* be paranoid; this should be impossible */
 144         fprintf(stderr,"fault; quantized value>max\n");
 145         exit(1);
 146       }
 147       if(q->sequencep)last=(now*delta)+mindel+last;
 148     }
 149   }
 150 }
 151
 152 /* much easier :-) */
 153 void vqgen_unquantize(vqgen *v,quant_meta *q){
 154   long j,k;
 155   double mindel=float24_unpack(q->min);
 156   double delta=float24_unpack(q->delta);
 157
 158   for(j=0;j<v->entries;j++){
 159     double last=0.;
 160     for(k=0;k<v->elements;k++){
 161       double now=_now(v,j)[k]*delta+last+mindel;
 162       _now(v,j)[k]=now;
 163       if(q->sequencep)last=now;
 164
 165     }
 166   }
 167 }
 168
 169 void vqgen_init(vqgen *v,int elements,int aux,int entries,
 170                 double  (*metric)(vqgen *,double *, double *),
 171                 double *(*weight)(vqgen *,double *)){
 172   memset(v,0,sizeof(vqgen));
 173
 174   v->elements=elements;
 175   v->aux=aux;
 176   v->allocated=32768;
 177   v->pointlist=malloc(v->allocated*(v->elements+v->aux)*sizeof(double));
 178
 179   v->entries=entries;
 180   v->entrylist=malloc(v->entries*v->elements*sizeof(double));
 181   v->assigned=malloc(v->entries*sizeof(long));
 182   v->bias=calloc(v->entries,sizeof(double));
 183   if(metric)
 184     v->metric_func=metric;
 185   else
 186     v->metric_func=_dist;
 187   if(weight)
 188     v->weight_func=weight;
 189   else
 190     v->weight_func=_weight_null;
 191
 192 }
 193
 194 void vqgen_addpoint(vqgen *v, double *p,double *a){
 195   if(v->points>=v->allocated){
 196     v->allocated*=2;
 197     v->pointlist=realloc(v->pointlist,v->allocated*(v->elements+v->aux)*
 198                          sizeof(double));
 199   }
 200
 201   memcpy(_point(v,v->points),p,sizeof(double)*v->elements);
 202   if(v->aux)memcpy(_point(v,v->points)+v->elements,a,sizeof(double)*v->aux);
 203   v->points++;
 204   if(v->points==v->entries)_vqgen_seed(v);
 205 }
 206
 207 int directdsort(const void *a, const void *b){
 208   double av=*((double *)a);
 209   double bv=*((double *)b);
 210   if(av>bv)return(-1);
 211   return(1);
 212 }
 213
 214 double vqgen_iterate(vqgen *v){
 215   long   i,j,k;
 216   double fdesired=(double)v->points/v->entries;
 217   long  desired=fdesired;
 218   long  desired2=desired*2;
 219   double asserror=0.;
 220   double meterror=0.;
 221   double *new=malloc(sizeof(double)*v->entries*v->elements);
 222   long   *nearcount=malloc(v->entries*sizeof(long));
 223   double *nearbias=malloc(v->entries*desired2*sizeof(double));
 224  #ifdef NOISY
 225    char buff[80];
 226    FILE *assig;
 227    FILE *bias;
 228    FILE *cells;
 229    sprintf(buff,"cells%d.m",v->it);
 230    cells=fopen(buff,"w");
 231    sprintf(buff,"assig%d.m",v->it);
 232    assig=fopen(buff,"w");
 233    sprintf(buff,"bias%d.m",v->it);
 234    bias=fopen(buff,"w");
 235  #endif
 236
 237
 238   if(v->entries<2){
 239     fprintf(stderr,"generation requires at least two entries\n");
 240     exit(1);
 241   }
 242
 243   /* fill in nearest points for entry biasing */
 244   /*memset(v->bias,0,sizeof(double)*v->entries);*/
 245   memset(nearcount,0,sizeof(long)*v->entries);
 246   memset(v->assigned,0,sizeof(long)*v->entries);
 247   for(i=0;i<v->points;i++){
 248     double *ppt=v->weight_func(v,_point(v,i));
 249     double firstmetric=v->metric_func(v,_now(v,0),ppt)+v->bias[0];
 250     double secondmetric=v->metric_func(v,_now(v,1),ppt)+v->bias[1];
 251     long   firstentry=0;
 252     long   secondentry=1;
 253
 254     if(i%100)spinnit("biasing... ",v->points+v->points+v->entries-i);
 255
 256     if(firstmetric>secondmetric){
 257       double temp=firstmetric;
 258       firstmetric=secondmetric;
 259       secondmetric=temp;
 260       firstentry=1;
 261       secondentry=0;
 262     }
 263
 264     for(j=2;j<v->entries;j++){
 265       double thismetric=v->metric_func(v,_now(v,j),ppt)+v->bias[j];
 266       if(thismetric<secondmetric){
 267         if(thismetric<firstmetric){
 268           secondmetric=firstmetric;
 269           secondentry=firstentry;
 270           firstmetric=thismetric;
 271           firstentry=j;
 272         }else{
 273           secondmetric=thismetric;
 274           secondentry=j;
 275         }
 276       }
 277     }
 278
 279     j=firstentry;
 280     for(j=0;j<v->entries;j++){
 281
 282       double thismetric;
 283       double *nearbiasptr=nearbias+desired2*j;
 284       long k=nearcount[j];
 285
 286       /* 'thismetric' is to be the bias value necessary in the current
 287          arrangement for entry j to capture point i */
 288       if(firstentry==j){
 289         /* use the secondary entry as the threshhold */
 290         thismetric=secondmetric-v->metric_func(v,_now(v,j),ppt);
 291       }else{
 292         /* use the primary entry as the threshhold */
 293         thismetric=firstmetric-v->metric_func(v,_now(v,j),ppt);
 294       }
 295
 296       /* a cute two-stage delayed sorting hack */
 297       if(k<desired){
 298         nearbiasptr[k]=thismetric;
 299         k++;
 300         if(k==desired){
 301           spinnit("biasing... ",v->points+v->points+v->entries-i);
 302           qsort(nearbiasptr,desired,sizeof(double),directdsort);
 303         }
 304
 305       }else if(thismetric>nearbiasptr[desired-1]){
 306         nearbiasptr[k]=thismetric;
 307         k++;
 308         if(k==desired2){
 309           spinnit("biasing... ",v->points+v->points+v->entries-i);
 310           qsort(nearbiasptr,desired2,sizeof(double),directdsort);
 311           k=desired;
 312         }
 313       }
 314       nearcount[j]=k;
 315     }
 316   }
 317
 318   /* inflate/deflate */
 319   for(i=0;i<v->entries;i++){
 320     double *nearbiasptr=nearbias+desired2*i;
 321
 322     spinnit("biasing... ",v->points+v->entries-i);
 323
 324     /* due to the delayed sorting, we likely need to finish it off....*/
 325     if(nearcount[i]>desired)
 326       qsort(nearbiasptr,nearcount[i],sizeof(double),directdsort);
 327
 328     v->bias[i]=nearbiasptr[desired-1];
 329   }
 330
 331   /* Now assign with new bias and find new midpoints */
 332   for(i=0;i<v->points;i++){
 333     double *ppt=v->weight_func(v,_point(v,i));
 334     double firstmetric=v->metric_func(v,_now(v,0),ppt)+v->bias[0];
 335     long   firstentry=0;
 336
 337     if(i%100)spinnit("centering... ",v->points-i);
 338
 339     for(j=0;j<v->entries;j++){
 340       double thismetric=v->metric_func(v,_now(v,j),ppt)+v->bias[j];
 341       if(thismetric<firstmetric){
 342         firstmetric=thismetric;
 343         firstentry=j;
 344       }
 345     }
 346
 347     j=firstentry;
 348
 349 #ifdef NOISY
 350     fprintf(cells,"%g %g\n%g %g\n\n",
 351           _now(v,j)[0],_now(v,j)[1],
 352           ppt[0],ppt[1]);
 353 #endif
 354
 355     meterror+=firstmetric-v->bias[firstentry];
 356     /* set up midpoints for next iter */
 357     if(v->assigned[j]++)
 358       for(k=0;k<v->elements;k++)
 359         vN(new,j)[k]+=ppt[k];
 360     else
 361       for(k=0;k<v->elements;k++)
 362         vN(new,j)[k]=ppt[k];
 363
 364   }
 365
 366   /* assign midpoints */
 367
 368   for(j=0;j<v->entries;j++){
 369 #ifdef NOISY
 370     fprintf(assig,"%ld\n",v->assigned[j]);
 371     fprintf(bias,"%g\n",v->bias[j]);
 372 #endif
 373     asserror+=fabs(v->assigned[j]-fdesired);
 374     if(v->assigned[j])
 375       for(k=0;k<v->elements;k++)
 376         _now(v,j)[k]=vN(new,j)[k]/v->assigned[j];
 377   }
 378
 379   asserror/=(v->entries*fdesired);
 380
 381   fprintf(stderr,"Pass #%d... ",v->it);
 382   fprintf(stderr,": dist %g(%g) metric error=%g \n",
 383           asserror,fdesired,meterror/v->points);
 384   v->it++;
 385
 386   free(new);
 387   free(nearcount);
 388   free(nearbias);
 389 #ifdef NOISY
 390   fclose(assig);
 391   fclose(bias);
 392   fclose(cells);
 393 #endif
 394   return(asserror);
 395 }
 396