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 09 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 typedef struct vqgen{
  33   int    elements;
  34
  35   /* point cache */
  36   double *pointlist;
  37   long   points;
  38   long   allocated;
  39
  40   /* entries */
  41   double *entry_now;
  42   double *entry_next;
  43   long   *assigned;
  44   double *error;
  45   double *bias;
  46   long   entries;
  47
  48   double (*error_func)   (struct vqgen *v,double *a,double *b);
  49 } vqgen;
  50
  51
  52 /*************************************************************************
  53  * Generating a vector quantization codebook is a *bit* like
  54  * simulating a weird m-dimensional foam, building a number of bubbles
  55  * with a [supposedly] constant, minimum error.  We train the 'foam'
  56  * with data points like a monte-carlo simulation. */
  57
  58 /* internal helpers *****************************************************/
  59 double *_point(vqgen *v,long ptr){
  60   return v->pointlist+(v->elements*ptr);
  61 }
  62
  63 double *_now(vqgen *v,long ptr){
  64   return v->entry_now+(v->elements*ptr);
  65 }
  66
  67 double *_next(vqgen *v,long ptr){
  68   return v->entry_next+(v->elements*ptr);
  69 }
  70
  71 double _error_func(vqgen *v,double *a, double *b){
  72   int i;
  73   int el=v->elements;
  74   double acc=0.;
  75   for(i=0;i<el;i++)
  76     acc+=(a[i]-b[i])*(a[i]-b[i]);
  77   return acc;
  78 }
  79
  80 void vqgen_init(vqgen *v,int elements,int entries){
  81   memset(v,0,sizeof(vqgen));
  82
  83   v->elements=elements;
  84   v->allocated=8192;
  85   v->pointlist=malloc(v->allocated*v->elements*sizeof(double));
  86
  87   v->entries=entries;
  88   v->entry_now=malloc(v->entries*v->elements*sizeof(double));
  89   v->entry_next=malloc(v->entries*v->elements*sizeof(double));
  90   v->assigned=malloc(v->entries*sizeof(long));
  91   v->error=malloc(v->entries*sizeof(double));
  92   v->bias=malloc(v->entries*sizeof(double));
  93   {
  94     int i;
  95     for(i=0;i<v->entries;i++)
  96       v->bias[i]=0.;
  97   }
  98
  99   /*v->lasterror=-1;*/
 100
 101   v->error_func=_error_func;
 102 }
 103
 104 void _vqgen_seed(vqgen *v){
 105   memcpy(v->entry_now,v->pointlist,sizeof(double)*v->entries*v->elements);
 106 }
 107
 108 void vqgen_addpoint(vqgen *v, double *p){
 109   if(v->points>=v->allocated){
 110     v->allocated*=2;
 111     v->pointlist=realloc(v->pointlist,v->allocated*v->elements*sizeof(double));
 112   }
 113
 114   memcpy(_point(v,v->points),p,sizeof(double)*v->elements);
 115   v->points++;
 116   if(v->points==v->entries)_vqgen_seed(v);
 117 }
 118
 119 void vqgen_iterate(vqgen *v,int biasp){
 120   static int iteration=0;
 121   long i,j;
 122   double averror=0.;
 123
 124   FILE *graph;
 125   FILE *err;
 126   FILE *bias;
 127   char name[80];
 128
 129
 130   sprintf(name,"space%d.m",iteration);
 131   graph=fopen(name,"w");
 132
 133   sprintf(name,"err%d.m",iteration);
 134   err=fopen(name,"w");
 135
 136   sprintf(name,"bias%d.m",iteration);
 137   bias=fopen(name,"w");
 138
 139
 140   /* init */
 141   memset(v->entry_next,0,sizeof(double)*v->elements*v->entries);
 142   memset(v->assigned,0,sizeof(long)*v->entries);
 143   memset(v->error,0,sizeof(double)*v->entries);
 144
 145   /* assign all the points, accumulate error */
 146   for(i=0;i<v->points;i++){
 147     double besterror=v->error_func(v,_now(v,0),_point(v,i))+v->bias[0];
 148     long   bestentry=0;
 149     for(j=1;j<v->entries;j++){
 150       double thiserror=v->error_func(v,_now(v,j),_point(v,i))+v->bias[j];
 151       if(thiserror<besterror){
 152         besterror=thiserror;
 153         bestentry=j;
 154       }
 155     }
 156
 157     v->error[bestentry]+=v->error_func(v,_now(v,bestentry),_point(v,i));
 158     v->assigned[bestentry]++;
 159     {
 160       double *n=_next(v,bestentry);
 161       double *p=_point(v,i);
 162       for(j=0;j<v->elements;j++)
 163         n[j]+=p[j];
 164     }
 165
 166     {
 167       double *n=_now(v,bestentry);
 168       double *p=_point(v,i);
 169       fprintf(graph,"%g %g\n%g %g\n\n",p[0],p[1],n[0],n[1]);
 170     }
 171   }
 172
 173   /* new midpoints */
 174   for(i=0;i<v->entries;i++){
 175     double *next=_next(v,i);
 176     double *now=_now(v,i);
 177     if(v->assigned[i]){
 178       for(j=0;j<v->elements;j++)
 179         next[j]/=v->assigned[i];
 180     }else{
 181       for(j=0;j<v->elements;j++)
 182         next[j]=now[j];
 183     }
 184   }
 185
 186   /* average global error */
 187   for(i=0;i<v->entries;i++)
 188     averror+=v->error[i];
 189
 190   averror/=v->entries;
 191
 192   /* positive/negative 'pressure' */
 193
 194   if(biasp){
 195     for(i=0;i<v->entries;i++){
 196       double bias=0;
 197       if(v->error[i]){
 198         bias=(averror-v->error[i])/v->assigned[i]*.05;
 199         v->bias[i]-=bias;
 200       }else{
 201         fprintf(stderr,"de-biasing\n");
 202         memset(v->bias,0,sizeof(double)*v->entries);
 203         break;
 204       }
 205
 206       /*if(bias>.1)bias=.1;
 207         if(bias<-.1)bias=-.1;*/
 208     }
 209     fprintf(stderr,"\n");
 210   }else{
 211     fprintf(stderr,"de-biasing\n");
 212     memset(v->bias,0,sizeof(double)*v->entries);
 213   }
 214
 215   /* dump state, report error */
 216   for(i=0;i<v->entries;i++){
 217     fprintf(err,"%g\n",v->error[i]);
 218     fprintf(bias,"%g\n",v->bias[i]);
 219   }
 220
 221   fprintf(stderr,"average error: %g\n",averror);
 222
 223   fclose(err);
 224   fclose(bias);
 225   fclose(graph);
 226   memcpy(v->entry_now,v->entry_next,sizeof(double)*v->entries*v->elements);
 227   iteration++;
 228 }
 229
 230 int main(int argc,char *argv[]){
 231   FILE *in=fopen(argv[1],"r");
 232   vqgen v;
 233   char buffer[80];
 234   int i;
 235
 236   vqgen_init(&v,2,16);
 237
 238   while(fgets(buffer,80,in)){
 239     double a[2];
 240     if(sscanf(buffer,"%lf %lf",a,a+1)==2)
 241       vqgen_addpoint(&v,a);
 242   }
 243   fclose(in);
 244
 245     vqgen_iterate(&v,0);
 246     vqgen_iterate(&v,0);
 247     vqgen_iterate(&v,0);
 248
 249   for(i=0;i<100;i++)
 250     vqgen_iterate(&v,1);
 251
 252     vqgen_iterate(&v,0);
 253     vqgen_iterate(&v,0);
 254     vqgen_iterate(&v,0);
 255     vqgen_iterate(&v,0);
 256     vqgen_iterate(&v,0);
 257     vqgen_iterate(&v,0);
 258     vqgen_iterate(&v,0);
 259     vqgen_iterate(&v,0);
 260     vqgen_iterate(&v,0);
 261     vqgen_iterate(&v,0);
 262     vqgen_iterate(&v,0);
 263     vqgen_iterate(&v,0);
 264     vqgen_iterate(&v,0);
 265
 266   return(0);
 267 }
 268
 269