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