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 11 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
  48   /* point cache */
  49   double *pointlist;
  50   long   points;
  51   long   allocated;
  52
  53   /* entries */
  54   double *entry_now;
  55   double *entry_next;
  56   long   *assigned;
  57   double *error;
  58   double *bias;
  59   long   entries;
  60
  61   double (*error_func)   (struct vqgen *v,double *a,double *b);
  62   double (*bias_func)     (struct vqgen *v,int entry);
  63 } vqgen;
  64
  65 /* internal helpers *****************************************************/
  66 double *_point(vqgen *v,long ptr){
  67   return v->pointlist+(v->elements*ptr);
  68 }
  69
  70 double *_now(vqgen *v,long ptr){
  71   return v->entry_now+(v->elements*ptr);
  72 }
  73
  74 double *_next(vqgen *v,long ptr){
  75   return v->entry_next+(v->elements*ptr);
  76 }
  77
  78 double _error_func(vqgen *v,double *a, double *b){
  79   int i;
  80   int el=v->elements;
  81   double acc=0.;
  82   for(i=0;i<el;i++)
  83     acc+=(a[i]-b[i])*(a[i]-b[i]);
  84   return acc;
  85 }
  86
  87 double _vqgen_distance(vqgen *v,double *a, double *b){
  88   int i;
  89   double acc=0.;
  90   for(i=0;i<v->elements;i++)acc+=(a[i]-b[i])*(a[i]-b[i]);
  91   return sqrt(acc);
  92 }
  93
  94 void vqgen_init(vqgen *v,int elements,int entries){
  95   memset(v,0,sizeof(vqgen));
  96
  97   v->elements=elements;
  98   v->allocated=8192;
  99   v->pointlist=malloc(v->allocated*v->elements*sizeof(double));
 100
 101   v->entries=entries;
 102   v->entry_now=malloc(v->entries*v->elements*sizeof(double));
 103   v->entry_next=malloc(v->entries*v->elements*sizeof(double));
 104   v->assigned=malloc(v->entries*sizeof(long));
 105   v->error=malloc(v->entries*sizeof(double));
 106   v->bias=malloc(v->entries*sizeof(double));
 107   {
 108     int i;
 109     for(i=0;i<v->entries;i++)
 110       v->bias[i]=10.;
 111   }
 112
 113   /*v->lasterror=-1;*/
 114
 115   v->error_func=_error_func;
 116 }
 117
 118 void _vqgen_seed(vqgen *v){
 119   memcpy(v->entry_now,v->pointlist,sizeof(double)*v->entries*v->elements);
 120 }
 121
 122 void vqgen_addpoint(vqgen *v, double *p){
 123   if(v->points>=v->allocated){
 124     v->allocated*=2;
 125     v->pointlist=realloc(v->pointlist,v->allocated*v->elements*sizeof(double));
 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 void vqgen_iterate(vqgen *v,int biasp){
 134   static int iteration=0;
 135   long i,j;
 136   double averror=0.;
 137   double realerror=0.;
 138   double *distance=alloca(sizeof(double)*v->entries);
 139
 140   FILE *graph;
 141   FILE *err;
 142   FILE *bias;
 143   FILE *assi;
 144   char name[80];
 145
 146   memset(distance,0,sizeof(double)*v->entries);
 147
 148   sprintf(name,"space%d.m",iteration);
 149   graph=fopen(name,"w");
 150
 151   sprintf(name,"err%d.m",iteration);
 152   err=fopen(name,"w");
 153
 154   sprintf(name,"bias%d.m",iteration);
 155   bias=fopen(name,"w");
 156
 157   sprintf(name,"as%d.m",iteration);
 158   assi=fopen(name,"w");
 159
 160
 161   /* init */
 162   memset(v->entry_next,0,sizeof(double)*v->elements*v->entries);
 163   memset(v->assigned,0,sizeof(long)*v->entries);
 164   memset(v->error,0,sizeof(double)*v->entries);
 165
 166   /* assign all the points, accumulate error */
 167   for(i=0;i<v->points;i++){
 168     double besterror=v->error_func(v,_now(v,0),_point(v,i))+v->bias[0];
 169     long   bestentry=0;
 170     for(j=1;j<v->entries;j++){
 171       double thiserror=v->error_func(v,_now(v,j),_point(v,i))+v->bias[j];
 172       if(thiserror<besterror){
 173         besterror=thiserror;
 174         bestentry=j;
 175       }
 176     }
 177
 178     v->error[bestentry]+=v->error_func(v,_now(v,bestentry),_point(v,i));
 179     distance[bestentry]+=_vqgen_distance(v,_now(v,bestentry),_point(v,i));
 180     realerror+=sqrt(v->error_func(v,_now(v,bestentry),_point(v,i))/v->elements);
 181     v->assigned[bestentry]++;
 182     {
 183       double *n=_next(v,bestentry);
 184       double *p=_point(v,i);
 185       for(j=0;j<v->elements;j++)
 186         n[j]+=p[j];
 187     }
 188
 189     {
 190       double *n=_now(v,bestentry);
 191       double *p=_point(v,i);
 192       fprintf(graph,"%g %g\n%g %g\n\n",p[0],p[1],n[0],n[1]);
 193     }
 194   }
 195
 196   /* new midpoints */
 197   for(i=0;i<v->entries;i++){
 198     double *next=_next(v,i);
 199     double *now=_now(v,i);
 200     if(v->assigned[i]){
 201       for(j=0;j<v->elements;j++)
 202         next[j]/=v->assigned[i];
 203     }else{
 204       for(j=0;j<v->elements;j++)
 205         next[j]=now[j];
 206     }
 207   }
 208
 209   /* average global error */
 210   for(i=0;i<v->entries;i++){
 211     averror+=v->error[i];
 212     if(v->error[i]==0)fprintf(stderr,"%d ",i);
 213   }
 214   fprintf(stderr,"\n");
 215
 216   averror/=v->entries;
 217
 218   /* positive/negative 'pressure' */
 219
 220   if(biasp){
 221     for(i=0;i<v->entries;i++){
 222       double ratio=(float)v->assigned[i]*v->entries/v->points;
 223       double averr=v->error[i]/v->assigned[i];
 224       ratio=(ratio-1.)*averr*.1+1.;
 225       v->bias[i]*=ratio;
 226     }
 227   }else{
 228     int i;
 229     fprintf(stderr,"de-biasing\n");
 230     for(i=0;i<v->entries;i++)v->bias[i]=10.;
 231   }
 232
 233   /* dump state, report error */
 234   for(i=0;i<v->entries;i++){
 235     fprintf(err,"%g\n",v->error[i]);
 236     fprintf(assi,"%ld\n",v->assigned[i]);
 237     fprintf(bias,"%g\n",v->bias[i]);
 238   }
 239
 240   fprintf(stderr,"average error: %g\n",realerror/v->points);
 241
 242   fclose(assi);
 243   fclose(err);
 244   fclose(bias);
 245   fclose(graph);
 246   memcpy(v->entry_now,v->entry_next,sizeof(double)*v->entries*v->elements);
 247   iteration++;
 248 }
 249
 250 int main(int argc,char *argv[]){
 251   FILE *in=fopen(argv[1],"r");
 252   vqgen v;
 253   char buffer[160];
 254   int i;
 255
 256   vqgen_init(&v,4,128);
 257
 258   while(fgets(buffer,160,in)){
 259     double a[8];
 260     if(sscanf(buffer,"%lf %lf %lf %lf",
 261               a,a+1,a+2,a+3)==4)
 262       vqgen_addpoint(&v,a);
 263   }
 264   fclose(in);
 265
 266     vqgen_iterate(&v,0);
 267     vqgen_iterate(&v,0);
 268     vqgen_iterate(&v,0);
 269
 270   for(i=0;i<100;i++)
 271     vqgen_iterate(&v,i%20);
 272
 273     vqgen_iterate(&v,0);
 274     vqgen_iterate(&v,0);
 275     vqgen_iterate(&v,0);
 276     vqgen_iterate(&v,0);
 277     vqgen_iterate(&v,0);
 278     vqgen_iterate(&v,0);
 279     vqgen_iterate(&v,0);
 280     vqgen_iterate(&v,0);
 281     vqgen_iterate(&v,0);
 282     vqgen_iterate(&v,0);
 283     vqgen_iterate(&v,0);
 284     vqgen_iterate(&v,0);
 285     vqgen_iterate(&v,0);
 286
 287   return(0);
 288 }
 289
 290