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.14 2000/01/28 09:05:21 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(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->c->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(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(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<16.*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(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 i,j;
 395   static_codebook *c=(static_codebook *)b->c;
 396   encode_aux *t;
 397
 398   memset(b,0,sizeof(codebook));
 399   memset(c,0,sizeof(static_codebook));
 400   b->c=c;
 401   t=c->encode_tree=calloc(1,sizeof(encode_aux));
 402
 403   for(i=0;i<v->entries;i++)entryindex[i]=i;
 404   for(i=0;i<v->points;i++)pointindex[i]=i;
 405
 406   c->dim=v->elements;
 407   c->entries=v->entries;
 408   c->lengthlist=calloc(c->entries,sizeof(long));
 409
 410   t->alloc=4096;
 411   t->ptr0=malloc(sizeof(long)*t->alloc);
 412   t->ptr1=malloc(sizeof(long)*t->alloc);
 413   t->p=malloc(sizeof(long)*t->alloc);
 414   t->q=malloc(sizeof(long)*t->alloc);
 415
 416   /* first, generate the encoding decision heirarchy */
 417   {
 418     long pointsofar=0;
 419     fprintf(stderr,"Total leaves: %d            \n",
 420             lp_split(v,b,entryindex,v->entries,
 421                      pointindex,v->points,0,&pointsofar));
 422   }
 423
 424   free(entryindex);
 425   free(pointindex);
 426
 427   fprintf(stderr,"Paring and rerouting redundant branches:\n");
 428   /* The tree is likely big and redundant.  Pare and reroute branches */
 429   {
 430     int changedflag=1;
 431     long *unique_entries=alloca(t->aux*sizeof(long));
 432
 433     while(changedflag){
 434       int nodes=0;
 435       changedflag=0;
 436
 437       fprintf(stderr,"\t...");
 438       /* span the tree node by node; list unique decision nodes and
 439          short circuit redundant branches */
 440
 441       for(i=0;i<t->aux;){
 442         int k;
 443
 444         /* check list of unique decisions */
 445         for(j=0;j<nodes;j++)
 446           if(_node_eq(t,i,unique_entries[j]))break;
 447
 448         if(j==nodes){
 449           /* a new entry */
 450           unique_entries[nodes++]=i++;
 451         }else{
 452           /* a redundant entry; find all higher nodes referencing it and
 453              short circuit them to the previously noted unique entry */
 454           changedflag=1;
 455           for(k=0;k<t->aux;k++){
 456             if(t->ptr0[k]==-i)t->ptr0[k]=-unique_entries[j];
 457             if(t->ptr1[k]==-i)t->ptr1[k]=-unique_entries[j];
 458           }
 459
 460           /* Now, we need to fill in the hole from this redundant
 461              entry in the listing.  Insert the last entry in the list.
 462              Fix the forward pointers to that last entry */
 463           t->aux--;
 464           t->ptr0[i]=t->ptr0[t->aux];
 465           t->ptr1[i]=t->ptr1[t->aux];
 466           t->p[i]=t->p[t->aux];
 467           t->q[i]=t->q[t->aux];
 468           for(k=0;k<t->aux;k++){
 469             if(t->ptr0[k]==-t->aux)t->ptr0[k]=-i;
 470             if(t->ptr1[k]==-t->aux)t->ptr1[k]=-i;
 471           }
 472           /* hole plugged */
 473
 474         }
 475       }
 476
 477       fprintf(stderr," %ld remaining\n",t->aux);
 478     }
 479   }
 480
 481   /* run all training points through the decision tree to get a final
 482      probability count */
 483   {
 484     long *probability=malloc(c->entries*sizeof(long));
 485     long *membership=malloc(c->entries*sizeof(long));
 486
 487     for(i=0;i<c->entries;i++)probability[i]=1; /* trivial guard */
 488
 489     b->valuelist=v->entrylist; /* temporary for vqenc_entry; replaced later */
 490     for(i=0;i<v->points;i++){
 491       int ret=codebook_entry(b,v->pointlist+i*v->elements);
 492       probability[ret]++;
 493     }
 494     for(i=0;i<v->entries;i++)membership[i]=i;
 495
 496     /* find codeword lengths */
 497     /* much more elegant means exist.  Brute force n^2, minimum thought */
 498     for(i=v->entries;i>1;i--){
 499       int first=-1,second=-1;
 500       long least=v->points+1;
 501
 502       /* find the two nodes to join */
 503       for(j=0;j<v->entries;j++)
 504         if(probability[j]<least){
 505           least=probability[j];
 506           first=membership[j];
 507         }
 508       least=v->points+1;
 509       for(j=0;j<v->entries;j++)
 510         if(probability[j]<least && membership[j]!=first){
 511           least=probability[j];
 512           second=membership[j];
 513         }
 514       if(first==-1 || second==-1){
 515         fprintf(stderr,"huffman fault; no free branch\n");
 516         exit(1);
 517       }
 518
 519       /* join them */
 520       least=probability[first]+probability[second];
 521       for(j=0;j<v->entries;j++)
 522         if(membership[j]==first || membership[j]==second){
 523           membership[j]=first;
 524           probability[j]=least;
 525           c->lengthlist[j]++;
 526         }
 527     }
 528     for(i=0;i<v->entries-1;i++)
 529       if(membership[i]!=membership[i+1]){
 530         fprintf(stderr,"huffman fault; failed to build single tree\n");
 531         exit(1);
 532       }
 533
 534     /* unneccessary metric */
 535     memset(probability,0,sizeof(long)*v->entries);
 536     for(i=0;i<v->points;i++){
 537       int ret=codebook_entry(b,v->pointlist+i*v->elements);
 538       probability[ret]++;
 539     }
 540
 541     free(probability);
 542     free(membership);
 543   }
 544
 545   /* Sort the entries by codeword length, short to long (eases
 546      assignment and packing to do it now) */
 547   {
 548     long *wordlen=c->lengthlist;
 549     long *index=malloc(c->entries*sizeof(long));
 550     long *revindex=malloc(c->entries*sizeof(long));
 551     int k;
 552     for(i=0;i<c->entries;i++)index[i]=i;
 553     isortvals=c->lengthlist;
 554     qsort(index,c->entries,sizeof(long),iascsort);
 555
 556     /* rearrange storage; ptr0/1 first as it needs a reverse index */
 557     /* n and c stay unchanged */
 558     for(i=0;i<c->entries;i++)revindex[index[i]]=i;
 559     for(i=0;i<t->aux;i++){
 560       if(t->ptr0[i]>=0)t->ptr0[i]=revindex[t->ptr0[i]];
 561       if(t->ptr1[i]>=0)t->ptr1[i]=revindex[t->ptr1[i]];
 562       t->p[i]=revindex[t->p[i]];
 563       t->q[i]=revindex[t->q[i]];
 564     }
 565     free(revindex);
 566
 567     /* map lengthlist and vallist with index */
 568     c->lengthlist=calloc(c->entries,sizeof(long));
 569     b->valuelist=malloc(sizeof(double)*c->entries*c->dim);
 570     c->quantlist=malloc(sizeof(long)*c->entries*c->dim);
 571     for(i=0;i<c->entries;i++){
 572       long e=index[i];
 573       for(k=0;k<c->dim;k++){
 574         b->valuelist[i*c->dim+k]=v->entrylist[e*c->dim+k];
 575         c->quantlist[i*c->dim+k]=quantlist[e*c->dim+k];
 576       }
 577       c->lengthlist[i]=wordlen[e];
 578     }
 579
 580     free(wordlen);
 581   }
 582
 583   fprintf(stderr,"Done.\n\n");
 584 }
 585