vq/vqgen.c

   1 /********************************************************************
   2  *                                                                  *
   3  * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE.   *
   4  * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS     *
   5  * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE *
   6  * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING.       *
   7  *                                                                  *
   8  * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2001             *
   9  * by the XIPHOPHORUS Company http://www.xiph.org/                  *
  10
  11  ********************************************************************
  12
  13  function: train a VQ codebook
  14  last mod: $Id: vqgen.c,v 1.39 2001/02/26 03:51:13 xiphmont Exp $
  15
  16  ********************************************************************/
  17
  18 /* This code is *not* part of libvorbis.  It is used to generate
  19    trained codebooks offline and then spit the results into a
  20    pregenerated codebook that is compiled into libvorbis.  It is an
  21    expensive (but good) algorithm.  Run it on big iron. */
  22
  23 /* There are so many optimizations to explore in *both* stages that
  24    considering the undertaking is almost withering.  For now, we brute
  25    force it all */
  26
  27 #include <stdlib.h>
  28 #include <stdio.h>
  29 #include <math.h>
  30 #include <string.h>
  31
  32 #include "vqgen.h"
  33 #include "bookutil.h"
  34
  35 /* Codebook generation happens in two steps:
  36
  37    1) Train the codebook with data collected from the encoder: We use
  38    one of a few error metrics (which represent the distance between a
  39    given data point and a candidate point in the training set) to
  40    divide the training set up into cells representing roughly equal
  41    probability of occurring.
  42
  43    2) Generate the codebook and auxiliary data from the trained data set
  44 */
  45
  46 /* Codebook training ****************************************************
  47  *
  48  * The basic idea here is that a VQ codebook is like an m-dimensional
  49  * foam with n bubbles.  The bubbles compete for space/volume and are
  50  * 'pressurized' [biased] according to some metric.  The basic alg
  51  * iterates through allowing the bubbles to compete for space until
  52  * they converge (if the damping is dome properly) on a steady-state
  53  * solution. Individual input points, collected from libvorbis, are
  54  * used to train the algorithm monte-carlo style.  */
  55
  56 /* internal helpers *****************************************************/
  57 #define vN(data,i) (data+v->elements*i)
  58
  59 /* default metric; squared 'distance' from desired value. */
  60 float _dist(vqgen *v,float *a, float *b){
  61   int i;
  62   int el=v->elements;
  63   float acc=0.f;
  64   for(i=0;i<el;i++){
  65     float val=(a[i]-b[i]);
  66     acc+=val*val;
  67   }
  68   return sqrt(acc);
  69 }
  70
  71 float *_weight_null(vqgen *v,float *a){
  72   return a;
  73 }
  74
  75 /* *must* be beefed up. */
  76 void _vqgen_seed(vqgen *v){
  77   long i;
  78   for(i=0;i<v->entries;i++)
  79     memcpy(_now(v,i),_point(v,i),sizeof(float)*v->elements);
  80   v->seeded=1;
  81 }
  82
  83 int directdsort(const void *a, const void *b){
  84   float av=*((float *)a);
  85   float bv=*((float *)b);
  86   if(av>bv)return(-1);
  87   return(1);
  88 }
  89
  90 void vqgen_cellmetric(vqgen *v){
  91   int j,k;
  92   float min=-1.f,max=-1.f,mean=0.f,acc=0.f;
  93   long dup=0,unused=0;
  94  #ifdef NOISY
  95   int i;
  96    char buff[80];
  97    float spacings[v->entries];
  98    int count=0;
  99    FILE *cells;
 100    sprintf(buff,"cellspace%d.m",v->it);
 101    cells=fopen(buff,"w");
 102 #endif
 103
 104   /* minimum, maximum, cell spacing */
 105   for(j=0;j<v->entries;j++){
 106     float localmin=-1.;
 107
 108     for(k=0;k<v->entries;k++){
 109       if(j!=k){
 110         float this=_dist(v,_now(v,j),_now(v,k));
 111         if(this>0){
 112           if(v->assigned[k] && (localmin==-1 || this<localmin))
 113             localmin=this;
 114         }else{
 115           if(k<j){
 116             dup++;
 117             break;
 118           }
 119         }
 120       }
 121     }
 122     if(k<v->entries)continue;
 123
 124     if(v->assigned[j]==0){
 125       unused++;
 126       continue;
 127     }
 128
 129     localmin=v->max[j]+localmin/2; /* this gives us rough diameter */
 130     if(min==-1 || localmin<min)min=localmin;
 131     if(max==-1 || localmin>max)max=localmin;
 132     mean+=localmin;
 133     acc++;
 134 #ifdef NOISY
 135     spacings[count++]=localmin;
 136 #endif
 137   }
 138
 139   fprintf(stderr,"cell diameter: %.03g::%.03g::%.03g (%ld unused/%ld dup)\n",
 140           min,mean/acc,max,unused,dup);
 141
 142 #ifdef NOISY
 143   qsort(spacings,count,sizeof(float),directdsort);
 144   for(i=0;i<count;i++)
 145     fprintf(cells,"%g\n",spacings[i]);
 146   fclose(cells);
 147 #endif
 148
 149 }
 150
 151 /* External calls *******************************************************/
 152
 153 /* We have two forms of quantization; in the first, each vector
 154    element in the codebook entry is orthogonal.  Residues would use this
 155    quantization for example.
 156
 157    In the second, we have a sequence of monotonically increasing
 158    values that we wish to quantize as deltas (to save space).  We
 159    still need to quantize so that absolute values are accurate. For
 160    example, LSP quantizes all absolute values, but the book encodes
 161    distance between values because each successive value is larger
 162    than the preceeding value.  Thus the desired quantibits apply to
 163    the encoded (delta) values, not abs positions. This requires minor
 164    additional encode-side trickery. */
 165
 166 void vqgen_quantize(vqgen *v,quant_meta *q){
 167
 168   float maxdel;
 169   float mindel;
 170
 171   float delta;
 172   float maxquant=((1<<q->quant)-1);
 173
 174   int j,k;
 175
 176   mindel=maxdel=_now(v,0)[0];
 177
 178   for(j=0;j<v->entries;j++){
 179     float last=0.f;
 180     for(k=0;k<v->elements;k++){
 181       if(mindel>_now(v,j)[k]-last)mindel=_now(v,j)[k]-last;
 182       if(maxdel<_now(v,j)[k]-last)maxdel=_now(v,j)[k]-last;
 183       if(q->sequencep)last=_now(v,j)[k];
 184     }
 185   }
 186
 187
 188   /* first find the basic delta amount from the maximum span to be
 189      encoded.  Loosen the delta slightly to allow for additional error
 190      during sequence quantization */
 191
 192   delta=(maxdel-mindel)/((1<<q->quant)-1.5f);
 193
 194   q->min=_float32_pack(mindel);
 195   q->delta=_float32_pack(delta);
 196
 197   mindel=_float32_unpack(q->min);
 198   delta=_float32_unpack(q->delta);
 199
 200   for(j=0;j<v->entries;j++){
 201     float last=0;
 202     for(k=0;k<v->elements;k++){
 203       float val=_now(v,j)[k];
 204       float now=rint((val-last-mindel)/delta);
 205
 206       _now(v,j)[k]=now;
 207       if(now<0){
 208         /* be paranoid; this should be impossible */
 209         fprintf(stderr,"fault; quantized value<0\n");
 210         exit(1);
 211       }
 212
 213       if(now>maxquant){
 214         /* be paranoid; this should be impossible */
 215         fprintf(stderr,"fault; quantized value>max\n");
 216         exit(1);
 217       }
 218       if(q->sequencep)last=(now*delta)+mindel+last;
 219     }
 220   }
 221 }
 222
 223 /* much easier :-).  Unlike in the codebook, we don't un-log log
 224    scales; we just make sure they're properly offset. */
 225 void vqgen_unquantize(vqgen *v,quant_meta *q){
 226   long j,k;
 227   float mindel=_float32_unpack(q->min);
 228   float delta=_float32_unpack(q->delta);
 229
 230   for(j=0;j<v->entries;j++){
 231     float last=0.f;
 232     for(k=0;k<v->elements;k++){
 233       float now=_now(v,j)[k];
 234       now=fabs(now)*delta+last+mindel;
 235       if(q->sequencep)last=now;
 236       _now(v,j)[k]=now;
 237     }
 238   }
 239 }
 240
 241 void vqgen_init(vqgen *v,int elements,int aux,int entries,float mindist,
 242                 float  (*metric)(vqgen *,float *, float *),
 243                 float *(*weight)(vqgen *,float *),int centroid){
 244   memset(v,0,sizeof(vqgen));
 245
 246   v->centroid=centroid;
 247   v->elements=elements;
 248   v->aux=aux;
 249   v->mindist=mindist;
 250   v->allocated=32768;
 251   v->pointlist=_ogg_malloc(v->allocated*(v->elements+v->aux)*sizeof(float));
 252
 253   v->entries=entries;
 254   v->entrylist=_ogg_malloc(v->entries*v->elements*sizeof(float));
 255   v->assigned=_ogg_malloc(v->entries*sizeof(long));
 256   v->bias=_ogg_calloc(v->entries,sizeof(float));
 257   v->max=_ogg_calloc(v->entries,sizeof(float));
 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   v->asciipoints=tmpfile();
 268
 269 }
 270
 271 void vqgen_addpoint(vqgen *v, float *p,float *a){
 272   int k;
 273   for(k=0;k<v->elements;k++)
 274     fprintf(v->asciipoints,"%.12g\n",p[k]);
 275   for(k=0;k<v->aux;k++)
 276     fprintf(v->asciipoints,"%.12g\n",a[k]);
 277
 278   if(v->points>=v->allocated){
 279     v->allocated*=2;
 280     v->pointlist=_ogg_realloc(v->pointlist,v->allocated*(v->elements+v->aux)*
 281                          sizeof(float));
 282   }
 283
 284   memcpy(_point(v,v->points),p,sizeof(float)*v->elements);
 285   if(v->aux)memcpy(_point(v,v->points)+v->elements,a,sizeof(float)*v->aux);
 286
 287   /* quantize to the density mesh if it's selected */
 288   if(v->mindist>0.f){
 289     /* quantize to the mesh */
 290     for(k=0;k<v->elements+v->aux;k++)
 291       _point(v,v->points)[k]=
 292         rint(_point(v,v->points)[k]/v->mindist)*v->mindist;
 293   }
 294   v->points++;
 295   if(!(v->points&0xff))spinnit("loading... ",v->points);
 296 }
 297
 298 /* yes, not threadsafe.  These utils aren't */
 299 static int sortit=0;
 300 static int sortsize=0;
 301 static int meshcomp(const void *a,const void *b){
 302   if(((sortit++)&0xfff)==0)spinnit("sorting mesh...",sortit);
 303   return(memcmp(a,b,sortsize));
 304 }
 305
 306 void vqgen_sortmesh(vqgen *v){
 307   sortit=0;
 308   if(v->mindist>0.f){
 309     long i,march=1;
 310
 311     /* sort to make uniqueness detection trivial */
 312     sortsize=(v->elements+v->aux)*sizeof(float);
 313     qsort(v->pointlist,v->points,sortsize,meshcomp);
 314
 315     /* now march through and eliminate dupes */
 316     for(i=1;i<v->points;i++){
 317       if(memcmp(_point(v,i),_point(v,i-1),sortsize)){
 318         /* a new, unique entry.  march it down */
 319         if(i>march)memcpy(_point(v,march),_point(v,i),sortsize);
 320         march++;
 321       }
 322       spinnit("eliminating density... ",v->points-i);
 323     }
 324
 325     /* we're done */
 326     fprintf(stderr,"\r%ld training points remining out of %ld"
 327             " after density mesh (%ld%%)\n",march,v->points,march*100/v->points);
 328     v->points=march;
 329
 330   }
 331   v->sorted=1;
 332 }
 333
 334 float vqgen_iterate(vqgen *v,int biasp){
 335   long   i,j,k;
 336
 337   float fdesired;
 338   long  desired;
 339   long  desired2;
 340
 341   float asserror=0.f;
 342   float meterror=0.f;
 343   float *new;
 344   float *new2;
 345   long   *nearcount;
 346   float *nearbias;
 347  #ifdef NOISY
 348    char buff[80];
 349    FILE *assig;
 350    FILE *bias;
 351    FILE *cells;
 352    sprintf(buff,"cells%d.m",v->it);
 353    cells=fopen(buff,"w");
 354    sprintf(buff,"assig%d.m",v->it);
 355    assig=fopen(buff,"w");
 356    sprintf(buff,"bias%d.m",v->it);
 357    bias=fopen(buff,"w");
 358  #endif
 359
 360
 361   if(v->entries<2){
 362     fprintf(stderr,"generation requires at least two entries\n");
 363     exit(1);
 364   }
 365
 366   if(!v->sorted)vqgen_sortmesh(v);
 367   if(!v->seeded)_vqgen_seed(v);
 368
 369   fdesired=(float)v->points/v->entries;
 370   desired=fdesired;
 371   desired2=desired*2;
 372   new=_ogg_malloc(sizeof(float)*v->entries*v->elements);
 373   new2=_ogg_malloc(sizeof(float)*v->entries*v->elements);
 374   nearcount=_ogg_malloc(v->entries*sizeof(long));
 375   nearbias=_ogg_malloc(v->entries*desired2*sizeof(float));
 376
 377   /* fill in nearest points for entry biasing */
 378   /*memset(v->bias,0,sizeof(float)*v->entries);*/
 379   memset(nearcount,0,sizeof(long)*v->entries);
 380   memset(v->assigned,0,sizeof(long)*v->entries);
 381   if(biasp){
 382     for(i=0;i<v->points;i++){
 383       float *ppt=v->weight_func(v,_point(v,i));
 384       float firstmetric=v->metric_func(v,_now(v,0),ppt)+v->bias[0];
 385       float secondmetric=v->metric_func(v,_now(v,1),ppt)+v->bias[1];
 386       long   firstentry=0;
 387       long   secondentry=1;
 388
 389       if(!(i&0xff))spinnit("biasing... ",v->points+v->points+v->entries-i);
 390
 391       if(firstmetric>secondmetric){
 392         float temp=firstmetric;
 393         firstmetric=secondmetric;
 394         secondmetric=temp;
 395         firstentry=1;
 396         secondentry=0;
 397       }
 398
 399       for(j=2;j<v->entries;j++){
 400         float thismetric=v->metric_func(v,_now(v,j),ppt)+v->bias[j];
 401         if(thismetric<secondmetric){
 402           if(thismetric<firstmetric){
 403             secondmetric=firstmetric;
 404             secondentry=firstentry;
 405             firstmetric=thismetric;
 406             firstentry=j;
 407           }else{
 408             secondmetric=thismetric;
 409             secondentry=j;
 410           }
 411         }
 412       }
 413
 414       j=firstentry;
 415       for(j=0;j<v->entries;j++){
 416
 417         float thismetric,localmetric;
 418         float *nearbiasptr=nearbias+desired2*j;
 419         long k=nearcount[j];
 420
 421         localmetric=v->metric_func(v,_now(v,j),ppt);
 422         /* 'thismetric' is to be the bias value necessary in the current
 423            arrangement for entry j to capture point i */
 424         if(firstentry==j){
 425           /* use the secondary entry as the threshhold */
 426           thismetric=secondmetric-localmetric;
 427         }else{
 428           /* use the primary entry as the threshhold */
 429           thismetric=firstmetric-localmetric;
 430         }
 431
 432         /* support the idea of 'minimum distance'... if we want the
 433            cells in a codebook to be roughly some minimum size (as with
 434            the low resolution residue books) */
 435
 436         /* a cute two-stage delayed sorting hack */
 437         if(k<desired){
 438           nearbiasptr[k]=thismetric;
 439           k++;
 440           if(k==desired){
 441             spinnit("biasing... ",v->points+v->points+v->entries-i);
 442             qsort(nearbiasptr,desired,sizeof(float),directdsort);
 443           }
 444
 445         }else if(thismetric>nearbiasptr[desired-1]){
 446           nearbiasptr[k]=thismetric;
 447           k++;
 448           if(k==desired2){
 449             spinnit("biasing... ",v->points+v->points+v->entries-i);
 450             qsort(nearbiasptr,desired2,sizeof(float),directdsort);
 451             k=desired;
 452           }
 453         }
 454         nearcount[j]=k;
 455       }
 456     }
 457
 458     /* inflate/deflate */
 459
 460     for(i=0;i<v->entries;i++){
 461       float *nearbiasptr=nearbias+desired2*i;
 462
 463       spinnit("biasing... ",v->points+v->entries-i);
 464
 465       /* due to the delayed sorting, we likely need to finish it off....*/
 466       if(nearcount[i]>desired)
 467         qsort(nearbiasptr,nearcount[i],sizeof(float),directdsort);
 468
 469       v->bias[i]=nearbiasptr[desired-1];
 470
 471     }
 472   }else{
 473     memset(v->bias,0,v->entries*sizeof(float));
 474   }
 475
 476   /* Now assign with new bias and find new midpoints */
 477   for(i=0;i<v->points;i++){
 478     float *ppt=v->weight_func(v,_point(v,i));
 479     float firstmetric=v->metric_func(v,_now(v,0),ppt)+v->bias[0];
 480     long   firstentry=0;
 481
 482     if(!(i&0xff))spinnit("centering... ",v->points-i);
 483
 484     for(j=0;j<v->entries;j++){
 485       float thismetric=v->metric_func(v,_now(v,j),ppt)+v->bias[j];
 486       if(thismetric<firstmetric){
 487         firstmetric=thismetric;
 488         firstentry=j;
 489       }
 490     }
 491
 492     j=firstentry;
 493
 494 #ifdef NOISY
 495     fprintf(cells,"%g %g\n%g %g\n\n",
 496           _now(v,j)[0],_now(v,j)[1],
 497           ppt[0],ppt[1]);
 498 #endif
 499
 500     firstmetric-=v->bias[j];
 501     meterror+=firstmetric;
 502
 503     if(v->centroid==0){
 504       /* set up midpoints for next iter */
 505       if(v->assigned[j]++){
 506         for(k=0;k<v->elements;k++)
 507           vN(new,j)[k]+=ppt[k];
 508         if(firstmetric>v->max[j])v->max[j]=firstmetric;
 509       }else{
 510         for(k=0;k<v->elements;k++)
 511           vN(new,j)[k]=ppt[k];
 512         v->max[j]=firstmetric;
 513       }
 514     }else{
 515       /* centroid */
 516       if(v->assigned[j]++){
 517         for(k=0;k<v->elements;k++){
 518           if(vN(new,j)[k]>ppt[k])vN(new,j)[k]=ppt[k];
 519           if(vN(new2,j)[k]<ppt[k])vN(new2,j)[k]=ppt[k];
 520         }
 521         if(firstmetric>v->max[firstentry])v->max[j]=firstmetric;
 522       }else{
 523         for(k=0;k<v->elements;k++){
 524           vN(new,j)[k]=ppt[k];
 525           vN(new2,j)[k]=ppt[k];
 526         }
 527         v->max[firstentry]=firstmetric;
 528       }
 529     }
 530   }
 531
 532   /* assign midpoints */
 533
 534   for(j=0;j<v->entries;j++){
 535 #ifdef NOISY
 536     fprintf(assig,"%ld\n",v->assigned[j]);
 537     fprintf(bias,"%g\n",v->bias[j]);
 538 #endif
 539     asserror+=fabs(v->assigned[j]-fdesired);
 540     if(v->assigned[j]){
 541       if(v->centroid==0){
 542         for(k=0;k<v->elements;k++)
 543           _now(v,j)[k]=vN(new,j)[k]/v->assigned[j];
 544       }else{
 545         for(k=0;k<v->elements;k++)
 546           _now(v,j)[k]=(vN(new,j)[k]+vN(new2,j)[k])/2.f;
 547       }
 548     }
 549   }
 550
 551   asserror/=(v->entries*fdesired);
 552
 553   fprintf(stderr,"Pass #%d... ",v->it);
 554   fprintf(stderr,": dist %g(%g) metric error=%g \n",
 555           asserror,fdesired,meterror/v->points);
 556   v->it++;
 557
 558   free(new);
 559   free(nearcount);
 560   free(nearbias);
 561 #ifdef NOISY
 562   fclose(assig);
 563   fclose(bias);
 564   fclose(cells);
 565 #endif
 566   return(asserror);
 567 }
 568