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