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.29 2000/02/16 16:18:40 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 int directdsort(const void *a, const void *b){
  84   double av=*((double *)a);
  85   double bv=*((double *)b);
  86   if(av>bv)return(-1);
  87   return(1);
  88 }
  89
  90 void vqgen_cellmetric(vqgen *v){
  91   int i,j,k;
  92   double min=-1.,max=-1.,mean=0.,acc=0.;
  93   long dup=0,unused=0;
  94  #ifdef NOISY
  95    char buff[80];
  96    double 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     double localmin=-1.;
 106
 107     for(k=0;k<v->entries;k++){
 108       if(j!=k){
 109         double 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 (%d unused/%d dup)\n",
 139           min,mean/acc,max,unused,dup);
 140
 141 #ifdef NOISY
 142   qsort(spacings,count,sizeof(double),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   double maxdel;
 168   double mindel;
 169
 170   double delta;
 171   double 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     double last=0.;
 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.5);
 192
 193   q->min=float24_pack(mindel);
 194   q->delta=float24_pack(delta);
 195
 196   for(j=0;j<v->entries;j++){
 197     double last=0;
 198     for(k=0;k<v->elements;k++){
 199       double val=_now(v,j)[k];
 200       double now=rint((val-last-mindel)/delta);
 201
 202       _now(v,j)[k]=now;
 203       if(now<0){
 204         /* be paranoid; this should be impossible */
 205         fprintf(stderr,"fault; quantized value<0\n");
 206         exit(1);
 207       }
 208
 209       if(now>maxquant){
 210         /* be paranoid; this should be impossible */
 211         fprintf(stderr,"fault; quantized value>max\n");
 212         exit(1);
 213       }
 214       if(q->sequencep)last=(now*delta)+mindel+last;
 215     }
 216   }
 217 }
 218
 219 /* much easier :-) */
 220 void vqgen_unquantize(vqgen *v,quant_meta *q){
 221   long j,k;
 222   double mindel=float24_unpack(q->min);
 223   double delta=float24_unpack(q->delta);
 224
 225   for(j=0;j<v->entries;j++){
 226     double last=0.;
 227     for(k=0;k<v->elements;k++){
 228       double now=_now(v,j)[k]*delta+last+mindel;
 229       _now(v,j)[k]=now;
 230       if(q->sequencep)last=now;
 231
 232     }
 233   }
 234 }
 235
 236 void vqgen_init(vqgen *v,int elements,int aux,int entries,double mindist,
 237                 double  (*metric)(vqgen *,double *, double *),
 238                 double *(*weight)(vqgen *,double *)){
 239   memset(v,0,sizeof(vqgen));
 240
 241   v->elements=elements;
 242   v->aux=aux;
 243   v->mindist=mindist;
 244   v->allocated=32768;
 245   v->pointlist=malloc(v->allocated*(v->elements+v->aux)*sizeof(double));
 246
 247   v->entries=entries;
 248   v->entrylist=malloc(v->entries*v->elements*sizeof(double));
 249   v->assigned=malloc(v->entries*sizeof(long));
 250   v->bias=calloc(v->entries,sizeof(double));
 251   v->max=calloc(v->entries,sizeof(double));
 252   if(metric)
 253     v->metric_func=metric;
 254   else
 255     v->metric_func=_dist;
 256   if(weight)
 257     v->weight_func=weight;
 258   else
 259     v->weight_func=_weight_null;
 260
 261 }
 262
 263 void vqgen_addpoint(vqgen *v, double *p,double *a){
 264   if(v->points>=v->allocated){
 265     v->allocated*=2;
 266     v->pointlist=realloc(v->pointlist,v->allocated*(v->elements+v->aux)*
 267                          sizeof(double));
 268   }
 269
 270   memcpy(_point(v,v->points),p,sizeof(double)*v->elements);
 271   if(v->aux)memcpy(_point(v,v->points)+v->elements,a,sizeof(double)*v->aux);
 272   v->points++;
 273   if(v->points==v->entries)_vqgen_seed(v);
 274   if(!(v->points&0xff))spinnit("loading... ",v->points);
 275 }
 276
 277 double vqgen_iterate(vqgen *v){
 278   long   i,j,k;
 279   long   biasable;
 280
 281   double fdesired=(double)v->points/v->entries;
 282   long  desired=fdesired;
 283   long  desired2=desired*2;
 284
 285   double asserror=0.;
 286   double meterror=0.;
 287   double *new=malloc(sizeof(double)*v->entries*v->elements);
 288   long   *nearcount=malloc(v->entries*sizeof(long));
 289   double *nearbias=malloc(v->entries*desired2*sizeof(double));
 290  #ifdef NOISY
 291    char buff[80];
 292    FILE *assig;
 293    FILE *bias;
 294    FILE *cells;
 295    sprintf(buff,"cells%d.m",v->it);
 296    cells=fopen(buff,"w");
 297    sprintf(buff,"assig%d.m",v->it);
 298    assig=fopen(buff,"w");
 299    sprintf(buff,"bias%d.m",v->it);
 300    bias=fopen(buff,"w");
 301  #endif
 302
 303
 304   if(v->entries<2){
 305     fprintf(stderr,"generation requires at least two entries\n");
 306     exit(1);
 307   }
 308
 309   /* fill in nearest points for entry biasing */
 310   /*memset(v->bias,0,sizeof(double)*v->entries);*/
 311   memset(nearcount,0,sizeof(long)*v->entries);
 312   memset(v->assigned,0,sizeof(long)*v->entries);
 313   biasable=0;
 314   for(i=0;i<v->points;i++){
 315     double *ppt=v->weight_func(v,_point(v,i));
 316     double firstmetric=v->metric_func(v,_now(v,0),ppt)+v->bias[0];
 317     double secondmetric=v->metric_func(v,_now(v,1),ppt)+v->bias[1];
 318     long   firstentry=0;
 319     long   secondentry=1;
 320     int    biasflag=1;
 321
 322     if(!(i&0xff))spinnit("biasing... ",v->points+v->points+v->entries-i);
 323
 324     if(firstmetric>secondmetric){
 325       double temp=firstmetric;
 326       firstmetric=secondmetric;
 327       secondmetric=temp;
 328       firstentry=1;
 329       secondentry=0;
 330     }
 331
 332     for(j=2;j<v->entries;j++){
 333       double thismetric=v->metric_func(v,_now(v,j),ppt)+v->bias[j];
 334       if(thismetric<secondmetric){
 335         if(thismetric<firstmetric){
 336           secondmetric=firstmetric;
 337           secondentry=firstentry;
 338           firstmetric=thismetric;
 339           firstentry=j;
 340         }else{
 341           secondmetric=thismetric;
 342           secondentry=j;
 343         }
 344       }
 345     }
 346
 347     j=firstentry;
 348     for(j=0;j<v->entries;j++){
 349
 350       double thismetric,localmetric;
 351       double *nearbiasptr=nearbias+desired2*j;
 352       long k=nearcount[j];
 353
 354       localmetric=v->metric_func(v,_now(v,j),ppt);
 355       /* 'thismetric' is to be the bias value necessary in the current
 356          arrangement for entry j to capture point i */
 357       if(firstentry==j){
 358         /* use the secondary entry as the threshhold */
 359         thismetric=secondmetric-localmetric;
 360       }else{
 361         /* use the primary entry as the threshhold */
 362         thismetric=firstmetric-localmetric;
 363       }
 364
 365       /* support the idea of 'minimum distance'... if we want the
 366          cells in a codebook to be roughly some minimum size (as with
 367          the low resolution residue books) */
 368
 369       if(localmetric>=v->mindist){
 370
 371         /* a cute two-stage delayed sorting hack */
 372         if(k<desired){
 373           nearbiasptr[k]=thismetric;
 374           k++;
 375           if(k==desired){
 376             spinnit("biasing... ",v->points+v->points+v->entries-i);
 377             qsort(nearbiasptr,desired,sizeof(double),directdsort);
 378           }
 379
 380         }else if(thismetric>nearbiasptr[desired-1]){
 381           nearbiasptr[k]=thismetric;
 382           k++;
 383           if(k==desired2){
 384             spinnit("biasing... ",v->points+v->points+v->entries-i);
 385             qsort(nearbiasptr,desired2,sizeof(double),directdsort);
 386             k=desired;
 387           }
 388         }
 389         nearcount[j]=k;
 390       }else
 391         biasflag=0;
 392     }
 393     biasable+=biasflag;
 394   }
 395
 396   /* inflate/deflate */
 397
 398   for(i=0;i<v->entries;i++){
 399     double *nearbiasptr=nearbias+desired2*i;
 400
 401     spinnit("biasing... ",v->points+v->entries-i);
 402
 403     /* due to the delayed sorting, we likely need to finish it off....*/
 404     if(nearcount[i]>desired)
 405       qsort(nearbiasptr,nearcount[i],sizeof(double),directdsort);
 406
 407     /* desired is the *maximum* bias queue size.  If we're using
 408        minimum distance, we're not interested in the max size... we're
 409        interested in the biasable number of points */
 410     {
 411       long localdesired=(double)biasable/v->entries;
 412       if(localdesired)
 413         v->bias[i]=nearbiasptr[localdesired-1];
 414       else
 415         v->bias[i]=nearbiasptr[0];
 416     }
 417   }
 418
 419   /* Now assign with new bias and find new midpoints */
 420   for(i=0;i<v->points;i++){
 421     double *ppt=v->weight_func(v,_point(v,i));
 422     double firstmetric=v->metric_func(v,_now(v,0),ppt)+v->bias[0];
 423     long   firstentry=0;
 424
 425     if(!(i&0xff))spinnit("centering... ",v->points-i);
 426
 427     for(j=0;j<v->entries;j++){
 428       double thismetric=v->metric_func(v,_now(v,j),ppt)+v->bias[j];
 429       if(thismetric<firstmetric){
 430         firstmetric=thismetric;
 431         firstentry=j;
 432       }
 433     }
 434
 435     j=firstentry;
 436
 437 #ifdef NOISY
 438     fprintf(cells,"%g %g\n%g %g\n\n",
 439           _now(v,j)[0],_now(v,j)[1],
 440           ppt[0],ppt[1]);
 441 #endif
 442
 443     firstmetric-=v->bias[firstentry];
 444     meterror+=firstmetric;
 445     /* set up midpoints for next iter */
 446     if(v->assigned[j]++){
 447       for(k=0;k<v->elements;k++)
 448         vN(new,j)[k]+=ppt[k];
 449       if(firstmetric>v->max[firstentry])v->max[firstentry]=firstmetric;
 450     }else{
 451       for(k=0;k<v->elements;k++)
 452         vN(new,j)[k]=ppt[k];
 453       v->max[firstentry]=firstmetric;
 454     }
 455   }
 456
 457   /* assign midpoints */
 458
 459   for(j=0;j<v->entries;j++){
 460 #ifdef NOISY
 461     fprintf(assig,"%ld\n",v->assigned[j]);
 462     fprintf(bias,"%g\n",v->bias[j]);
 463 #endif
 464     asserror+=fabs(v->assigned[j]-fdesired);
 465     if(v->assigned[j])
 466       for(k=0;k<v->elements;k++)
 467         _now(v,j)[k]=vN(new,j)[k]/v->assigned[j];
 468   }
 469
 470   asserror/=(v->entries*fdesired);
 471
 472   fprintf(stderr,"Pass #%d... ",v->it);
 473   fprintf(stderr,": dist %g(%g) metric error=%g \n",
 474           asserror,fdesired,meterror/v->points);
 475   v->it++;
 476
 477   free(new);
 478   free(nearcount);
 479   free(nearbias);
 480 #ifdef NOISY
 481   fclose(assig);
 482   fclose(bias);
 483   fclose(cells);
 484 #endif
 485   return(asserror);
 486 }
 487