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.34 2000/10/12 03:13:02 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 float _dist(vqgen *v,float *a, float *b){
  63   int i;
  64   int el=v->elements;
  65   float acc=0.;
  66   for(i=0;i<el;i++){
  67     float val=(a[i]-b[i]);
  68     acc+=val*val;
  69   }
  70   return sqrt(acc);
  71 }
  72
  73 float *_weight_null(vqgen *v,float *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(float)*v->elements);
  82   v->seeded=1;
  83 }
  84
  85 int directdsort(const void *a, const void *b){
  86   float av=*((float *)a);
  87   float bv=*((float *)b);
  88   if(av>bv)return(-1);
  89   return(1);
  90 }
  91
  92 void vqgen_cellmetric(vqgen *v){
  93   int j,k;
  94   float min=-1.,max=-1.,mean=0.,acc=0.;
  95   long dup=0,unused=0;
  96  #ifdef NOISY
  97   int i;
  98    char buff[80];
  99    float spacings[v->entries];
 100    int count=0;
 101    FILE *cells;
 102    sprintf(buff,"cellspace%d.m",v->it);
 103    cells=fopen(buff,"w");
 104 #endif
 105
 106   /* minimum, maximum, cell spacing */
 107   for(j=0;j<v->entries;j++){
 108     float localmin=-1.;
 109
 110     for(k=0;k<v->entries;k++){
 111       if(j!=k){
 112         float this=_dist(v,_now(v,j),_now(v,k));
 113         if(this>0){
 114           if(v->assigned[k] && (localmin==-1 || this<localmin))
 115             localmin=this;
 116         }else{
 117           if(k<j){
 118             dup++;
 119             break;
 120           }
 121         }
 122       }
 123     }
 124     if(k<v->entries)continue;
 125
 126     if(v->assigned[j]==0){
 127       unused++;
 128       continue;
 129     }
 130
 131     localmin=v->max[j]+localmin/2; /* this gives us rough diameter */
 132     if(min==-1 || localmin<min)min=localmin;
 133     if(max==-1 || localmin>max)max=localmin;
 134     mean+=localmin;
 135     acc++;
 136 #ifdef NOISY
 137     spacings[count++]=localmin;
 138 #endif
 139   }
 140
 141   fprintf(stderr,"cell diameter: %.03g::%.03g::%.03g (%ld unused/%ld dup)\n",
 142           min,mean/acc,max,unused,dup);
 143
 144 #ifdef NOISY
 145   qsort(spacings,count,sizeof(float),directdsort);
 146   for(i=0;i<count;i++)
 147     fprintf(cells,"%g\n",spacings[i]);
 148   fclose(cells);
 149 #endif
 150
 151 }
 152
 153 /* External calls *******************************************************/
 154
 155 /* We have two forms of quantization; in the first, each vector
 156    element in the codebook entry is orthogonal.  Residues would use this
 157    quantization for example.
 158
 159    In the second, we have a sequence of monotonically increasing
 160    values that we wish to quantize as deltas (to save space).  We
 161    still need to quantize so that absolute values are accurate. For
 162    example, LSP quantizes all absolute values, but the book encodes
 163    distance between values because each successive value is larger
 164    than the preceeding value.  Thus the desired quantibits apply to
 165    the encoded (delta) values, not abs positions. This requires minor
 166    additional encode-side trickery. */
 167
 168 void vqgen_quantize(vqgen *v,quant_meta *q){
 169
 170   float maxdel;
 171   float mindel;
 172
 173   float delta;
 174   float maxquant=((1<<q->quant)-1);
 175
 176   int j,k;
 177
 178   mindel=maxdel=_now(v,0)[0];
 179
 180   for(j=0;j<v->entries;j++){
 181     float last=0.;
 182     for(k=0;k<v->elements;k++){
 183       if(mindel>_now(v,j)[k]-last)mindel=_now(v,j)[k]-last;
 184       if(maxdel<_now(v,j)[k]-last)maxdel=_now(v,j)[k]-last;
 185       if(q->sequencep)last=_now(v,j)[k];
 186     }
 187   }
 188
 189
 190   /* first find the basic delta amount from the maximum span to be
 191      encoded.  Loosen the delta slightly to allow for additional error
 192      during sequence quantization */
 193
 194   delta=(maxdel-mindel)/((1<<q->quant)-1.5);
 195
 196   q->min=_float32_pack(mindel);
 197   q->delta=_float32_pack(delta);
 198
 199   mindel=_float32_unpack(q->min);
 200   delta=_float32_unpack(q->delta);
 201
 202   for(j=0;j<v->entries;j++){
 203     float last=0;
 204     for(k=0;k<v->elements;k++){
 205       float val=_now(v,j)[k];
 206       float now=rint((val-last-mindel)/delta);
 207
 208       _now(v,j)[k]=now;
 209       if(now<0){
 210         /* be paranoid; this should be impossible */
 211         fprintf(stderr,"fault; quantized value<0\n");
 212         exit(1);
 213       }
 214
 215       if(now>maxquant){
 216         /* be paranoid; this should be impossible */
 217         fprintf(stderr,"fault; quantized value>max\n");
 218         exit(1);
 219       }
 220       if(q->sequencep)last=(now*delta)+mindel+last;
 221     }
 222   }
 223 }
 224
 225 /* much easier :-).  Unlike in the codebook, we don't un-log log
 226    scales; we just make sure they're properly offset. */
 227 void vqgen_unquantize(vqgen *v,quant_meta *q){
 228   long j,k;
 229   float mindel=_float32_unpack(q->min);
 230   float delta=_float32_unpack(q->delta);
 231
 232   for(j=0;j<v->entries;j++){
 233     float last=0.;
 234     for(k=0;k<v->elements;k++){
 235       float now=_now(v,j)[k];
 236       now=fabs(now)*delta+last+mindel;
 237       if(q->sequencep)last=now;
 238       _now(v,j)[k]=now;
 239     }
 240   }
 241 }
 242
 243 void vqgen_init(vqgen *v,int elements,int aux,int entries,float mindist,
 244                 float  (*metric)(vqgen *,float *, float *),
 245                 float *(*weight)(vqgen *,float *),int centroid){
 246   memset(v,0,sizeof(vqgen));
 247
 248   v->centroid=centroid;
 249   v->elements=elements;
 250   v->aux=aux;
 251   v->mindist=mindist;
 252   v->allocated=32768;
 253   v->pointlist=malloc(v->allocated*(v->elements+v->aux)*sizeof(float));
 254
 255   v->entries=entries;
 256   v->entrylist=malloc(v->entries*v->elements*sizeof(float));
 257   v->assigned=malloc(v->entries*sizeof(long));
 258   v->bias=calloc(v->entries,sizeof(float));
 259   v->max=calloc(v->entries,sizeof(float));
 260   if(metric)
 261     v->metric_func=metric;
 262   else
 263     v->metric_func=_dist;
 264   if(weight)
 265     v->weight_func=weight;
 266   else
 267     v->weight_func=_weight_null;
 268
 269   v->asciipoints=tmpfile();
 270
 271 }
 272
 273 void vqgen_addpoint(vqgen *v, float *p,float *a){
 274   int k;
 275   for(k=0;k<v->elements;k++)
 276     fprintf(v->asciipoints,"%.12g\n",p[k]);
 277   for(k=0;k<v->aux;k++)
 278     fprintf(v->asciipoints,"%.12g\n",a[k]);
 279
 280   if(v->points>=v->allocated){
 281     v->allocated*=2;
 282     v->pointlist=realloc(v->pointlist,v->allocated*(v->elements+v->aux)*
 283                          sizeof(float));
 284   }
 285
 286   memcpy(_point(v,v->points),p,sizeof(float)*v->elements);
 287   if(v->aux)memcpy(_point(v,v->points)+v->elements,a,sizeof(float)*v->aux);
 288
 289   /* quantize to the density mesh if it's selected */
 290   if(v->mindist>0.){
 291     /* quantize to the mesh */
 292     for(k=0;k<v->elements+v->aux;k++)
 293       _point(v,v->points)[k]=
 294         rint(_point(v,v->points)[k]/v->mindist)*v->mindist;
 295   }
 296   v->points++;
 297   if(!(v->points&0xff))spinnit("loading... ",v->points);
 298 }
 299
 300 /* yes, not threadsafe.  These utils aren't */
 301 static int sortit=0;
 302 static int sortsize=0;
 303 static int meshcomp(const void *a,const void *b){
 304   if(((sortit++)&0xfff)==0)spinnit("sorting mesh...",sortit);
 305   return(memcmp(a,b,sortsize));
 306 }
 307
 308 void vqgen_sortmesh(vqgen *v){
 309   sortit=0;
 310   if(v->mindist>0.){
 311     long i,march=1;
 312
 313     /* sort to make uniqueness detection trivial */
 314     sortsize=(v->elements+v->aux)*sizeof(float);
 315     qsort(v->pointlist,v->points,sortsize,meshcomp);
 316
 317     /* now march through and eliminate dupes */
 318     for(i=1;i<v->points;i++){
 319       if(memcmp(_point(v,i),_point(v,i-1),sortsize)){
 320         /* a new, unique entry.  march it down */
 321         if(i>march)memcpy(_point(v,march),_point(v,i),sortsize);
 322         march++;
 323       }
 324       spinnit("eliminating density... ",v->points-i);
 325     }
 326
 327     /* we're done */
 328     fprintf(stderr,"\r%ld training points remining out of %ld"
 329             " after density mesh (%ld%%)\n",march,v->points,march*100/v->points);
 330     v->points=march;
 331
 332   }
 333   v->sorted=1;
 334 }
 335
 336 float vqgen_iterate(vqgen *v,int biasp){
 337   long   i,j,k;
 338
 339   float fdesired;
 340   long  desired;
 341   long  desired2;
 342
 343   float asserror=0.;
 344   float meterror=0.;
 345   float *new;
 346   float *new2;
 347   long   *nearcount;
 348   float *nearbias;
 349  #ifdef NOISY
 350    char buff[80];
 351    FILE *assig;
 352    FILE *bias;
 353    FILE *cells;
 354    sprintf(buff,"cells%d.m",v->it);
 355    cells=fopen(buff,"w");
 356    sprintf(buff,"assig%d.m",v->it);
 357    assig=fopen(buff,"w");
 358    sprintf(buff,"bias%d.m",v->it);
 359    bias=fopen(buff,"w");
 360  #endif
 361
 362
 363   if(v->entries<2){
 364     fprintf(stderr,"generation requires at least two entries\n");
 365     exit(1);
 366   }
 367
 368   if(!v->sorted)vqgen_sortmesh(v);
 369   if(!v->seeded)_vqgen_seed(v);
 370
 371   fdesired=(float)v->points/v->entries;
 372   desired=fdesired;
 373   desired2=desired*2;
 374   new=malloc(sizeof(float)*v->entries*v->elements);
 375   new2=malloc(sizeof(float)*v->entries*v->elements);
 376   nearcount=malloc(v->entries*sizeof(long));
 377   nearbias=malloc(v->entries*desired2*sizeof(float));
 378
 379   /* fill in nearest points for entry biasing */
 380   /*memset(v->bias,0,sizeof(float)*v->entries);*/
 381   memset(nearcount,0,sizeof(long)*v->entries);
 382   memset(v->assigned,0,sizeof(long)*v->entries);
 383   if(biasp){
 384     for(i=0;i<v->points;i++){
 385       float *ppt=v->weight_func(v,_point(v,i));
 386       float firstmetric=v->metric_func(v,_now(v,0),ppt)+v->bias[0];
 387       float secondmetric=v->metric_func(v,_now(v,1),ppt)+v->bias[1];
 388       long   firstentry=0;
 389       long   secondentry=1;
 390
 391       if(!(i&0xff))spinnit("biasing... ",v->points+v->points+v->entries-i);
 392
 393       if(firstmetric>secondmetric){
 394         float temp=firstmetric;
 395         firstmetric=secondmetric;
 396         secondmetric=temp;
 397         firstentry=1;
 398         secondentry=0;
 399       }
 400
 401       for(j=2;j<v->entries;j++){
 402         float thismetric=v->metric_func(v,_now(v,j),ppt)+v->bias[j];
 403         if(thismetric<secondmetric){
 404           if(thismetric<firstmetric){
 405             secondmetric=firstmetric;
 406             secondentry=firstentry;
 407             firstmetric=thismetric;
 408             firstentry=j;
 409           }else{
 410             secondmetric=thismetric;
 411             secondentry=j;
 412           }
 413         }
 414       }
 415
 416       j=firstentry;
 417       for(j=0;j<v->entries;j++){
 418
 419         float thismetric,localmetric;
 420         float *nearbiasptr=nearbias+desired2*j;
 421         long k=nearcount[j];
 422
 423         localmetric=v->metric_func(v,_now(v,j),ppt);
 424         /* 'thismetric' is to be the bias value necessary in the current
 425            arrangement for entry j to capture point i */
 426         if(firstentry==j){
 427           /* use the secondary entry as the threshhold */
 428           thismetric=secondmetric-localmetric;
 429         }else{
 430           /* use the primary entry as the threshhold */
 431           thismetric=firstmetric-localmetric;
 432         }
 433
 434         /* support the idea of 'minimum distance'... if we want the
 435            cells in a codebook to be roughly some minimum size (as with
 436            the low resolution residue books) */
 437
 438         /* a cute two-stage delayed sorting hack */
 439         if(k<desired){
 440           nearbiasptr[k]=thismetric;
 441           k++;
 442           if(k==desired){
 443             spinnit("biasing... ",v->points+v->points+v->entries-i);
 444             qsort(nearbiasptr,desired,sizeof(float),directdsort);
 445           }
 446
 447         }else if(thismetric>nearbiasptr[desired-1]){
 448           nearbiasptr[k]=thismetric;
 449           k++;
 450           if(k==desired2){
 451             spinnit("biasing... ",v->points+v->points+v->entries-i);
 452             qsort(nearbiasptr,desired2,sizeof(float),directdsort);
 453             k=desired;
 454           }
 455         }
 456         nearcount[j]=k;
 457       }
 458     }
 459
 460     /* inflate/deflate */
 461
 462     for(i=0;i<v->entries;i++){
 463       float *nearbiasptr=nearbias+desired2*i;
 464
 465       spinnit("biasing... ",v->points+v->entries-i);
 466
 467       /* due to the delayed sorting, we likely need to finish it off....*/
 468       if(nearcount[i]>desired)
 469         qsort(nearbiasptr,nearcount[i],sizeof(float),directdsort);
 470
 471       v->bias[i]=nearbiasptr[desired-1];
 472
 473     }
 474   }else{
 475     memset(v->bias,0,v->entries*sizeof(float));
 476   }
 477
 478   /* Now assign with new bias and find new midpoints */
 479   for(i=0;i<v->points;i++){
 480     float *ppt=v->weight_func(v,_point(v,i));
 481     float firstmetric=v->metric_func(v,_now(v,0),ppt)+v->bias[0];
 482     long   firstentry=0;
 483
 484     if(!(i&0xff))spinnit("centering... ",v->points-i);
 485
 486     for(j=0;j<v->entries;j++){
 487       float thismetric=v->metric_func(v,_now(v,j),ppt)+v->bias[j];
 488       if(thismetric<firstmetric){
 489         firstmetric=thismetric;
 490         firstentry=j;
 491       }
 492     }
 493
 494     j=firstentry;
 495
 496 #ifdef NOISY
 497     fprintf(cells,"%g %g\n%g %g\n\n",
 498           _now(v,j)[0],_now(v,j)[1],
 499           ppt[0],ppt[1]);
 500 #endif
 501
 502     firstmetric-=v->bias[j];
 503     meterror+=firstmetric;
 504
 505     if(v->centroid==0){
 506       /* set up midpoints for next iter */
 507       if(v->assigned[j]++){
 508         for(k=0;k<v->elements;k++)
 509           vN(new,j)[k]+=ppt[k];
 510         if(firstmetric>v->max[j])v->max[j]=firstmetric;
 511       }else{
 512         for(k=0;k<v->elements;k++)
 513           vN(new,j)[k]=ppt[k];
 514         v->max[j]=firstmetric;
 515       }
 516     }else{
 517       /* centroid */
 518       if(v->assigned[j]++){
 519         for(k=0;k<v->elements;k++){
 520           if(vN(new,j)[k]>ppt[k])vN(new,j)[k]=ppt[k];
 521           if(vN(new2,j)[k]<ppt[k])vN(new2,j)[k]=ppt[k];
 522         }
 523         if(firstmetric>v->max[firstentry])v->max[j]=firstmetric;
 524       }else{
 525         for(k=0;k<v->elements;k++){
 526           vN(new,j)[k]=ppt[k];
 527           vN(new2,j)[k]=ppt[k];
 528         }
 529         v->max[firstentry]=firstmetric;
 530       }
 531     }
 532   }
 533
 534   /* assign midpoints */
 535
 536   for(j=0;j<v->entries;j++){
 537 #ifdef NOISY
 538     fprintf(assig,"%ld\n",v->assigned[j]);
 539     fprintf(bias,"%g\n",v->bias[j]);
 540 #endif
 541     asserror+=fabs(v->assigned[j]-fdesired);
 542     if(v->assigned[j]){
 543       if(v->centroid==0){
 544         for(k=0;k<v->elements;k++)
 545           _now(v,j)[k]=vN(new,j)[k]/v->assigned[j];
 546       }else{
 547         for(k=0;k<v->elements;k++)
 548           _now(v,j)[k]=(vN(new,j)[k]+vN(new2,j)[k])/2.;
 549       }
 550     }
 551   }
 552
 553   asserror/=(v->entries*fdesired);
 554
 555   fprintf(stderr,"Pass #%d... ",v->it);
 556   fprintf(stderr,": dist %g(%g) metric error=%g \n",
 557           asserror,fdesired,meterror/v->points);
 558   v->it++;
 559
 560   free(new);
 561   free(nearcount);
 562   free(nearbias);
 563 #ifdef NOISY
 564   fclose(assig);
 565   fclose(bias);
 566   fclose(cells);
 567 #endif
 568   return(asserror);
 569 }
 570