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-2000             *
   9  * by Monty <monty@xiph.org> and The XIPHOPHORUS Company            *
  10  * http://www.xiph.org/                                             *
  11  *                                                                  *
  12  ********************************************************************
  13
  14  function: build a VQ codebook and the encoding decision 'tree'
  15  last mod: $Id: vqsplit.c,v 1.8 2000/01/05 03:11:12 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 #include "vqgen.h"
  33 #include <sys/time.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 /* Building a codebook from trained set **********************************
  47
  48    The codebook in raw form is technically finished once it's trained.
  49    However, we want to finalize the representative codebook values for
  50    each entry and generate auxiliary information to optimize encoding.
  51    We generate the auxiliary coding tree using collected data,
  52    probably the same data as in the original training */
  53
  54 /* At each recursion, the data set is split in half.  Cells with data
  55    points on side A go into set A, same with set B.  The sets may
  56    overlap.  If the cell overlaps the deviding line only very slightly
  57    (provided parameter), we may choose to ignore the overlap in order
  58    to pare the tree down */
  59
  60 long *isortvals;
  61 int iascsort(const void *a,const void *b){
  62   long av=isortvals[*((long *)a)];
  63   long bv=isortvals[*((long *)b)];
  64   return(av-bv);
  65 }
  66
  67 extern double _dist_sq(vqgen *v,double *a, double *b);
  68
  69 /* goes through the split, but just counts it and returns a metric*/
  70 void vqsp_count(vqgen *v,long *membership,
  71                 long *entryindex,long entries,
  72                 long *pointindex,long points,
  73                 long *entryA,long *entryB,
  74                 double *n, double c,
  75                 long *entriesA,long *entriesB,long *entriesC){
  76   long i,j,k;
  77   long A=0,B=0,C=0;
  78
  79   memset(entryA,0,sizeof(long)*entries);
  80   memset(entryB,0,sizeof(long)*entries);
  81
  82   /* Do the points belonging to this cell occur on sideA, sideB or
  83      both? */
  84
  85   for(i=0;i<points;i++){
  86     double *ppt=_point(v,pointindex[i]);
  87     long   firstentry=membership[i];
  88     double position=-c;
  89
  90     if(!entryA[firstentry] || !entryB[firstentry]){
  91       for(k=0;k<v->elements;k++)
  92       position+=ppt[k]*n[k];
  93       if(position>0.)
  94         entryA[firstentry]=1;
  95       else
  96         entryB[firstentry]=1;
  97     }
  98   }
  99
 100   /* The entry splitting isn't total, so that storage has to be
 101      allocated for recursion.  Reuse the entryA/entryB vectors */
 102   /* keep the entries in ascending order (relative to the original
 103      list); we rely on that stability when ordering p/q choice */
 104   for(j=0;j<entries;j++){
 105     if(entryA[j] && entryB[j])C++;
 106     if(entryA[j])entryA[A++]=entryindex[j];
 107     if(entryB[j])entryB[B++]=entryindex[j];
 108   }
 109   *entriesA=A;
 110   *entriesB=B;
 111   *entriesC=C;
 112 }
 113
 114 void pq_in_out(vqgen *v,double *n,double *c,double *p,double *q){
 115   int k;
 116   *c=0.;
 117   for(k=0;k<v->elements;k++){
 118     double center=(p[k]+q[k])/2.;
 119     n[k]=(center-q[k])*2.;
 120     *c+=center*n[k];
 121   }
 122 }
 123
 124 static void spinnit(int n){
 125   static int p=0;
 126   static long lasttime=0;
 127   long test;
 128   struct timeval thistime;
 129
 130   gettimeofday(&thistime,NULL);
 131   test=thistime.tv_sec*10+thistime.tv_usec/100000;
 132   if(lasttime!=test){
 133     lasttime=test;
 134
 135     fprintf(stderr," %d ",n);
 136
 137     p++;if(p>3)p=0;
 138     switch(p){
 139     case 0:
 140       fprintf(stderr,"|    \r");
 141       break;
 142     case 1:
 143       fprintf(stderr,"/    \r");
 144       break;
 145     case 2:
 146       fprintf(stderr,"-    \r");
 147       break;
 148     case 3:
 149       fprintf(stderr,"\\    \r");
 150       break;
 151     }
 152     fflush(stderr);
 153   }
 154 }
 155
 156 int lp_split(vqgen *v,vqbook *b,
 157              long *entryindex,long entries,
 158              long *pointindex,long points,
 159              long depth){
 160
 161   /* The encoder, regardless of book, will be using a straight
 162      euclidian distance-to-point metric to determine closest point.
 163      Thus we split the cells using the same (we've already trained the
 164      codebook set spacing and distribution using special metrics and
 165      even a midpoint division won't disturb the basic properties) */
 166
 167   long ret;
 168   double *n=alloca(sizeof(double)*v->elements);
 169   double c;
 170   long *entryA=calloc(entries,sizeof(long));
 171   long *entryB=calloc(entries,sizeof(long));
 172   long entriesA=0;
 173   long entriesB=0;
 174   long entriesC=0;
 175   long pointsA=0;
 176   long i,j,k;
 177
 178   long *membership=malloc(sizeof(long)*points);
 179
 180   long   besti=-1;
 181   long   bestj=-1;
 182
 183   /* which cells do points belong to?  Do this before n^2 best pair chooser. */
 184
 185   for(i=0;i<points;i++){
 186     double *ppt=_point(v,pointindex[i]);
 187     long   firstentry=0;
 188     double firstmetric=_dist_sq(v,_now(v,entryindex[0]),ppt);
 189
 190     if(points*entries>64*1024)spinnit(entries);
 191
 192     for(j=1;j<entries;j++){
 193       double thismetric=_dist_sq(v,_now(v,entryindex[j]),ppt);
 194       if(thismetric<firstmetric){
 195         firstmetric=thismetric;
 196         firstentry=j;
 197       }
 198     }
 199
 200     membership[i]=firstentry;
 201   }
 202
 203   /* We need to find the dividing hyperplane. find the median of each
 204      axis as the centerpoint and the normal facing farthest point */
 205
 206   /* more than one way to do this part.  For small sets, we can brute
 207      force it. */
 208
 209   if(entries<8 || (double)points*entries*entries<128.*1024*1024){
 210     /* try every pair possibility */
 211     double best=0;
 212     double this;
 213     for(i=0;i<entries-1;i++){
 214       for(j=i+1;j<entries;j++){
 215         spinnit(entries-i);
 216         pq_in_out(v,n,&c,_now(v,entryindex[i]),_now(v,entryindex[j]));
 217         vqsp_count(v,membership,
 218                    entryindex,entries,
 219                    pointindex,points,
 220                    entryA,entryB,
 221                    n, c,
 222                    &entriesA,&entriesB,&entriesC);
 223         this=(entriesA-entriesC)*(entriesB-entriesC);
 224
 225         /* when choosing best, we also want some form of stability to
 226            make sure more branches are pared later; secondary
 227            weighting isn;t needed as the entry lists are in ascending
 228            order, and we always try p/q in the same sequence */
 229
 230         if( (besti==-1) ||
 231             (this>best) ){
 232
 233           best=this;
 234           besti=entryindex[i];
 235           bestj=entryindex[j];
 236
 237         }
 238       }
 239     }
 240   }else{
 241     double *p=alloca(v->elements*sizeof(double));
 242     double *q=alloca(v->elements*sizeof(double));
 243     double best=0.;
 244
 245     /* try COG/normal and furthest pairs */
 246     /* meanpoint */
 247     /* eventually, we want to select the closest entry and figure n/c
 248        from p/q (because storing n/c is too large */
 249     for(k=0;k<v->elements;k++){
 250       spinnit(entries);
 251
 252       p[k]=0.;
 253       for(j=0;j<entries;j++)
 254         p[k]+=v->entrylist[entryindex[j]*v->elements+k];
 255       p[k]/=entries;
 256
 257     }
 258
 259     /* we go through the entries one by one, looking for the entry on
 260        the other side closest to the point of reflection through the
 261        center */
 262
 263     for(i=0;i<entries;i++){
 264       double *ppi=_now(v,entryindex[i]);
 265       double ref_best=0.;
 266       double ref_j=-1;
 267       double this;
 268       spinnit(entries-i);
 269
 270       for(k=0;k<v->elements;k++)
 271         q[k]=2*p[k]-ppi[k];
 272
 273       for(j=0;j<entries;j++){
 274         if(j!=i){
 275           double this=_dist_sq(v,q,_now(v,entryindex[j]));
 276           if(ref_j==-1 || this<ref_best){
 277             ref_best=this;
 278             ref_j=entryindex[j];
 279           }
 280         }
 281       }
 282
 283       pq_in_out(v,n,&c,ppi,_now(v,ref_j));
 284       vqsp_count(v,membership,
 285                  entryindex,entries,
 286                  pointindex,points,
 287                  entryA,entryB,
 288                  n, c,
 289                  &entriesA,&entriesB,&entriesC);
 290       this=(entriesA-entriesC)*(entriesB-entriesC);
 291
 292         /* when choosing best, we also want some form of stability to
 293            make sure more branches are pared later; secondary
 294            weighting isn;t needed as the entry lists are in ascending
 295            order, and we always try p/q in the same sequence */
 296
 297       if( (besti==-1) ||
 298           (this>best) ){
 299
 300         best=this;
 301         besti=entryindex[i];
 302         bestj=ref_j;
 303
 304       }
 305     }
 306     if(besti>bestj){
 307       long temp=besti;
 308       besti=bestj;
 309       bestj=temp;
 310     }
 311
 312   }
 313
 314   /* find cells enclosing points */
 315   /* count A/B points */
 316
 317   pq_in_out(v,n,&c,_now(v,besti),_now(v,bestj));
 318   vqsp_count(v,membership,
 319              entryindex,entries,
 320              pointindex,points,
 321              entryA,entryB,
 322              n, c,
 323              &entriesA,&entriesB,&entriesC);
 324
 325   free(membership);
 326
 327   /* the point index is split, so we do an Order n rearrangement into
 328      A first/B last and just pass it on */
 329   {
 330     long Aptr=0;
 331     long Bptr=points-1;
 332     while(Aptr<=Bptr){
 333       while(Aptr<=Bptr){
 334         double position=-c;
 335         for(k=0;k<v->elements;k++)
 336           position+=_point(v,pointindex[Aptr])[k]*n[k];
 337         if(position<=0.)break; /* not in A */
 338         Aptr++;
 339       }
 340       while(Aptr<=Bptr){
 341         double position=-c;
 342         for(k=0;k<v->elements;k++)
 343           position+=_point(v,pointindex[Bptr])[k]*n[k];
 344         if(position>0.)break; /* not in B */
 345         Bptr--;
 346       }
 347       if(Aptr<Bptr){
 348         long temp=pointindex[Aptr];
 349         pointindex[Aptr]=pointindex[Bptr];
 350         pointindex[Bptr]=temp;
 351       }
 352       pointsA=Aptr;
 353     }
 354   }
 355
 356   fprintf(stderr,"split: total=%ld depth=%ld set A=%ld:%ld:%ld=B\n",
 357           entries,depth,entriesA-entriesC,entriesC,entriesB-entriesC);
 358   {
 359     long thisaux=b->aux++;
 360     if(b->aux>=b->alloc){
 361       b->alloc*=2;
 362       b->ptr0=realloc(b->ptr0,sizeof(long)*b->alloc);
 363       b->ptr1=realloc(b->ptr1,sizeof(long)*b->alloc);
 364       b->p=realloc(b->p,sizeof(long)*b->alloc);
 365       b->q=realloc(b->q,sizeof(long)*b->alloc);
 366     }
 367
 368     b->p[thisaux]=besti;
 369     b->q[thisaux]=bestj;
 370
 371     if(entriesA==1){
 372       ret=1;
 373       b->ptr0[thisaux]=entryA[0];
 374     }else{
 375       b->ptr0[thisaux]= -b->aux;
 376       ret=lp_split(v,b,entryA,entriesA,pointindex,pointsA,depth+1);
 377     }
 378     if(entriesB==1){
 379       ret++;
 380       b->ptr1[thisaux]=entryB[0];
 381     }else{
 382       b->ptr1[thisaux]= -b->aux;
 383       ret+=lp_split(v,b,entryB,entriesB,pointindex+pointsA,points-pointsA,
 384                     depth+1);
 385     }
 386   }
 387   free(entryA);
 388   free(entryB);
 389   return(ret);
 390 }
 391
 392 static int _node_eq(vqbook *v, long a, long b){
 393   long    Aptr0=v->ptr0[a];
 394   long    Aptr1=v->ptr1[a];
 395   long    Ap   =v->p[a];
 396   long    Aq   =v->q[a];
 397   long    Bptr0=v->ptr0[b];
 398   long    Bptr1=v->ptr1[b];
 399   long    Bp   =v->p[b];
 400   long    Bq   =v->q[b];
 401   int i;
 402
 403   /* the possibility of choosing the same p and q, but switched, can;t
 404      happen because we always look for the best p/q in the same search
 405      order and the search is stable */
 406
 407   if(Aptr0==Bptr0 && Aptr1==Bptr1 && Ap==Bp && Aq==Bq)
 408     return(1);
 409
 410   return(0);
 411 }
 412
 413 void vqsp_book(vqgen *v, vqbook *b, long *quantlist){
 414   long *entryindex=malloc(sizeof(long)*v->entries);
 415   long *pointindex=malloc(sizeof(long)*v->points);
 416   long *membership=malloc(sizeof(long)*v->points);
 417   long i,j;
 418
 419   memset(b,0,sizeof(vqbook));
 420   for(i=0;i<v->entries;i++)entryindex[i]=i;
 421   for(i=0;i<v->points;i++)pointindex[i]=i;
 422
 423   b->dim=v->elements;
 424   b->entries=v->entries;
 425   b->alloc=4096;
 426
 427   b->ptr0=malloc(sizeof(long)*b->alloc);
 428   b->ptr1=malloc(sizeof(long)*b->alloc);
 429   b->p=malloc(sizeof(long)*b->alloc);
 430   b->q=malloc(sizeof(long)*b->alloc);
 431   b->lengthlist=calloc(b->entries,sizeof(long));
 432
 433   /* which cells do points belong to?  Only do this once. */
 434
 435   for(i=0;i<v->points;i++){
 436     double *ppt=_point(v,i);
 437     long   firstentry=0;
 438     double firstmetric=_dist_sq(v,_now(v,0),ppt);
 439
 440     for(j=1;j<v->entries;j++){
 441       double thismetric=_dist_sq(v,_now(v,j),ppt);
 442       if(thismetric<firstmetric){
 443         firstmetric=thismetric;
 444         firstentry=j;
 445       }
 446     }
 447
 448     membership[i]=firstentry;
 449   }
 450
 451
 452   /* first, generate the encoding decision heirarchy */
 453   fprintf(stderr,"Total leaves: %d\n",
 454           lp_split(v,b,entryindex,v->entries,
 455                    pointindex,v->points,0));
 456
 457   free(entryindex);
 458   free(pointindex);
 459
 460   fprintf(stderr,"Paring and rerouting redundant branches:\n");
 461   /* The tree is likely big and redundant.  Pare and reroute branches */
 462   {
 463     int changedflag=1;
 464     long *unique_entries=alloca(b->aux*sizeof(long));
 465
 466     while(changedflag){
 467       int nodes=0;
 468       changedflag=0;
 469
 470       fprintf(stderr,"\t...");
 471       /* span the tree node by node; list unique decision nodes and
 472          short circuit redundant branches */
 473
 474       for(i=0;i<b->aux;){
 475         int k;
 476
 477         /* check list of unique decisions */
 478         for(j=0;j<nodes;j++)
 479           if(_node_eq(b,i,unique_entries[j]))break;
 480
 481         if(j==nodes){
 482           /* a new entry */
 483           unique_entries[nodes++]=i++;
 484         }else{
 485           /* a redundant entry; find all higher nodes referencing it and
 486              short circuit them to the previously noted unique entry */
 487           changedflag=1;
 488           for(k=0;k<b->aux;k++){
 489             if(b->ptr0[k]==-i)b->ptr0[k]=-unique_entries[j];
 490             if(b->ptr1[k]==-i)b->ptr1[k]=-unique_entries[j];
 491           }
 492
 493           /* Now, we need to fill in the hole from this redundant
 494              entry in the listing.  Insert the last entry in the list.
 495              Fix the forward pointers to that last entry */
 496           b->aux--;
 497           b->ptr0[i]=b->ptr0[b->aux];
 498           b->ptr1[i]=b->ptr1[b->aux];
 499           b->p[i]=b->p[b->aux];
 500           b->q[i]=b->q[b->aux];
 501           for(k=0;k<b->aux;k++){
 502             if(b->ptr0[k]==-b->aux)b->ptr0[k]=-i;
 503             if(b->ptr1[k]==-b->aux)b->ptr1[k]=-i;
 504           }
 505           /* hole plugged */
 506
 507         }
 508       }
 509
 510       fprintf(stderr," %ld remaining\n",b->aux);
 511     }
 512   }
 513
 514   /* run all training points through the decision tree to get a final
 515      probability count */
 516   {
 517     long *probability=malloc(b->entries*sizeof(long));
 518     long *membership=malloc(b->entries*sizeof(long));
 519
 520     for(i=0;i<b->entries;i++)probability[i]=1; /* trivial guard */
 521
 522     b->valuelist=v->entrylist; /* temporary for vqenc_entry; replaced later */
 523     for(i=0;i<v->points;i++){
 524       int ret=vqenc_entry(b,v->pointlist+i*v->elements);
 525       probability[ret]++;
 526     }
 527     for(i=0;i<v->entries;i++)membership[i]=i;
 528
 529     /* find codeword lengths */
 530     /* much more elegant means exist.  Brute force n^2, minimum thought */
 531     for(i=v->entries;i>1;i--){
 532       int first=-1,second=-1;
 533       long least=v->points+1;
 534
 535       /* find the two nodes to join */
 536       for(j=0;j<v->entries;j++)
 537         if(probability[j]<least){
 538           least=probability[j];
 539           first=membership[j];
 540         }
 541       least=v->points+1;
 542       for(j=0;j<v->entries;j++)
 543         if(probability[j]<least && membership[j]!=first){
 544           least=probability[j];
 545           second=membership[j];
 546         }
 547       if(first==-1 || second==-1){
 548         fprintf(stderr,"huffman fault; no free branch\n");
 549         exit(1);
 550       }
 551
 552       /* join them */
 553       least=probability[first]+probability[second];
 554       for(j=0;j<v->entries;j++)
 555         if(membership[j]==first || membership[j]==second){
 556           membership[j]=first;
 557           probability[j]=least;
 558           b->lengthlist[j]++;
 559         }
 560     }
 561     for(i=0;i<v->entries-1;i++)
 562       if(membership[i]!=membership[i+1]){
 563         fprintf(stderr,"huffman fault; failed to build single tree\n");
 564         exit(1);
 565       }
 566
 567     /* unneccessary metric */
 568     memset(probability,0,sizeof(long)*v->entries);
 569     for(i=0;i<v->points;i++){
 570       int ret=vqenc_entry(b,v->pointlist+i*v->elements);
 571       probability[ret]++;
 572     }
 573
 574     free(probability);
 575     free(membership);
 576   }
 577
 578   /* Sort the entries by codeword length, short to long (eases
 579      assignment and packing to do it now) */
 580   {
 581     long *wordlen=b->lengthlist;
 582     long *index=malloc(b->entries*sizeof(long));
 583     long *revindex=malloc(b->entries*sizeof(long));
 584     int k;
 585     for(i=0;i<b->entries;i++)index[i]=i;
 586     isortvals=b->lengthlist;
 587     qsort(index,b->entries,sizeof(long),iascsort);
 588
 589     /* rearrange storage; ptr0/1 first as it needs a reverse index */
 590     /* n and c stay unchanged */
 591     for(i=0;i<b->entries;i++)revindex[index[i]]=i;
 592     for(i=0;i<b->aux;i++){
 593       if(b->ptr0[i]>=0)b->ptr0[i]=revindex[b->ptr0[i]];
 594       if(b->ptr1[i]>=0)b->ptr1[i]=revindex[b->ptr1[i]];
 595       b->p[i]=revindex[b->p[i]];
 596       b->q[i]=revindex[b->q[i]];
 597     }
 598     free(revindex);
 599
 600     /* map lengthlist and vallist with index */
 601     b->lengthlist=calloc(b->entries,sizeof(long));
 602     b->valuelist=malloc(sizeof(double)*b->entries*b->dim);
 603     b->quantlist=malloc(sizeof(long)*b->entries*b->dim);
 604     for(i=0;i<b->entries;i++){
 605       long e=index[i];
 606       for(k=0;k<b->dim;k++){
 607         b->valuelist[i*b->dim+k]=v->entrylist[e*b->dim+k];
 608         b->quantlist[i*b->dim+k]=quantlist[e*b->dim+k];
 609       }
 610       b->lengthlist[i]=wordlen[e];
 611     }
 612
 613     free(wordlen);
 614   }
 615
 616   /* generate the codewords (short to long) */
 617   {
 618     long current=0;
 619     long length=0;
 620     b->codelist=malloc(sizeof(long)*b->entries);
 621     for(i=0;i<b->entries;i++){
 622       if(length != b->lengthlist[i]){
 623         current<<=(b->lengthlist[i]-length);
 624         length=b->lengthlist[i];
 625       }
 626       b->codelist[i]=current;
 627       current++;
 628     }
 629   }
 630
 631   /* sanity check the codewords */
 632   for(i=0;i<b->entries;i++){
 633     for(j=i+1;j<b->entries;j++){
 634       if(b->codelist[i]==b->codelist[j]){
 635         fprintf(stderr,"Error; codewords for %ld and %ld both equal %lx\n",
 636                 i, j, b->codelist[i]);
 637         exit(1);
 638       }
 639     }
 640   }
 641   fprintf(stderr,"Done.\n\n");
 642 }
 643
 644 /* slow version for use here that does not use a preexpanded n/c. */
 645 int vqenc_entry(vqbook *b,double *val){
 646   int ptr=0,k;
 647   double *n=alloca(b->dim*sizeof(double));
 648
 649   while(1){
 650     double c=0.;
 651     double *p=b->valuelist+b->p[ptr]*b->dim;
 652     double *q=b->valuelist+b->q[ptr]*b->dim;
 653
 654     for(k=0;k<b->dim;k++){
 655       n[k]=p[k]-q[k];
 656       c-=(p[k]+q[k])*n[k];
 657     }
 658     c/=2.;
 659
 660     for(k=0;k<b->dim;k++)
 661       c+=n[k]*val[k];
 662     if(c>0.) /* in A */
 663       ptr= -b->ptr0[ptr];
 664     else     /* in B */
 665       ptr= -b->ptr1[ptr];
 666     if(ptr<=0)break;
 667   }
 668   return(-ptr);
 669 }
 670
 671
 672
 673