vq/vqgen.c

   1 /********************************************************************
   2  *                                                                  *
   3  * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE.   *
   4  * USE, DISTRIBUTION AND REPRODUCTION OF THIS SOURCE IS GOVERNED BY *
   5  * THE GNU LESSER/LIBRARY PUBLIC LICENSE, WHICH IS INCLUDED WITH    *
   6  * THIS SOURCE. PLEASE READ THESE TERMS BEFORE DISTRIBUTING.        *
   7  *                                                                  *
   8  * THE OggVorbis 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.37 2000/12/21 21:04:49 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 float _dist(vqgen *v,float *a, float *b){
  62   int i;
  63   int el=v->elements;
  64   float acc=0.f;
  65   for(i=0;i<el;i++){
  66     float val=(a[i]-b[i]);
  67     acc+=val*val;
  68   }
  69   return sqrt(acc);
  70 }
  71
  72 float *_weight_null(vqgen *v,float *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(float)*v->elements);
  81   v->seeded=1;
  82 }
  83
  84 int directdsort(const void *a, const void *b){
  85   float av=*((float *)a);
  86   float bv=*((float *)b);
  87   if(av>bv)return(-1);
  88   return(1);
  89 }
  90
  91 void vqgen_cellmetric(vqgen *v){
  92   int j,k;
  93   float min=-1.f,max=-1.f,mean=0.f,acc=0.f;
  94   long dup=0,unused=0;
  95  #ifdef NOISY
  96   int i;
  97    char buff[80];
  98    float 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     float localmin=-1.;
 108
 109     for(k=0;k<v->entries;k++){
 110       if(j!=k){
 111         float 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(float),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   float maxdel;
 170   float mindel;
 171
 172   float delta;
 173   float 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     float last=0.f;
 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.5f);
 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     float last=0;
 203     for(k=0;k<v->elements;k++){
 204       float val=_now(v,j)[k];
 205       float 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   float mindel=_float32_unpack(q->min);
 229   float delta=_float32_unpack(q->delta);
 230
 231   for(j=0;j<v->entries;j++){
 232     float last=0.f;
 233     for(k=0;k<v->elements;k++){
 234       float 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,float mindist,
 243                 float  (*metric)(vqgen *,float *, float *),
 244                 float *(*weight)(vqgen *,float *),int centroid){
 245   memset(v,0,sizeof(vqgen));
 246
 247   v->centroid=centroid;
 248   v->elements=elements;
 249   v->aux=aux;
 250   v->mindist=mindist;
 251   v->allocated=32768;
 252   v->pointlist=_ogg_malloc(v->allocated*(v->elements+v->aux)*sizeof(float));
 253
 254   v->entries=entries;
 255   v->entrylist=_ogg_malloc(v->entries*v->elements*sizeof(float));
 256   v->assigned=_ogg_malloc(v->entries*sizeof(long));
 257   v->bias=_ogg_calloc(v->entries,sizeof(float));
 258   v->max=_ogg_calloc(v->entries,sizeof(float));
 259   if(metric)
 260     v->metric_func=metric;
 261   else
 262     v->metric_func=_dist;
 263   if(weight)
 264     v->weight_func=weight;
 265   else
 266     v->weight_func=_weight_null;
 267
 268   v->asciipoints=tmpfile();
 269
 270 }
 271
 272 void vqgen_addpoint(vqgen *v, float *p,float *a){
 273   int k;
 274   for(k=0;k<v->elements;k++)
 275     fprintf(v->asciipoints,"%.12g\n",p[k]);
 276   for(k=0;k<v->aux;k++)
 277     fprintf(v->asciipoints,"%.12g\n",a[k]);
 278
 279   if(v->points>=v->allocated){
 280     v->allocated*=2;
 281     v->pointlist=_ogg_realloc(v->pointlist,v->allocated*(v->elements+v->aux)*
 282                          sizeof(float));
 283   }
 284
 285   memcpy(_point(v,v->points),p,sizeof(float)*v->elements);
 286   if(v->aux)memcpy(_point(v,v->points)+v->elements,a,sizeof(float)*v->aux);
 287
 288   /* quantize to the density mesh if it's selected */
 289   if(v->mindist>0.f){
 290     /* quantize to the mesh */
 291     for(k=0;k<v->elements+v->aux;k++)
 292       _point(v,v->points)[k]=
 293         rint(_point(v,v->points)[k]/v->mindist)*v->mindist;
 294   }
 295   v->points++;
 296   if(!(v->points&0xff))spinnit("loading... ",v->points);
 297 }
 298
 299 /* yes, not threadsafe.  These utils aren't */
 300 static int sortit=0;
 301 static int sortsize=0;
 302 static int meshcomp(const void *a,const void *b){
 303   if(((sortit++)&0xfff)==0)spinnit("sorting mesh...",sortit);
 304   return(memcmp(a,b,sortsize));
 305 }
 306
 307 void vqgen_sortmesh(vqgen *v){
 308   sortit=0;
 309   if(v->mindist>0.f){
 310     long i,march=1;
 311
 312     /* sort to make uniqueness detection trivial */
 313     sortsize=(v->elements+v->aux)*sizeof(float);
 314     qsort(v->pointlist,v->points,sortsize,meshcomp);
 315
 316     /* now march through and eliminate dupes */
 317     for(i=1;i<v->points;i++){
 318       if(memcmp(_point(v,i),_point(v,i-1),sortsize)){
 319         /* a new, unique entry.  march it down */
 320         if(i>march)memcpy(_point(v,march),_point(v,i),sortsize);
 321         march++;
 322       }
 323       spinnit("eliminating density... ",v->points-i);
 324     }
 325
 326     /* we're done */
 327     fprintf(stderr,"\r%ld training points remining out of %ld"
 328             " after density mesh (%ld%%)\n",march,v->points,march*100/v->points);
 329     v->points=march;
 330
 331   }
 332   v->sorted=1;
 333 }
 334
 335 float vqgen_iterate(vqgen *v,int biasp){
 336   long   i,j,k;
 337
 338   float fdesired;
 339   long  desired;
 340   long  desired2;
 341
 342   float asserror=0.f;
 343   float meterror=0.f;
 344   float *new;
 345   float *new2;
 346   long   *nearcount;
 347   float *nearbias;
 348  #ifdef NOISY
 349    char buff[80];
 350    FILE *assig;
 351    FILE *bias;
 352    FILE *cells;
 353    sprintf(buff,"cells%d.m",v->it);
 354    cells=fopen(buff,"w");
 355    sprintf(buff,"assig%d.m",v->it);
 356    assig=fopen(buff,"w");
 357    sprintf(buff,"bias%d.m",v->it);
 358    bias=fopen(buff,"w");
 359  #endif
 360
 361
 362   if(v->entries<2){
 363     fprintf(stderr,"generation requires at least two entries\n");
 364     exit(1);
 365   }
 366
 367   if(!v->sorted)vqgen_sortmesh(v);
 368   if(!v->seeded)_vqgen_seed(v);
 369
 370   fdesired=(float)v->points/v->entries;
 371   desired=fdesired;
 372   desired2=desired*2;
 373   new=_ogg_malloc(sizeof(float)*v->entries*v->elements);
 374   new2=_ogg_malloc(sizeof(float)*v->entries*v->elements);
 375   nearcount=_ogg_malloc(v->entries*sizeof(long));
 376   nearbias=_ogg_malloc(v->entries*desired2*sizeof(float));
 377
 378   /* fill in nearest points for entry biasing */
 379   /*memset(v->bias,0,sizeof(float)*v->entries);*/
 380   memset(nearcount,0,sizeof(long)*v->entries);
 381   memset(v->assigned,0,sizeof(long)*v->entries);
 382   if(biasp){
 383     for(i=0;i<v->points;i++){
 384       float *ppt=v->weight_func(v,_point(v,i));
 385       float firstmetric=v->metric_func(v,_now(v,0),ppt)+v->bias[0];
 386       float secondmetric=v->metric_func(v,_now(v,1),ppt)+v->bias[1];
 387       long   firstentry=0;
 388       long   secondentry=1;
 389
 390       if(!(i&0xff))spinnit("biasing... ",v->points+v->points+v->entries-i);
 391
 392       if(firstmetric>secondmetric){
 393         float temp=firstmetric;
 394         firstmetric=secondmetric;
 395         secondmetric=temp;
 396         firstentry=1;
 397         secondentry=0;
 398       }
 399
 400       for(j=2;j<v->entries;j++){
 401         float thismetric=v->metric_func(v,_now(v,j),ppt)+v->bias[j];
 402         if(thismetric<secondmetric){
 403           if(thismetric<firstmetric){
 404             secondmetric=firstmetric;
 405             secondentry=firstentry;
 406             firstmetric=thismetric;
 407             firstentry=j;
 408           }else{
 409             secondmetric=thismetric;
 410             secondentry=j;
 411           }
 412         }
 413       }
 414
 415       j=firstentry;
 416       for(j=0;j<v->entries;j++){
 417
 418         float thismetric,localmetric;
 419         float *nearbiasptr=nearbias+desired2*j;
 420         long k=nearcount[j];
 421
 422         localmetric=v->metric_func(v,_now(v,j),ppt);
 423         /* 'thismetric' is to be the bias value necessary in the current
 424            arrangement for entry j to capture point i */
 425         if(firstentry==j){
 426           /* use the secondary entry as the threshhold */
 427           thismetric=secondmetric-localmetric;
 428         }else{
 429           /* use the primary entry as the threshhold */
 430           thismetric=firstmetric-localmetric;
 431         }
 432
 433         /* support the idea of 'minimum distance'... if we want the
 434            cells in a codebook to be roughly some minimum size (as with
 435            the low resolution residue books) */
 436
 437         /* a cute two-stage delayed sorting hack */
 438         if(k<desired){
 439           nearbiasptr[k]=thismetric;
 440           k++;
 441           if(k==desired){
 442             spinnit("biasing... ",v->points+v->points+v->entries-i);
 443             qsort(nearbiasptr,desired,sizeof(float),directdsort);
 444           }
 445
 446         }else if(thismetric>nearbiasptr[desired-1]){
 447           nearbiasptr[k]=thismetric;
 448           k++;
 449           if(k==desired2){
 450             spinnit("biasing... ",v->points+v->points+v->entries-i);
 451             qsort(nearbiasptr,desired2,sizeof(float),directdsort);
 452             k=desired;
 453           }
 454         }
 455         nearcount[j]=k;
 456       }
 457     }
 458
 459     /* inflate/deflate */
 460
 461     for(i=0;i<v->entries;i++){
 462       float *nearbiasptr=nearbias+desired2*i;
 463
 464       spinnit("biasing... ",v->points+v->entries-i);
 465
 466       /* due to the delayed sorting, we likely need to finish it off....*/
 467       if(nearcount[i]>desired)
 468         qsort(nearbiasptr,nearcount[i],sizeof(float),directdsort);
 469
 470       v->bias[i]=nearbiasptr[desired-1];
 471
 472     }
 473   }else{
 474     memset(v->bias,0,v->entries*sizeof(float));
 475   }
 476
 477   /* Now assign with new bias and find new midpoints */
 478   for(i=0;i<v->points;i++){
 479     float *ppt=v->weight_func(v,_point(v,i));
 480     float firstmetric=v->metric_func(v,_now(v,0),ppt)+v->bias[0];
 481     long   firstentry=0;
 482
 483     if(!(i&0xff))spinnit("centering... ",v->points-i);
 484
 485     for(j=0;j<v->entries;j++){
 486       float thismetric=v->metric_func(v,_now(v,j),ppt)+v->bias[j];
 487       if(thismetric<firstmetric){
 488         firstmetric=thismetric;
 489         firstentry=j;
 490       }
 491     }
 492
 493     j=firstentry;
 494
 495 #ifdef NOISY
 496     fprintf(cells,"%g %g\n%g %g\n\n",
 497           _now(v,j)[0],_now(v,j)[1],
 498           ppt[0],ppt[1]);
 499 #endif
 500
 501     firstmetric-=v->bias[j];
 502     meterror+=firstmetric;
 503
 504     if(v->centroid==0){
 505       /* set up midpoints for next iter */
 506       if(v->assigned[j]++){
 507         for(k=0;k<v->elements;k++)
 508           vN(new,j)[k]+=ppt[k];
 509         if(firstmetric>v->max[j])v->max[j]=firstmetric;
 510       }else{
 511         for(k=0;k<v->elements;k++)
 512           vN(new,j)[k]=ppt[k];
 513         v->max[j]=firstmetric;
 514       }
 515     }else{
 516       /* centroid */
 517       if(v->assigned[j]++){
 518         for(k=0;k<v->elements;k++){
 519           if(vN(new,j)[k]>ppt[k])vN(new,j)[k]=ppt[k];
 520           if(vN(new2,j)[k]<ppt[k])vN(new2,j)[k]=ppt[k];
 521         }
 522         if(firstmetric>v->max[firstentry])v->max[j]=firstmetric;
 523       }else{
 524         for(k=0;k<v->elements;k++){
 525           vN(new,j)[k]=ppt[k];
 526           vN(new2,j)[k]=ppt[k];
 527         }
 528         v->max[firstentry]=firstmetric;
 529       }
 530     }
 531   }
 532
 533   /* assign midpoints */
 534
 535   for(j=0;j<v->entries;j++){
 536 #ifdef NOISY
 537     fprintf(assig,"%ld\n",v->assigned[j]);
 538     fprintf(bias,"%g\n",v->bias[j]);
 539 #endif
 540     asserror+=fabs(v->assigned[j]-fdesired);
 541     if(v->assigned[j]){
 542       if(v->centroid==0){
 543         for(k=0;k<v->elements;k++)
 544           _now(v,j)[k]=vN(new,j)[k]/v->assigned[j];
 545       }else{
 546         for(k=0;k<v->elements;k++)
 547           _now(v,j)[k]=(vN(new,j)[k]+vN(new2,j)[k])/2.f;
 548       }
 549     }
 550   }
 551
 552   asserror/=(v->entries*fdesired);
 553
 554   fprintf(stderr,"Pass #%d... ",v->it);
 555   fprintf(stderr,": dist %g(%g) metric error=%g \n",
 556           asserror,fdesired,meterror/v->points);
 557   v->it++;
 558
 559   free(new);
 560   free(nearcount);
 561   free(nearbias);
 562 #ifdef NOISY
 563   fclose(assig);
 564   fclose(bias);
 565   fclose(cells);
 566 #endif
 567   return(asserror);
 568 }
 569