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