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  * Vector quantization is a *bit* like an n-body problem in
  54  * m-dimensional space implemented as a monte-carlo simulation.  But
  55  * not really.  We're first-order only, and the 'forces' are
  56  * attractive or repulsive. Think 'xspringies' but without velocity or
  57  * mass. */
  58
  59 /* internal helpers *****************************************************/
  60 double *_point(vqgen *v,long ptr){
  61   return v->pointlist+(v->elements*ptr);
  62 }
  63
  64 double *_now(vqgen *v,long ptr){
  65   return v->entry_now+(v->elements*ptr);
  66 }
  67
  68 double *_next(vqgen *v,long ptr){
  69   return v->entry_next+(v->elements*ptr);
  70 }
  71
  72 double _error_func(vqgen *v,double *a, double *b){
  73   int i;
  74   int el=v->elements;
  75   double acc=0.;
  76   for(i=0;i<el;i++)
  77     acc+=(a[i]-b[i])*(a[i]-b[i]);
  78   return acc;
  79 }
  80
  81 void vqgen_init(vqgen *v,int elements,int entries){
  82   memset(v,0,sizeof(vqgen));
  83
  84   v->elements=elements;
  85   v->allocated=8192;
  86   v->pointlist=malloc(v->allocated*v->elements*sizeof(double));
  87
  88   v->entries=entries;
  89   v->entry_now=malloc(v->entries*v->elements*sizeof(double));
  90   v->entry_next=malloc(v->entries*v->elements*sizeof(double));
  91   v->assigned=malloc(v->entries*sizeof(long));
  92   v->error=malloc(v->entries*sizeof(double));
  93   v->bias=malloc(v->entries*sizeof(double));
  94   {
  95     int i;
  96     for(i=0;i<v->entries;i++)
  97       v->bias[i]=1.;
  98   }
  99
 100   /*v->lasterror=-1;*/
 101
 102   v->error_func=_error_func;
 103 }
 104
 105 void _vqgen_seed(vqgen *v){
 106   memcpy(v->entry_now,v->pointlist,sizeof(double)*v->entries*v->elements);
 107 }
 108
 109 void vqgen_addpoint(vqgen *v, double *p){
 110   if(v->points>=v->allocated){
 111     v->allocated*=2;
 112     v->pointlist=realloc(v->pointlist,v->allocated*v->elements*sizeof(double));
 113   }
 114
 115   memcpy(_point(v,v->points),p,sizeof(double)*v->elements);
 116   v->points++;
 117   if(v->points==v->entries)_vqgen_seed(v);
 118 }
 119
 120 void vqgen_iterate(vqgen *v){
 121   static int iteration=0;
 122   long i,j;
 123   double averror=0.;
 124
 125   FILE *graph;
 126   FILE *err;
 127   FILE *bias;
 128   char name[80];
 129
 130
 131   sprintf(name,"space%d.m",iteration);
 132   graph=fopen(name,"w");
 133
 134   sprintf(name,"err%d.m",iteration);
 135   err=fopen(name,"w");
 136
 137   sprintf(name,"bias%d.m",iteration);
 138   bias=fopen(name,"w");
 139
 140
 141   /* init */
 142   memset(v->entry_next,0,sizeof(double)*v->elements*v->entries);
 143   memset(v->assigned,0,sizeof(long)*v->entries);
 144   memset(v->error,0,sizeof(double)*v->entries);
 145
 146   /* assign all the points, accumulate error */
 147   for(i=0;i<v->points;i++){
 148     double besterror=v->error_func(v,_now(v,0),_point(v,i))*v->bias[0];
 149     long   bestentry=0;
 150     for(j=1;j<v->entries;j++){
 151       double thiserror=v->error_func(v,_now(v,j),_point(v,i))*v->bias[j];
 152       if(thiserror<besterror){
 153         besterror=thiserror;
 154         bestentry=j;
 155       }
 156     }
 157
 158     v->error[bestentry]+=v->error_func(v,_now(v,bestentry),_point(v,i));
 159     v->assigned[bestentry]++;
 160     {
 161       double *n=_next(v,bestentry);
 162       double *p=_point(v,i);
 163       for(j=0;j<v->elements;j++)
 164         n[j]+=p[j];
 165     }
 166
 167     {
 168       double *n=_now(v,bestentry);
 169       double *p=_point(v,i);
 170       fprintf(graph,"%g %g\n%g %g\n\n",p[0],p[1],n[0],n[1]);
 171     }
 172   }
 173
 174   /* new midpoints */
 175   for(i=0;i<v->entries;i++){
 176     double *next=_next(v,i);
 177     double *now=_now(v,i);
 178     if(v->assigned[i]){
 179       for(j=0;j<v->elements;j++)
 180         next[j]/=v->assigned[i];
 181     }else{
 182       for(j=0;j<v->elements;j++)
 183         next[j]=now[j];
 184     }
 185   }
 186
 187   /* average global error */
 188   for(i=0;i<v->entries;i++)
 189     averror+=v->error[i];
 190
 191   averror/=v->entries;
 192
 193   /* positive/negative 'pressure' */
 194
 195   if(iteration%10){
 196     for(i=0;i<v->entries;i++){
 197       double bias=0;
 198       if(v->error[i])
 199         bias=1.+ (averror-v->error[i])/v->error[i];
 200
 201       if(bias>5)bias=5;
 202       if(bias<.02)bias=.2;
 203       v->bias[i]/=bias;
 204     }
 205   }else{
 206     int i;
 207     for(i=0;i<v->entries;i++)
 208       v->bias[i]=1.;
 209   }
 210
 211   /* dump state, report error */
 212   for(i=0;i<v->entries;i++){
 213     fprintf(err,"%g\n",v->error[i]);
 214     fprintf(bias,"%g\n",v->bias[i]);
 215   }
 216
 217   fprintf(stderr,"average error: %g\n",averror);
 218
 219   fclose(err);
 220   fclose(bias);
 221   fclose(graph);
 222   memcpy(v->entry_now,v->entry_next,sizeof(double)*v->entries*v->elements);
 223   iteration++;
 224 }
 225
 226 int main(int argc,char *argv[]){
 227   FILE *in=fopen(argv[1],"r");
 228   vqgen v;
 229   char buffer[80];
 230   int i;
 231
 232   vqgen_init(&v,2,128);
 233
 234   while(fgets(buffer,80,in)){
 235     double a[2];
 236     if(sscanf(buffer,"%lf %lf",a,a+1)==2)
 237       vqgen_addpoint(&v,a);
 238   }
 239   fclose(in);
 240
 241   for(i=0;i<100;i++)
 242     vqgen_iterate(&v);
 243
 244   return(0);
 245 }