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]=0.;
 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
 139   FILE *graph;
 140   FILE *err;
 141   FILE *bias;
 142   char name[80];
 143
 144
 145   sprintf(name,"space%d.m",iteration);
 146   graph=fopen(name,"w");
 147
 148   sprintf(name,"err%d.m",iteration);
 149   err=fopen(name,"w");
 150
 151   sprintf(name,"bias%d.m",iteration);
 152   bias=fopen(name,"w");
 153
 154
 155   /* init */
 156   memset(v->entry_next,0,sizeof(double)*v->elements*v->entries);
 157   memset(v->assigned,0,sizeof(long)*v->entries);
 158   memset(v->error,0,sizeof(double)*v->entries);
 159
 160   /* assign all the points, accumulate error */
 161   for(i=0;i<v->points;i++){
 162     double besterror=v->error_func(v,_now(v,0),_point(v,i))+v->bias[0];
 163     long   bestentry=0;
 164     for(j=1;j<v->entries;j++){
 165       double thiserror=v->error_func(v,_now(v,j),_point(v,i))+v->bias[j];
 166       if(thiserror<besterror){
 167         besterror=thiserror;
 168         bestentry=j;
 169       }
 170     }
 171
 172     v->error[bestentry]+=v->error_func(v,_now(v,bestentry),_point(v,i));
 173     realerror+=sqrt(v->error_func(v,_now(v,bestentry),_point(v,i))/v->elements);
 174
 175     v->assigned[bestentry]++;
 176     {
 177       double *n=_next(v,bestentry);
 178       double *p=_point(v,i);
 179       for(j=0;j<v->elements;j++)
 180         n[j]+=p[j];
 181     }
 182
 183     {
 184       double *n=_now(v,bestentry);
 185       double *p=_point(v,i);
 186       fprintf(graph,"%g %g\n%g %g\n\n",p[0],p[1],n[0],n[1]);
 187     }
 188   }
 189
 190   /* new midpoints */
 191   for(i=0;i<v->entries;i++){
 192     double *next=_next(v,i);
 193     double *now=_now(v,i);
 194     if(v->assigned[i]){
 195       for(j=0;j<v->elements;j++)
 196         next[j]/=v->assigned[i];
 197     }else{
 198       for(j=0;j<v->elements;j++)
 199         next[j]=now[j];
 200     }
 201   }
 202
 203   /* average global error */
 204   for(i=0;i<v->entries;i++){
 205     averror+=v->error[i];
 206     if(v->error[i]==0)fprintf(stderr,"%d ",i);
 207   }
 208   fprintf(stderr,"\n",i);
 209
 210   averror/=v->entries;
 211
 212   /* positive/negative 'pressure' */
 213
 214   if(biasp){
 215     for(i=0;i<v->entries;i++){
 216       double bias=0;
 217       if(v->error[i]){
 218         bias=(averror-v->error[i])/v->assigned[i]*.2;
 219         v->bias[i]-=bias;
 220       }else{
 221         fprintf(stderr,"de-biasing\n");
 222         memset(v->bias,0,sizeof(double)*v->entries);
 223         break;
 224       }
 225
 226       /*if(bias>.1)bias=.1;
 227         if(bias<-.1)bias=-.1;*/
 228     }
 229     fprintf(stderr,"\n");
 230   }else{
 231     fprintf(stderr,"de-biasing\n");
 232     memset(v->bias,0,sizeof(double)*v->entries);
 233   }
 234
 235   /* dump state, report error */
 236   for(i=0;i<v->entries;i++){
 237     fprintf(err,"%g\n",v->error[i]);
 238     fprintf(bias,"%g\n",v->bias[i]);
 239   }
 240
 241   fprintf(stderr,"average error: %g\n",realerror/v->points);
 242
 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%10);
 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