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