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.31 2000/05/08 20:49:51 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 #include "../lib/sharedbook.h"
  36
  37 /* Codebook generation happens in two steps:
  38
  39    1) Train the codebook with data collected from the encoder: We use
  40    one of a few error metrics (which represent the distance between a
  41    given data point and a candidate point in the training set) to
  42    divide the training set up into cells representing roughly equal
  43    probability of occurring.
  44
  45    2) Generate the codebook and auxiliary data from the trained data set
  46 */
  47
  48 /* Codebook training ****************************************************
  49  *
  50  * The basic idea here is that a VQ codebook is like an m-dimensional
  51  * foam with n bubbles.  The bubbles compete for space/volume and are
  52  * 'pressurized' [biased] according to some metric.  The basic alg
  53  * iterates through allowing the bubbles to compete for space until
  54  * they converge (if the damping is dome properly) on a steady-state
  55  * solution. Individual input points, collected from libvorbis, are
  56  * used to train the algorithm monte-carlo style.  */
  57
  58 /* internal helpers *****************************************************/
  59 #define vN(data,i) (data+v->elements*i)
  60
  61 /* default metric; squared 'distance' from desired value. */
  62 double _dist(vqgen *v,double *a, double *b){
  63   int i;
  64   int el=v->elements;
  65   double acc=0.;
  66   for(i=0;i<el;i++){
  67     double val=(a[i]-b[i]);
  68     acc+=val*val;
  69   }
  70   return sqrt(acc);
  71 }
  72
  73 double *_weight_null(vqgen *v,double *a){
  74   return a;
  75 }
  76
  77 /* *must* be beefed up. */
  78 void _vqgen_seed(vqgen *v){
  79   long i;
  80   for(i=0;i<v->entries;i++)
  81     memcpy(_now(v,i),_point(v,i),sizeof(double)*v->elements);
  82 }
  83
  84 int directdsort(const void *a, const void *b){
  85   double av=*((double *)a);
  86   double bv=*((double *)b);
  87   if(av>bv)return(-1);
  88   return(1);
  89 }
  90
  91 void vqgen_cellmetric(vqgen *v){
  92   int j,k;
  93   double min=-1.,max=-1.,mean=0.,acc=0.;
  94   long dup=0,unused=0;
  95  #ifdef NOISY
  96   int i;
  97    char buff[80];
  98    double spacings[v->entries];
  99    int count=0;
 100    FILE *cells;
 101    sprintf(buff,"cellspace%d.m",v->it);
 102    cells=fopen(buff,"w");
 103 #endif
 104
 105   /* minimum, maximum, cell spacing */
 106   for(j=0;j<v->entries;j++){
 107     double localmin=-1.;
 108
 109     for(k=0;k<v->entries;k++){
 110       if(j!=k){
 111         double this=_dist(v,_now(v,j),_now(v,k));
 112         if(this>0){
 113           if(v->assigned[k] && (localmin==-1 || this<localmin))
 114             localmin=this;
 115         }else{
 116           if(k<j){
 117             dup++;
 118             break;
 119           }
 120         }
 121       }
 122     }
 123     if(k<v->entries)continue;
 124
 125     if(v->assigned[j]==0){
 126       unused++;
 127       continue;
 128     }
 129
 130     localmin=v->max[j]+localmin/2; /* this gives us rough diameter */
 131     if(min==-1 || localmin<min)min=localmin;
 132     if(max==-1 || localmin>max)max=localmin;
 133     mean+=localmin;
 134     acc++;
 135 #ifdef NOISY
 136     spacings[count++]=localmin;
 137 #endif
 138   }
 139
 140   fprintf(stderr,"cell diameter: %.03g::%.03g::%.03g (%ld unused/%ld dup)\n",
 141           min,mean/acc,max,unused,dup);
 142
 143 #ifdef NOISY
 144   qsort(spacings,count,sizeof(double),directdsort);
 145   for(i=0;i<count;i++)
 146     fprintf(cells,"%g\n",spacings[i]);
 147   fclose(cells);
 148 #endif
 149
 150 }
 151
 152 /* External calls *******************************************************/
 153
 154 /* We have two forms of quantization; in the first, each vector
 155    element in the codebook entry is orthogonal.  Residues would use this
 156    quantization for example.
 157
 158    In the second, we have a sequence of monotonically increasing
 159    values that we wish to quantize as deltas (to save space).  We
 160    still need to quantize so that absolute values are accurate. For
 161    example, LSP quantizes all absolute values, but the book encodes
 162    distance between values because each successive value is larger
 163    than the preceeding value.  Thus the desired quantibits apply to
 164    the encoded (delta) values, not abs positions. This requires minor
 165    additional encode-side trickery. */
 166
 167 void vqgen_quantize(vqgen *v,quant_meta *q){
 168
 169   double maxdel;
 170   double mindel;
 171
 172   double delta;
 173   double maxquant=((1<<q->quant)-1);
 174
 175   int j,k;
 176
 177   mindel=maxdel=_now(v,0)[0];
 178
 179   for(j=0;j<v->entries;j++){
 180     double last=0.;
 181     for(k=0;k<v->elements;k++){
 182       if(mindel>_now(v,j)[k]-last)mindel=_now(v,j)[k]-last;
 183       if(maxdel<_now(v,j)[k]-last)maxdel=_now(v,j)[k]-last;
 184       if(q->sequencep)last=_now(v,j)[k];
 185     }
 186   }
 187
 188
 189   /* first find the basic delta amount from the maximum span to be
 190      encoded.  Loosen the delta slightly to allow for additional error
 191      during sequence quantization */
 192
 193   delta=(maxdel-mindel)/((1<<q->quant)-1.5);
 194
 195   q->min=_float32_pack(mindel);
 196   q->delta=_float32_pack(delta);
 197
 198   mindel=_float32_unpack(q->min);
 199   delta=_float32_unpack(q->delta);
 200
 201   for(j=0;j<v->entries;j++){
 202     double last=0;
 203     for(k=0;k<v->elements;k++){
 204       double val=_now(v,j)[k];
 205       double now=rint((val-last-mindel)/delta);
 206
 207       _now(v,j)[k]=now;
 208       if(now<0){
 209         /* be paranoid; this should be impossible */
 210         fprintf(stderr,"fault; quantized value<0\n");
 211         exit(1);
 212       }
 213
 214       if(now>maxquant){
 215         /* be paranoid; this should be impossible */
 216         fprintf(stderr,"fault; quantized value>max\n");
 217         exit(1);
 218       }
 219       if(q->sequencep)last=(now*delta)+mindel+last;
 220     }
 221   }
 222 }
 223
 224 /* much easier :-).  Unlike in the codebook, we don't un-log log
 225    scales; we just make sure they're properly offset. */
 226 void vqgen_unquantize(vqgen *v,quant_meta *q){
 227   long j,k;
 228   double mindel=_float32_unpack(q->min);
 229   double delta=_float32_unpack(q->delta);
 230
 231   for(j=0;j<v->entries;j++){
 232     double last=0.;
 233     for(k=0;k<v->elements;k++){
 234       double now=_now(v,j)[k];
 235       now=fabs(now)*delta+last+mindel;
 236       if(q->sequencep)last=now;
 237       _now(v,j)[k]=now;
 238     }
 239   }
 240 }
 241
 242 void vqgen_init(vqgen *v,int elements,int aux,int entries,double mindist,
 243                 double  (*metric)(vqgen *,double *, double *),
 244                 double *(*weight)(vqgen *,double *)){
 245   memset(v,0,sizeof(vqgen));
 246
 247   v->elements=elements;
 248   v->aux=aux;
 249   v->mindist=mindist;
 250   v->allocated=32768;
 251   v->pointlist=malloc(v->allocated*(v->elements+v->aux)*sizeof(double));
 252
 253   v->entries=entries;
 254   v->entrylist=malloc(v->entries*v->elements*sizeof(double));
 255   v->assigned=malloc(v->entries*sizeof(long));
 256   v->bias=calloc(v->entries,sizeof(double));
 257   v->max=calloc(v->entries,sizeof(double));
 258   if(metric)
 259     v->metric_func=metric;
 260   else
 261     v->metric_func=_dist;
 262   if(weight)
 263     v->weight_func=weight;
 264   else
 265     v->weight_func=_weight_null;
 266
 267 }
 268
 269 void vqgen_addpoint(vqgen *v, double *p,double *a){
 270   if(v->points>=v->allocated){
 271     v->allocated*=2;
 272     v->pointlist=realloc(v->pointlist,v->allocated*(v->elements+v->aux)*
 273                          sizeof(double));
 274   }
 275
 276   memcpy(_point(v,v->points),p,sizeof(double)*v->elements);
 277   if(v->aux)memcpy(_point(v,v->points)+v->elements,a,sizeof(double)*v->aux);
 278   v->points++;
 279   if(v->points==v->entries)_vqgen_seed(v);
 280   if(!(v->points&0xff))spinnit("loading... ",v->points);
 281 }
 282
 283 double vqgen_iterate(vqgen *v,int biasp){
 284   long   i,j,k;
 285   long   biasable;
 286
 287   double fdesired=(double)v->points/v->entries;
 288   long  desired=fdesired;
 289   long  desired2=desired*2;
 290
 291   double asserror=0.;
 292   double meterror=0.;
 293   double *new=malloc(sizeof(double)*v->entries*v->elements);
 294   long   *nearcount=malloc(v->entries*sizeof(long));
 295   double *nearbias=malloc(v->entries*desired2*sizeof(double));
 296  #ifdef NOISY
 297    char buff[80];
 298    FILE *assig;
 299    FILE *bias;
 300    FILE *cells;
 301    sprintf(buff,"cells%d.m",v->it);
 302    cells=fopen(buff,"w");
 303    sprintf(buff,"assig%d.m",v->it);
 304    assig=fopen(buff,"w");
 305    sprintf(buff,"bias%d.m",v->it);
 306    bias=fopen(buff,"w");
 307  #endif
 308
 309
 310   if(v->entries<2){
 311     fprintf(stderr,"generation requires at least two entries\n");
 312     exit(1);
 313   }
 314
 315   /* fill in nearest points for entry biasing */
 316   /*memset(v->bias,0,sizeof(double)*v->entries);*/
 317   memset(nearcount,0,sizeof(long)*v->entries);
 318   memset(v->assigned,0,sizeof(long)*v->entries);
 319   biasable=0;
 320   if(biasp){
 321     for(i=0;i<v->points;i++){
 322       double *ppt=v->weight_func(v,_point(v,i));
 323       double firstmetric=v->metric_func(v,_now(v,0),ppt)+v->bias[0];
 324       double secondmetric=v->metric_func(v,_now(v,1),ppt)+v->bias[1];
 325       long   firstentry=0;
 326       long   secondentry=1;
 327       int    biasflag=1;
 328
 329       if(!(i&0xff))spinnit("biasing... ",v->points+v->points+v->entries-i);
 330
 331       if(firstmetric>secondmetric){
 332         double temp=firstmetric;
 333         firstmetric=secondmetric;
 334         secondmetric=temp;
 335         firstentry=1;
 336         secondentry=0;
 337       }
 338
 339       for(j=2;j<v->entries;j++){
 340         double thismetric=v->metric_func(v,_now(v,j),ppt)+v->bias[j];
 341         if(thismetric<secondmetric){
 342           if(thismetric<firstmetric){
 343             secondmetric=firstmetric;
 344             secondentry=firstentry;
 345             firstmetric=thismetric;
 346             firstentry=j;
 347           }else{
 348             secondmetric=thismetric;
 349             secondentry=j;
 350           }
 351         }
 352       }
 353
 354       j=firstentry;
 355       for(j=0;j<v->entries;j++){
 356
 357         double thismetric,localmetric;
 358         double *nearbiasptr=nearbias+desired2*j;
 359         long k=nearcount[j];
 360
 361         localmetric=v->metric_func(v,_now(v,j),ppt);
 362         /* 'thismetric' is to be the bias value necessary in the current
 363            arrangement for entry j to capture point i */
 364         if(firstentry==j){
 365           /* use the secondary entry as the threshhold */
 366           thismetric=secondmetric-localmetric;
 367         }else{
 368           /* use the primary entry as the threshhold */
 369           thismetric=firstmetric-localmetric;
 370         }
 371
 372         /* support the idea of 'minimum distance'... if we want the
 373            cells in a codebook to be roughly some minimum size (as with
 374            the low resolution residue books) */
 375
 376         if(localmetric>=v->mindist){
 377
 378           /* a cute two-stage delayed sorting hack */
 379           if(k<desired){
 380             nearbiasptr[k]=thismetric;
 381             k++;
 382             if(k==desired){
 383               spinnit("biasing... ",v->points+v->points+v->entries-i);
 384               qsort(nearbiasptr,desired,sizeof(double),directdsort);
 385             }
 386
 387           }else if(thismetric>nearbiasptr[desired-1]){
 388             nearbiasptr[k]=thismetric;
 389             k++;
 390             if(k==desired2){
 391               spinnit("biasing... ",v->points+v->points+v->entries-i);
 392               qsort(nearbiasptr,desired2,sizeof(double),directdsort);
 393               k=desired;
 394             }
 395           }
 396           nearcount[j]=k;
 397         }else
 398           biasflag=0;
 399       }
 400       biasable+=biasflag;
 401     }
 402
 403     /* inflate/deflate */
 404
 405     for(i=0;i<v->entries;i++){
 406       double *nearbiasptr=nearbias+desired2*i;
 407
 408       spinnit("biasing... ",v->points+v->entries-i);
 409
 410       /* due to the delayed sorting, we likely need to finish it off....*/
 411       if(nearcount[i]>desired)
 412         qsort(nearbiasptr,nearcount[i],sizeof(double),directdsort);
 413
 414       /* desired is the *maximum* bias queue size.  If we're using
 415          minimum distance, we're not interested in the max size... we're
 416          interested in the biasable number of points */
 417       {
 418         long localdesired=(double)biasable/v->entries;
 419         if(localdesired)
 420           v->bias[i]=nearbiasptr[localdesired-1];
 421         else
 422           v->bias[i]=nearbiasptr[0];
 423       }
 424     }
 425   }else{
 426     memset(v->bias,0,v->entries*sizeof(double));
 427   }
 428
 429   /* Now assign with new bias and find new midpoints */
 430   for(i=0;i<v->points;i++){
 431     double *ppt=v->weight_func(v,_point(v,i));
 432     double firstmetric=v->metric_func(v,_now(v,0),ppt)+v->bias[0];
 433     long   firstentry=0;
 434
 435     if(!(i&0xff))spinnit("centering... ",v->points-i);
 436
 437     for(j=0;j<v->entries;j++){
 438       double thismetric=v->metric_func(v,_now(v,j),ppt)+v->bias[j];
 439       if(thismetric<firstmetric){
 440         firstmetric=thismetric;
 441         firstentry=j;
 442       }
 443     }
 444
 445     j=firstentry;
 446
 447 #ifdef NOISY
 448     fprintf(cells,"%g %g\n%g %g\n\n",
 449           _now(v,j)[0],_now(v,j)[1],
 450           ppt[0],ppt[1]);
 451 #endif
 452
 453     firstmetric-=v->bias[firstentry];
 454     meterror+=firstmetric;
 455     /* set up midpoints for next iter */
 456     if(v->assigned[j]++){
 457       for(k=0;k<v->elements;k++)
 458         vN(new,j)[k]+=ppt[k];
 459       if(firstmetric>v->max[firstentry])v->max[firstentry]=firstmetric;
 460     }else{
 461       for(k=0;k<v->elements;k++)
 462         vN(new,j)[k]=ppt[k];
 463       v->max[firstentry]=firstmetric;
 464     }
 465   }
 466
 467   /* assign midpoints */
 468
 469   for(j=0;j<v->entries;j++){
 470 #ifdef NOISY
 471     fprintf(assig,"%ld\n",v->assigned[j]);
 472     fprintf(bias,"%g\n",v->bias[j]);
 473 #endif
 474     asserror+=fabs(v->assigned[j]-fdesired);
 475     if(v->assigned[j])
 476       for(k=0;k<v->elements;k++)
 477         _now(v,j)[k]=vN(new,j)[k]/v->assigned[j];
 478   }
 479
 480   asserror/=(v->entries*fdesired);
 481
 482   fprintf(stderr,"Pass #%d... ",v->it);
 483   fprintf(stderr,": dist %g(%g) metric error=%g \n",
 484           asserror,fdesired,meterror/v->points);
 485   v->it++;
 486
 487   free(new);
 488   free(nearcount);
 489   free(nearbias);
 490 #ifdef NOISY
 491   fclose(assig);
 492   fclose(bias);
 493   fclose(cells);
 494 #endif
 495   return(asserror);
 496 }
 497