vq/vqgen.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: Nov 18 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 #include <stdlib.h>
  27 #include <stdio.h>
  28 #include <math.h>
  29 #include <string.h>
  30
  31 /************************************************************************
  32  * The basic idea here is that a VQ codebook is like an m-dimensional
  33  * foam with n bubbles.  The bubbles compete for space/volume and are
  34  * 'pressurized' [biased] according to some metric.  The basic alg
  35  * iterates through allowing the bubbles to compete for space until
  36  * they converge (if the damping is dome properly) on a steady-state
  37  * solution.
  38  *
  39  * We use the ratio of local to average 'error' as the metric to bias a
  40  * variable-length word codebook, and probability of occurrence within
  41  * that bubble as the metric to bias fixed length word
  42  * codebooks. Individual input points, collected from libvorbis, are
  43  * used to train the algorithm monte-carlo style.  */
  44
  45 typedef struct vqgen{
  46   int    elements;
  47   double errspread;
  48
  49   /* point cache */
  50   double *pointlist;
  51   long   *first;
  52   long   *second;
  53   long   points;
  54   long   allocated;
  55
  56   /* entries */
  57   double *entry_now;
  58   double *entry_next;
  59   long   *assigned;
  60   double *metric;
  61   double *bias;
  62   long   entries;
  63
  64   double (*metric_func)   (struct vqgen *v,double *a,double *b);
  65 } vqgen;
  66
  67 /* internal helpers *****************************************************/
  68 double *_point(vqgen *v,long ptr){
  69   return v->pointlist+(v->elements*ptr);
  70 }
  71
  72 double *_now(vqgen *v,long ptr){
  73   return v->entry_now+(v->elements*ptr);
  74 }
  75
  76 double *_next(vqgen *v,long ptr){
  77   return v->entry_next+(v->elements*ptr);
  78 }
  79
  80 double _dist_sq(vqgen *v,double *a, double *b){
  81   int i;
  82   int el=v->elements;
  83   double acc=0.;
  84   for(i=0;i<el;i++){
  85     double val=(a[i]-b[i]);
  86     acc+=val*val;
  87   }
  88   return acc;
  89 }
  90
  91 void vqgen_init(vqgen *v,int elements,int entries,
  92                 double (*metric)(vqgen *,double *, double *),
  93                 double spread){
  94   memset(v,0,sizeof(vqgen));
  95
  96   v->elements=elements;
  97   v->errspread=spread;
  98   v->allocated=32768;
  99   v->pointlist=malloc(v->allocated*v->elements*sizeof(double));
 100   v->first=malloc(v->allocated*sizeof(long));
 101   v->second=malloc(v->allocated*sizeof(long));
 102
 103   v->entries=entries;
 104   v->entry_now=malloc(v->entries*v->elements*sizeof(double));
 105   v->entry_next=malloc(v->entries*v->elements*sizeof(double));
 106   v->assigned=malloc(v->entries*sizeof(long));
 107   v->metric=malloc(v->entries*sizeof(double));
 108   v->bias=calloc(v->entries,sizeof(double));
 109   if(metric)
 110     v->metric_func=metric;
 111   else
 112     v->metric_func=_dist_sq;
 113 }
 114
 115 /* *must* be beefed up.  Perhaps a Floyd-Steinberg like scattering? */
 116 void _vqgen_seed(vqgen *v){
 117   memcpy(v->entry_now,v->pointlist,sizeof(double)*v->entries*v->elements);
 118 }
 119
 120 void vqgen_addpoint(vqgen *v, double *p){
 121   if(v->points>=v->allocated){
 122     v->allocated*=2;
 123     v->pointlist=realloc(v->pointlist,v->allocated*v->elements*sizeof(double));
 124     v->first=realloc(v->first,v->allocated*sizeof(long));
 125     v->second=realloc(v->second,v->allocated*sizeof(long));
 126   }
 127
 128   memcpy(_point(v,v->points),p,sizeof(double)*v->elements);
 129   v->points++;
 130   if(v->points==v->entries)_vqgen_seed(v);
 131 }
 132
 133 double *sort_t_vals;
 134 int sort_t_compare(const void *ap,const void *bp){
 135   long a=*(long *)ap;
 136   long b=*(long *)bp;
 137   double val=sort_t_vals[a]-sort_t_vals[b];
 138   if(val<0)return(1);
 139   if(val>0)return(-1);
 140   return(0);
 141 }
 142
 143 void vqgen_iterate(vqgen *v,int biasp){
 144   static int iteration=0;
 145   long i,j;
 146   double avmetric=0.;
 147
 148   FILE *as;
 149   FILE *graph;
 150   FILE *err;
 151   FILE *bias;
 152   char name[80];
 153
 154
 155   sprintf(name,"space%d.m",iteration);
 156   graph=fopen(name,"w");
 157
 158   sprintf(name,"err%d.m",iteration);
 159   err=fopen(name,"w");
 160
 161   sprintf(name,"bias%d.m",iteration);
 162   bias=fopen(name,"w");
 163
 164   sprintf(name,"as%d.m",iteration);
 165   as=fopen(name,"w");
 166
 167
 168   /* init */
 169   memset(v->entry_next,0,sizeof(double)*v->elements*v->entries);
 170   memset(v->assigned,0,sizeof(long)*v->entries);
 171   memset(v->metric,0,sizeof(double)*v->entries);
 172
 173   if(v->entries<2)exit(1);
 174
 175   /* assign all the points, accumulate metric */
 176   for(i=0;i<v->points;i++){
 177     double firstmetric=v->metric_func(v,_now(v,0),_point(v,i))+v->bias[0];
 178     double secondmetric=v->metric_func(v,_now(v,1),_point(v,i))+v->bias[1];
 179     long   firstentry=0;
 180     long   secondentry=0;
 181
 182     if(firstmetric>secondmetric){
 183       double tempmetric=firstmetric;
 184       firstmetric=secondmetric;
 185       secondmetric=tempmetric;
 186       firstentry=1;
 187       secondentry=1;
 188     }
 189
 190     for(j=2;j<v->entries;j++){
 191       double thismetric=v->metric_func(v,_now(v,j),_point(v,i))+v->bias[j];
 192       if(thismetric<secondmetric){
 193         if(thismetric<firstmetric){
 194           secondmetric=firstmetric;
 195           secondentry=firstentry;
 196           firstmetric=thismetric;
 197           firstentry=j;
 198         }else{
 199           secondmetric=thismetric;
 200           secondentry=j;
 201         }
 202       }
 203     }
 204
 205     v->first[i]=firstentry;
 206     v->second[i]=secondentry;
 207     v->metric[firstentry]+=firstmetric-v->bias[firstentry];
 208     v->assigned[firstentry]++;
 209
 210     {
 211       double *no=_now(v,firstentry);
 212       double *ne=_next(v,firstentry);
 213       double *p=_point(v,i);
 214       for(j=0;j<v->elements;j++)
 215         ne[j]+=p[j];
 216       fprintf(graph,"%g %g\n%g %g\n\n",p[0],p[1],no[0],no[1]);
 217     }
 218   }
 219
 220   /* new midpoints */
 221   for(i=0;i<v->entries;i++){
 222     double *next=_next(v,i);
 223     double *now=_now(v,i);
 224     if(v->assigned[i]){
 225       for(j=0;j<v->elements;j++)
 226         next[j]/=v->assigned[i];
 227     }else{
 228       for(j=0;j<v->elements;j++)
 229         next[j]=now[j];
 230     }
 231   }
 232
 233   /* average global metric */
 234   for(i=0;i<v->entries;i++){
 235     avmetric+=v->metric[i];
 236   }
 237   avmetric/=v->entries;
 238   fprintf(stderr,"%ld... ",iteration);
 239
 240   /* positive/negative 'pressure' */
 241   if(biasp){
 242     /* Another round of n*m. Above we had to index by point.  Now, we
 243        have to index by entry */
 244     long   *index=alloca(v->points*sizeof(long));
 245     double *vals =alloca(v->points*sizeof(double));
 246     double *met  =alloca(v->points*sizeof(double));
 247     for(i=0;i<v->entries;i++){
 248
 249       /* sort all n-thousand points by effective metric */
 250
 251       for(j=0;j<v->points;j++){
 252         int threshentry;
 253         index[j]=j;
 254         if(v->first[j]==i){
 255           /* use the secondary entry as the threshhold */
 256           threshentry=v->second[j];
 257         }else{
 258           /* use the primary entry as the threshhold */
 259           threshentry=v->first[j];
 260         }
 261
 262         /* bias value right at the threshhold to grab this point */
 263         met[j]=v->metric_func(v,_now(v,i),_point(v,j));
 264         vals[j]=
 265           v->metric_func(v,_now(v,threshentry),_point(v,j))+
 266           v->bias[threshentry]-
 267           v->metric_func(v,_now(v,i),_point(v,j));
 268       }
 269
 270       sort_t_vals=vals;
 271       qsort(index,v->points,sizeof(long),sort_t_compare);
 272
 273       {
 274         /* find the desired number of points and level of bias.  We
 275            nominally want all entries to represent the same number of
 276            points, but we may also have a metric range restriction */
 277         int desired=(v->points+v->entries-1)/v->entries;
 278
 279         if(v->errspread>=1.){
 280           double acc=0.;
 281           for(j=0;j<desired;j++){
 282             acc+=met[index[j]];
 283             if(acc*v->errspread>avmetric)break;
 284           }
 285         }else
 286           j=desired;
 287         v->bias[i]=vals[index[j-1]];
 288       }
 289     }
 290   }else{
 291     memset(v->bias,0,sizeof(double)*v->entries);
 292   }
 293
 294   /* dump state, report metric */
 295   for(i=0;i<v->entries;i++){
 296     for(j=0;j<v->assigned[i];j++)
 297       fprintf(err,"%g\n",v->metric[i]/v->assigned[i]);
 298     fprintf(bias,"%g\n",v->bias[i]);
 299     fprintf(as,"%ld\n",v->assigned[i]);
 300   }
 301
 302   fprintf(stderr,"average metric: %g\n",avmetric/v->elements/v->points*v->entries);
 303
 304   fclose(as);
 305   fclose(err);
 306   fclose(bias);
 307   fclose(graph);
 308   memcpy(v->entry_now,v->entry_next,sizeof(double)*v->entries*v->elements);
 309   iteration++;
 310 }
 311
 312 int main(int argc,char *argv[]){
 313   FILE *in=fopen(argv[1],"r");
 314   vqgen v;
 315   char buffer[160];
 316   int i;
 317
 318   vqgen_init(&v,4,24,_dist_sq,1.);
 319
 320   while(fgets(buffer,160,in)){
 321     double a[8];
 322     if(sscanf(buffer,"%lf %lf %lf %lf",
 323               a,a+1,a+2,a+3)==4)
 324       vqgen_addpoint(&v,a);
 325   }
 326   fclose(in);
 327
 328   for(i=0;i<20;i++)
 329     vqgen_iterate(&v,0);
 330
 331   for(i=0;i<100;i++)
 332     vqgen_iterate(&v,1);
 333
 334   for(i=0;i<5;i++)
 335     vqgen_iterate(&v,0);
 336
 337   return(0);
 338 }
 339
 340