libvorbis-1.0.1/vq/vqsplit.c

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