vq/vqsplit.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-1999             *
   9  * by 1999 Monty <monty@xiph.org> and The XIPHOPHORUS Company       *
  10  * http://www.xiph.org/                                             *
  11  *                                                                  *
  12  ********************************************************************
  13
  14  function: build a VQ codebook
  15  author: Monty <xiphmont@mit.edu>
  16  modifications by: Monty
  17  last modification date: Dec 10 1999
  18
  19  ********************************************************************/
  20
  21 /* This code is *not* part of libvorbis.  It is used to generate
  22    trained codebooks offline and then spit the results into a
  23    pregenerated codebook that is compiled into libvorbis.  It is an
  24    expensive (but good) algorithm.  Run it on big iron. */
  25
  26 /* There are so many optimizations to explore in *both* stages that
  27    considering the undertaking is almost withering.  For now, we brute
  28    force it all */
  29
  30 #include <stdlib.h>
  31 #include <stdio.h>
  32 #include <math.h>
  33 #include <string.h>
  34 #include "vqgen.h"
  35 #include <sys/time.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 /* Building a codebook from trained set **********************************
  49
  50    The codebook in raw form is technically finished once it's trained.
  51    However, we want to finalize the representative codebook values for
  52    each entry and generate auxiliary information to optimize encoding.
  53    We generate the auxiliary coding tree using collected data,
  54    probably the same data as in the original training */
  55
  56 /* At each recursion, the data set is split in half.  Cells with data
  57    points on side A go into set A, same with set B.  The sets may
  58    overlap.  If the cell overlaps the deviding line only very slightly
  59    (provided parameter), we may choose to ignore the overlap in order
  60    to pare the tree down */
  61
  62 double *sortvals;
  63 int els;
  64 int dascsort(const void *a,const void *b){
  65   double av=sortvals[*((long *)a) * els];
  66   double bv=sortvals[*((long *)b) * els];
  67   if(av<bv)return(-1);
  68   return(1);
  69 }
  70
  71 long *isortvals;
  72 int iascsort(const void *a,const void *b){
  73   long av=isortvals[*((long *)a)];
  74   long bv=isortvals[*((long *)b)];
  75   return(av-bv);
  76 }
  77
  78 extern double _dist_sq(vqgen *v,double *a, double *b);
  79
  80 /* goes through the split, but just counts it and returns a metric*/
  81 void vqsp_count(vqgen *v,long *entryindex,long entries,
  82                 long *pointindex,long points,
  83                 long *entryA,long *entryB,
  84                 double *n, double c, double slack,
  85                 long *entriesA,long *entriesB,long *entriesC){
  86   long i,j,k;
  87   long A=0,B=0,C=0;
  88
  89   memset(entryA,0,sizeof(long)*entries);
  90   memset(entryB,0,sizeof(long)*entries);
  91
  92   for(i=0;i<points;i++){
  93     double *ppt=_point(v,pointindex[i]);
  94     long   firstentry=0;
  95     double firstmetric=_dist_sq(v,_now(v,entryindex[0]),ppt);
  96     double position=-c;
  97
  98     for(j=1;j<entries;j++){
  99       double thismetric=_dist_sq(v,_now(v,entryindex[j]),ppt);
 100       if(thismetric<firstmetric){
 101         firstmetric=thismetric;
 102         firstentry=j;
 103       }
 104     }
 105
 106     /* count point split */
 107     for(k=0;k<v->elements;k++)
 108       position+=ppt[k]*n[k];
 109     if(position>0.){
 110       entryA[firstentry]++;
 111     }else{
 112       entryB[firstentry]++;
 113     }
 114   }
 115
 116   /* look to see if entries are in the slack zone */
 117   /* The entry splitting isn't total, so that storage has to be
 118      allocated for recursion.  Reuse the entryA/entryB vectors */
 119   for(j=0;j<entries;j++){
 120     long total=entryA[j]+entryB[j];
 121     if((double)entryA[j]/total<slack){
 122       entryA[j]=0;
 123     }else if((double)entryB[j]/total<slack){
 124       entryB[j]=0;
 125     }
 126     if(entryA[j] && entryB[j])C++;
 127     if(entryA[j])entryA[A++]=entryindex[j];
 128     if(entryB[j])entryB[B++]=entryindex[j];
 129   }
 130   *entriesA=A;
 131   *entriesB=B;
 132   *entriesC=C;
 133 }
 134
 135 void pq_in_out(vqgen *v,double *n,double *c,double *p,double *q){
 136   int k;
 137   *c=0.;
 138   for(k=0;k<v->elements;k++){
 139     double center=(p[k]+q[k])/2.;
 140     n[k]=(center-q[k])*2.;
 141     *c+=center*n[k];
 142   }
 143 }
 144
 145 void pq_center_out(vqgen *v,double *n,double *c,double *center,double *q){
 146   int k;
 147   *c=0.;
 148   for(k=0;k<v->elements;k++){
 149     n[k]=(center[k]-q[k])*2.;
 150     *c+=center[k]*n[k];
 151   }
 152 }
 153
 154 static void spinnit(void){
 155   static int p=0;
 156   static long lasttime=0;
 157   long test;
 158   struct timeval thistime;
 159
 160   gettimeofday(&thistime,NULL);
 161   test=thistime.tv_sec*10+thistime.tv_usec/100000;
 162   if(lasttime!=test){
 163     lasttime=test;
 164
 165     p++;if(p>3)p=0;
 166     switch(p){
 167     case 0:
 168       fprintf(stderr,"|\b");
 169       break;
 170     case 1:
 171       fprintf(stderr,"/\b");
 172       break;
 173     case 2:
 174       fprintf(stderr,"-\b");
 175       break;
 176     case 3:
 177       fprintf(stderr,"\\\b");
 178       break;
 179     }
 180     fflush(stderr);
 181   }
 182 }
 183
 184 int lp_split(vqgen *v,vqbook *b,
 185              long *entryindex,long entries,
 186              long *pointindex,long points,
 187              long depth,double slack){
 188
 189   /* The encoder, regardless of book, will be using a straight
 190      euclidian distance-to-point metric to determine closest point.
 191      Thus we split the cells using the same (we've already trained the
 192      codebook set spacing and distribution using special metrics and
 193      even a midpoint division won't disturb the basic properties) */
 194
 195   long ret;
 196   double *p;
 197   double *q;
 198   double *n;
 199   double c;
 200   long *entryA=calloc(entries,sizeof(long));
 201   long *entryB=calloc(entries,sizeof(long));
 202   long entriesA=0;
 203   long entriesB=0;
 204   long entriesC=0;
 205   long pointsA=0;
 206   long i,j,k;
 207
 208   p=alloca(sizeof(double)*v->elements);
 209   q=alloca(sizeof(double)*v->elements);
 210   n=alloca(sizeof(double)*v->elements);
 211   memset(p,0,sizeof(double)*v->elements);
 212
 213   /* We need to find the dividing hyperplane. find the median of each
 214      axis as the centerpoint and the normal facing farthest point */
 215
 216   /* more than one way to do this part.  For small sets, we can brute
 217      force it. */
 218
 219   if(entries<32){
 220     /* try every pair possibility */
 221     double best=0;
 222     long   besti=0;
 223     long   bestj=0;
 224     double this;
 225     for(i=0;i<entries-1;i++){
 226       for(j=i+1;j<entries;j++){
 227         spinnit();
 228         pq_in_out(v,n,&c,_now(v,entryindex[i]),_now(v,entryindex[j]));
 229         vqsp_count(v,entryindex,entries,
 230                  pointindex,points,
 231                  entryA,entryB,
 232                  n, c, slack,
 233                  &entriesA,&entriesB,&entriesC);
 234         this=(entriesA-entriesC)*(entriesB-entriesC);
 235
 236         if(this>best){
 237           best=this;
 238           besti=i;
 239           bestj=j;
 240         }
 241       }
 242     }
 243     pq_in_out(v,n,&c,_now(v,entryindex[besti]),_now(v,entryindex[bestj]));
 244   }else{
 245     double best=0.;
 246     long   bestj=0;
 247
 248     /* try COG/normal and furthest pairs */
 249     /* medianpoint */
 250     for(k=0;k<v->elements;k++){
 251       spinnit();
 252
 253       /* just sort the index array */
 254       sortvals=v->entrylist+k;
 255       els=v->elements;
 256       qsort(entryindex,entries,sizeof(long),dascsort);
 257       if(entries&0x1){
 258         p[k]=v->entrylist[(entryindex[entries/2])*v->elements+k];
 259       }else{
 260         p[k]=(v->entrylist[(entryindex[entries/2])*v->elements+k]+
 261               v->entrylist[(entryindex[entries/2-1])*v->elements+k])/2.;
 262       }
 263     }
 264
 265     spinnit();
 266
 267     /* try every normal, but just for distance */
 268     for(j=0;j<entries;j++){
 269       double *ppj=_now(v,entryindex[j]);
 270       double this=_dist_sq(v,p,ppj);
 271       if(this>best){
 272         best=this;
 273         bestj=j;
 274       }
 275     }
 276
 277     pq_center_out(v,n,&c,p,_point(v,pointindex[bestj]));
 278
 279
 280   }
 281
 282   /* find cells enclosing points */
 283   /* count A/B points */
 284
 285   vqsp_count(v,entryindex,entries,
 286            pointindex,points,
 287            entryA,entryB,
 288            n, c, slack,
 289            &entriesA,&entriesB,&entriesC);
 290
 291   /* the point index is split evenly, so we do an Order n
 292      rearrangement into A first/B last and just pass it on */
 293   {
 294     long Aptr=0;
 295     long Bptr=points-1;
 296     while(Aptr<=Bptr){
 297       while(Aptr<=Bptr){
 298         double position=-c;
 299         for(k=0;k<v->elements;k++)
 300           position+=_point(v,pointindex[Aptr])[k]*n[k];
 301         if(position<=0.)break; /* not in A */
 302         Aptr++;
 303       }
 304       while(Aptr<=Bptr){
 305         double position=-c;
 306         for(k=0;k<v->elements;k++)
 307           position+=_point(v,pointindex[Bptr])[k]*n[k];
 308         if(position>0.)break; /* not in B */
 309         Bptr--;
 310       }
 311       if(Aptr<Bptr){
 312         long temp=pointindex[Aptr];
 313         pointindex[Aptr]=pointindex[Bptr];
 314         pointindex[Bptr]=temp;
 315       }
 316       pointsA=Aptr;
 317     }
 318   }
 319
 320   fprintf(stderr,"split: total=%ld depth=%ld set A=%ld:%ld:%ld=B\n",
 321           entries,depth,entriesA-entriesC,entriesC,entriesB-entriesC);
 322   {
 323     long thisaux=b->aux++;
 324     if(b->aux>=b->alloc){
 325       b->alloc*=2;
 326       b->ptr0=realloc(b->ptr0,sizeof(long)*b->alloc);
 327       b->ptr1=realloc(b->ptr1,sizeof(long)*b->alloc);
 328       b->n=realloc(b->n,sizeof(double)*b->dim*b->alloc);
 329       b->c=realloc(b->c,sizeof(double)*b->alloc);
 330     }
 331
 332     memcpy(b->n+b->dim*thisaux,n,sizeof(double)*v->elements);
 333     b->c[thisaux]=c;
 334
 335     if(entriesA==1){
 336       ret=1;
 337       b->ptr0[thisaux]=entryA[0];
 338     }else{
 339       b->ptr0[thisaux]= -b->aux;
 340       ret=lp_split(v,b,entryA,entriesA,pointindex,pointsA,depth+1,slack);
 341     }
 342     if(entriesB==1){
 343       ret++;
 344       b->ptr1[thisaux]=entryB[0];
 345     }else{
 346       b->ptr1[thisaux]= -b->aux;
 347       ret+=lp_split(v,b,entryB,entriesB,pointindex+pointsA,points-pointsA,
 348                    depth+1,slack);
 349     }
 350   }
 351   free(entryA);
 352   free(entryB);
 353   return(ret);
 354 }
 355
 356 void vqsp_book(vqgen *v, vqbook *b){
 357   long *entryindex=malloc(sizeof(double)*v->entries);
 358   long *pointindex=malloc(sizeof(double)*v->points);
 359   long i,j;
 360
 361   memset(b,0,sizeof(vqbook));
 362   for(i=0;i<v->entries;i++)entryindex[i]=i;
 363   for(i=0;i<v->points;i++)pointindex[i]=i;
 364
 365   b->dim=v->elements;
 366   b->entries=v->entries;
 367   b->alloc=4096;
 368
 369   b->ptr0=malloc(sizeof(long)*b->alloc);
 370   b->ptr1=malloc(sizeof(long)*b->alloc);
 371   b->n=malloc(sizeof(double)*b->dim*b->alloc);
 372   b->c=malloc(sizeof(double)*b->alloc);
 373   b->lengthlist=calloc(b->entries,sizeof(long));
 374
 375   /* first, generate the encoding decision heirarchy */
 376   fprintf(stderr,"Total leaves: %d\n",
 377           lp_split(v,b,entryindex,v->entries, pointindex,v->points,0,0));
 378
 379   free(entryindex);
 380   free(pointindex);
 381
 382   /* run all training points through the decision tree to get a final
 383      probability count */
 384   {
 385     long *probability=calloc(b->entries,sizeof(long));
 386     long *membership=malloc(b->entries*sizeof(long));
 387
 388     for(i=0;i<v->points;i++){
 389       int ret=vqenc_entry(b,v->pointlist+i*v->elements);
 390       probability[ret]++;
 391     }
 392     for(i=0;i<v->entries;i++)membership[i]=i;
 393
 394     /* find codeword lengths */
 395     /* much more elegant means exist.  Brute force n^2, minimum thought */
 396     for(i=v->entries;i>1;i--){
 397       int first=-1,second=-1;
 398       long least=v->points+1;
 399
 400       /* find the two nodes to join */
 401       for(j=0;j<v->entries;j++)
 402         if(probability[j]<least){
 403           least=probability[j];
 404           first=membership[j];
 405         }
 406       least=v->points+1;
 407       for(j=0;j<v->entries;j++)
 408         if(probability[j]<least && membership[j]!=first){
 409           least=probability[j];
 410           second=membership[j];
 411         }
 412       if(first==-1 || second==-1){
 413         fprintf(stderr,"huffman fault; no free branch\n");
 414         exit(1);
 415       }
 416
 417       /* join them */
 418       least=probability[first]+probability[second];
 419       for(j=0;j<v->entries;j++)
 420         if(membership[j]==first || membership[j]==second){
 421           membership[j]=first;
 422           probability[j]=least;
 423           b->lengthlist[j]++;
 424         }
 425     }
 426     for(i=0;i<v->entries-1;i++)
 427       if(membership[i]!=membership[i+1]){
 428         fprintf(stderr,"huffman fault; failed to build single tree\n");
 429         exit(1);
 430       }
 431
 432     /* unneccessary metric */
 433     memset(probability,0,sizeof(long)*v->entries);
 434     for(i=0;i<v->points;i++){
 435       int ret=vqenc_entry(b,v->pointlist+i*v->elements);
 436       probability[ret]++;
 437     }
 438
 439     /* print the results */
 440     for(i=0;i<v->entries;i++)
 441       fprintf(stderr,"%ld: count=%ld codelength=%ld\n",i,probability[i],b->lengthlist[i]);
 442
 443     free(probability);
 444     free(membership);
 445   }
 446
 447   /* Sort the entries by codeword length, short to long (eases
 448      assignment and packing to do it now) */
 449   {
 450     long *wordlen=b->lengthlist;
 451     long *index=malloc(b->entries*sizeof(long));
 452     long *revindex=malloc(b->entries*sizeof(long));
 453     int k;
 454     for(i=0;i<b->entries;i++)index[i]=i;
 455     isortvals=b->lengthlist;
 456     qsort(index,b->entries,sizeof(long),iascsort);
 457
 458     /* rearrange storage; ptr0/1 first as it needs a reverse index */
 459     /* n and c stay unchanged */
 460     for(i=0;i<b->entries;i++)revindex[index[i]]=i;
 461     for(i=0;i<b->aux;i++){
 462       if(b->ptr0[i]>=0)b->ptr0[i]=revindex[i];
 463       if(b->ptr1[i]>=0)b->ptr1[i]=revindex[i];
 464     }
 465     free(revindex);
 466
 467     /* map lengthlist and vallist with index */
 468     b->lengthlist=calloc(b->entries,sizeof(long));
 469     b->valuelist=malloc(sizeof(double)*b->entries*b->dim);
 470     for(i=0;i<b->entries;i++){
 471       long e=index[i];
 472       for(k=0;k<b->dim;k++)
 473         b->valuelist[i*b->dim+k]=v->entrylist[e*b->dim+k];
 474       b->lengthlist[i]=wordlen[e];
 475     }
 476
 477     free(wordlen);
 478   }
 479
 480   /* generate the codewords (short to long) */
 481   {
 482     long current=0;
 483     long length=0;
 484     b->codelist=malloc(sizeof(long)*b->entries);
 485     for(i=0;i<b->entries;i++){
 486       if(length != b->lengthlist[i]){
 487         current<<=(b->lengthlist[i]-length);
 488         length=b->lengthlist[i];
 489       }
 490       b->codelist[i]=current;
 491       fprintf(stderr,"codeword %ld: %ld (length %d)\n",
 492               i, current, b->lengthlist[i]);
 493       current++;
 494     }
 495   }
 496
 497   /* sanity check the codewords */
 498   for(i=0;i<b->entries;i++){
 499     for(j=i+1;j<b->entries;j++){
 500       if(b->codelist[i]==b->codelist[j]){
 501         fprintf(stderr,"Error; codewords for %ld and %ld both equal %lx\n",
 502                 i, j, b->codelist[i]);
 503         exit(1);
 504       }
 505     }
 506   }
 507
 508 }
 509
 510 int vqenc_entry(vqbook *b,double *val){
 511   int ptr=0,k;
 512   while(1){
 513     double c= -b->c[ptr];
 514     double *nptr=b->n+b->dim*ptr;
 515     for(k=0;k<b->dim;k++)
 516       c+=nptr[k]*val[k];
 517     if(c>0.) /* in A */
 518       ptr= -b->ptr0[ptr];
 519     else     /* in B */
 520       ptr= -b->ptr1[ptr];
 521     if(ptr<=0)break;
 522   }
 523   return(-ptr);
 524 }
 525
 526
 527
 528