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-1999             *
   9  * by 1999 Monty <monty@xiph.org> and The XIPHOPHORUS Company       *
  10  * http://www.xiph.org/                                             *
  11  *                                                                  *
  12  ********************************************************************
  13
  14  function: build a VQ codebook
  15  author: Monty <xiphmont@mit.edu>
  16  modifications by: Monty
  17  last modification date: Dec 10 1999
  18
  19  ********************************************************************/
  20
  21 /* This code is *not* part of libvorbis.  It is used to generate
  22    trained codebooks offline and then spit the results into a
  23    pregenerated codebook that is compiled into libvorbis.  It is an
  24    expensive (but good) algorithm.  Run it on big iron. */
  25
  26 /* There are so many optimizations to explore in *both* stages that
  27    considering the undertaking is almost withering.  For now, we brute
  28    force it all */
  29
  30 #include <stdlib.h>
  31 #include <stdio.h>
  32 #include <math.h>
  33 #include <string.h>
  34 #include "vqgen.h"
  35 #include <sys/time.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 double *sortvals;
  63 int els;
  64 int iascsort(const void *a,const void *b){
  65   double av=sortvals[*((long *)a) * els];
  66   double bv=sortvals[*((long *)b) * els];
  67   if(av<bv)return(-1);
  68   return(1);
  69 }
  70
  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 *entryindex,long entries,
  75                 long *pointindex,long points,
  76                 long *entryA,long *entryB,
  77                 double *n, double c, double slack,
  78                 long *entriesA,long *entriesB,long *entriesC){
  79   long i,j,k;
  80   long A=0,B=0,C=0;
  81
  82   memset(entryA,0,sizeof(long)*entries);
  83   memset(entryB,0,sizeof(long)*entries);
  84
  85   for(i=0;i<points;i++){
  86     double *ppt=_point(v,pointindex[i]);
  87     long   firstentry=0;
  88     double firstmetric=_dist_sq(v,_now(v,entryindex[0]),ppt);
  89     double position=-c;
  90
  91     for(j=1;j<entries;j++){
  92       double thismetric=_dist_sq(v,_now(v,entryindex[j]),ppt);
  93       if(thismetric<firstmetric){
  94         firstmetric=thismetric;
  95         firstentry=j;
  96       }
  97     }
  98
  99     /* count point split */
 100     for(k=0;k<v->elements;k++)
 101       position+=ppt[k]*n[k];
 102     if(position>0.){
 103       entryA[firstentry]++;
 104     }else{
 105       entryB[firstentry]++;
 106     }
 107   }
 108
 109   /* look to see if entries are in the slack zone */
 110   /* The entry splitting isn't total, so that storage has to be
 111      allocated for recursion.  Reuse the entryA/entryB vectors */
 112   for(j=0;j<entries;j++){
 113     long total=entryA[j]+entryB[j];
 114     if((double)entryA[j]/total<slack){
 115       entryA[j]=0;
 116     }else if((double)entryB[j]/total<slack){
 117       entryB[j]=0;
 118     }
 119     if(entryA[j] && entryB[j])C++;
 120     if(entryA[j])entryA[A++]=entryindex[j];
 121     if(entryB[j])entryB[B++]=entryindex[j];
 122   }
 123   *entriesA=A;
 124   *entriesB=B;
 125   *entriesC=C;
 126 }
 127
 128 void pq_in_out(vqgen *v,double *n,double *c,double *p,double *q){
 129   int k;
 130   *c=0.;
 131   for(k=0;k<v->elements;k++){
 132     double center=(p[k]+q[k])/2.;
 133     n[k]=(center-q[k])*2.;
 134     *c+=center*n[k];
 135   }
 136 }
 137
 138 void pq_center_out(vqgen *v,double *n,double *c,double *center,double *q){
 139   int k;
 140   *c=0.;
 141   for(k=0;k<v->elements;k++){
 142     n[k]=(center[k]-q[k])*2.;
 143     *c+=center[k]*n[k];
 144   }
 145 }
 146
 147 static void spinnit(void){
 148   static int p=0;
 149   static long lasttime=0;
 150   long test;
 151   struct timeval thistime;
 152
 153   gettimeofday(&thistime,NULL);
 154   test=thistime.tv_sec*10+thistime.tv_usec/100000;
 155   if(lasttime!=test){
 156     lasttime=test;
 157
 158     p++;if(p>3)p=0;
 159     switch(p){
 160     case 0:
 161       fprintf(stderr,"|\b");
 162       break;
 163     case 1:
 164       fprintf(stderr,"/\b");
 165       break;
 166     case 2:
 167       fprintf(stderr,"-\b");
 168       break;
 169     case 3:
 170       fprintf(stderr,"\\\b");
 171       break;
 172     }
 173     fflush(stderr);
 174   }
 175 }
 176
 177 int lp_split(vqgen *v,vqbook *b,
 178              long *entryindex,long entries,
 179              long *pointindex,long points,
 180              long depth,double slack){
 181
 182   /* The encoder, regardless of book, will be using a straight
 183      euclidian distance-to-point metric to determine closest point.
 184      Thus we split the cells using the same (we've already trained the
 185      codebook set spacing and distribution using special metrics and
 186      even a midpoint division won't disturb the basic properties) */
 187
 188   long ret;
 189   double *p;
 190   double *q;
 191   double *n;
 192   double c;
 193   long *entryA=calloc(entries,sizeof(long));
 194   long *entryB=calloc(entries,sizeof(long));
 195   long entriesA=0;
 196   long entriesB=0;
 197   long entriesC=0;
 198   long pointsA=0;
 199   long i,j,k;
 200
 201   p=alloca(sizeof(double)*v->elements);
 202   q=alloca(sizeof(double)*v->elements);
 203   n=alloca(sizeof(double)*v->elements);
 204   memset(p,0,sizeof(double)*v->elements);
 205
 206   /* We need to find the dividing hyperplane. find the median of each
 207      axis as the centerpoint and the normal facing farthest point */
 208
 209   /* more than one way to do this part.  For small sets, we can brute
 210      force it. */
 211
 212   if(entries<32){
 213     /* try every pair possibility */
 214     double best=0;
 215     long   besti=0;
 216     long   bestj=0;
 217     double this;
 218     for(i=0;i<entries-1;i++){
 219       for(j=i+1;j<entries;j++){
 220         spinnit();
 221         pq_in_out(v,n,&c,_now(v,entryindex[i]),_now(v,entryindex[j]));
 222         vqsp_count(v,entryindex,entries,
 223                  pointindex,points,
 224                  entryA,entryB,
 225                  n, c, slack,
 226                  &entriesA,&entriesB,&entriesC);
 227         this=(entriesA-entriesC)*(entriesB-entriesC);
 228
 229         if(this>best){
 230           best=this;
 231           besti=i;
 232           bestj=j;
 233         }
 234       }
 235     }
 236     pq_in_out(v,n,&c,_now(v,entryindex[besti]),_now(v,entryindex[bestj]));
 237   }else{
 238     double best=0.;
 239     long   bestj=0;
 240
 241     /* try COG/normal and furthest pairs */
 242     /* medianpoint */
 243     for(k=0;k<v->elements;k++){
 244       spinnit();
 245
 246       /* just sort the index array */
 247       sortvals=v->entrylist+k;
 248       els=v->elements;
 249       qsort(entryindex,entries,sizeof(long),iascsort);
 250       if(entries&0x1){
 251         p[k]=v->entrylist[(entryindex[entries/2])*v->elements+k];
 252       }else{
 253         p[k]=(v->entrylist[(entryindex[entries/2])*v->elements+k]+
 254               v->entrylist[(entryindex[entries/2-1])*v->elements+k])/2.;
 255       }
 256     }
 257
 258     spinnit();
 259
 260     /* try every normal, but just for distance */
 261     for(j=0;j<entries;j++){
 262       double *ppj=_now(v,entryindex[j]);
 263       double this=_dist_sq(v,p,ppj);
 264       if(this>best){
 265         best=this;
 266         bestj=j;
 267       }
 268     }
 269
 270     pq_center_out(v,n,&c,p,_point(v,pointindex[bestj]));
 271
 272
 273   }
 274
 275   /* find cells enclosing points */
 276   /* count A/B points */
 277
 278   vqsp_count(v,entryindex,entries,
 279            pointindex,points,
 280            entryA,entryB,
 281            n, c, slack,
 282            &entriesA,&entriesB,&entriesC);
 283
 284   /* the point index is split evenly, so we do an Order n
 285      rearrangement into A first/B last and just pass it on */
 286   {
 287     long Aptr=0;
 288     long Bptr=points-1;
 289     while(Aptr<=Bptr){
 290       while(Aptr<=Bptr){
 291         double position=-c;
 292         for(k=0;k<v->elements;k++)
 293           position+=_point(v,pointindex[Aptr])[k]*n[k];
 294         if(position<=0.)break; /* not in A */
 295         Aptr++;
 296       }
 297       while(Aptr<=Bptr){
 298         double position=-c;
 299         for(k=0;k<v->elements;k++)
 300           position+=_point(v,pointindex[Bptr])[k]*n[k];
 301         if(position>0.)break; /* not in B */
 302         Bptr--;
 303       }
 304       if(Aptr<Bptr){
 305         long temp=pointindex[Aptr];
 306         pointindex[Aptr]=pointindex[Bptr];
 307         pointindex[Bptr]=temp;
 308       }
 309       pointsA=Aptr;
 310     }
 311   }
 312
 313   fprintf(stderr,"split: total=%ld depth=%ld set A=%ld:%ld:%ld=B\n",
 314           entries,depth,entriesA-entriesC,entriesC,entriesB-entriesC);
 315   {
 316     long thisaux=b->aux++;
 317     if(b->aux>=b->alloc){
 318       b->alloc*=2;
 319       b->ptr0=realloc(b->ptr0,sizeof(long)*b->alloc);
 320       b->ptr1=realloc(b->ptr1,sizeof(long)*b->alloc);
 321       b->n=realloc(b->n,sizeof(double)*b->elements*b->alloc);
 322       b->c=realloc(b->c,sizeof(double)*b->alloc);
 323     }
 324
 325     memcpy(b->n+b->elements*thisaux,n,sizeof(double)*v->elements);
 326     b->c[thisaux]=c;
 327
 328     if(entriesA==1){
 329       ret=1;
 330       b->ptr0[thisaux]=entryA[0];
 331     }else{
 332       b->ptr0[thisaux]= -b->aux;
 333       ret=lp_split(v,b,entryA,entriesA,pointindex,pointsA,depth+1,slack);
 334     }
 335     if(entriesB==1){
 336       ret++;
 337       b->ptr1[thisaux]=entryB[0];
 338     }else{
 339       b->ptr1[thisaux]= -b->aux;
 340       ret+=lp_split(v,b,entryB,entriesB,pointindex+pointsA,points-pointsA,
 341                    depth+1,slack);
 342     }
 343   }
 344   free(entryA);
 345   free(entryB);
 346   return(ret);
 347 }
 348
 349 void vqsp_book(vqgen *v, vqbook *b){
 350   long *entryindex=malloc(sizeof(double)*v->entries);
 351   long *pointindex=malloc(sizeof(double)*v->points);
 352   long i;
 353
 354   memset(b,0,sizeof(vqbook));
 355   for(i=0;i<v->entries;i++)entryindex[i]=i;
 356   for(i=0;i<v->points;i++)pointindex[i]=i;
 357
 358   b->elements=v->elements;
 359   b->entries=v->entries;
 360   b->alloc=4096;
 361   b->ptr0=malloc(sizeof(long)*b->alloc);
 362   b->ptr1=malloc(sizeof(long)*b->alloc);
 363   b->n=malloc(sizeof(double)*b->elements*b->alloc);
 364   b->c=malloc(sizeof(double)*b->alloc);
 365
 366   b->valuelist=v->entrylist;
 367   b->codelist=malloc(sizeof(long)*b->entries);
 368   b->lengthlist=calloc(b->entries,sizeof(long));
 369
 370   /* first, generate the encoding decision heirarchy */
 371   fprintf(stderr,"Total leaves: %d\n",
 372           lp_split(v,b,entryindex,v->entries, pointindex,v->points,0,0));
 373
 374   /* run all training points through the decision tree to get a final
 375      probability count */
 376   {
 377     long *probability=calloc(b->entries,sizeof(long));
 378     long *membership=malloc(b->entries*sizeof(long));
 379     long j;
 380
 381     for(i=0;i<v->points;i++){
 382       int ret=vqenc_entry(b,v->pointlist+i*v->elements);
 383       probability[ret]++;
 384     }
 385     for(i=0;i<v->entries;i++)membership[i]=i;
 386
 387     /* find codeword lengths */
 388     /* much more elegant means exist.  Brute force n^2, minimum thought */
 389     for(i=v->entries;i>1;i--){
 390       int first=-1,second=-1;
 391       long least=v->points+1;
 392
 393       /* find the two nodes to join */
 394       for(j=0;j<v->entries;j++)
 395         if(probability[j]<least){
 396           least=probability[j];
 397           first=membership[j];
 398         }
 399       least=v->points+1;
 400       for(j=0;j<v->entries;j++)
 401         if(probability[j]<least && membership[j]!=first){
 402           least=probability[j];
 403           second=membership[j];
 404         }
 405       if(first==-1 || second==-1){
 406         fprintf(stderr,"huffman fault; no free branch\n");
 407         exit(1);
 408       }
 409
 410       /* join them */
 411       least=probability[first]+probability[second];
 412       for(j=0;j<v->entries;j++)
 413         if(membership[j]==first || membership[j]==second){
 414           membership[j]=first;
 415           probability[j]=least;
 416           b->lengthlist[j]++;
 417         }
 418     }
 419     for(i=0;i<v->entries-1;i++)
 420       if(membership[i]!=membership[i+1]){
 421         fprintf(stderr,"huffman fault; failed to build single tree\n");
 422         exit(1);
 423       }
 424
 425     /* unneccessary metric */
 426     memset(probability,0,sizeof(long)*v->entries);
 427     for(i=0;i<v->points;i++){
 428       int ret=vqenc_entry(b,v->pointlist+i*v->elements);
 429       probability[ret]++;
 430     }
 431
 432     /* print the results */
 433     for(i=0;i<v->entries;i++)
 434       fprintf(stderr,"%d: count=%ld codelength=%d\n",i,probability[i],b->lengthlist[i]);
 435
 436     free(probability);
 437     free(membership);
 438   }
 439
 440   /* rearrange, sort by codelength */
 441
 442
 443   /* generate the codewords (short to long) */
 444
 445
 446   free(entryindex);
 447   free(pointindex);
 448 }
 449
 450 int vqenc_entry(vqbook *b,double *val){
 451   int ptr=0,k;
 452   while(1){
 453     double c= -b->c[ptr];
 454     double *nptr=b->n+b->elements*ptr;
 455     for(k=0;k<b->elements;k++)
 456       c+=nptr[k]*val[k];
 457     if(c>0.) /* in A */
 458       ptr= -b->ptr0[ptr];
 459     else     /* in B */
 460       ptr= -b->ptr1[ptr];
 461     if(ptr<=0)break;
 462   }
 463   return(-ptr);
 464 }
 465
 466
 467
 468