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: Dec 10 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 /* There are so many optimizations to explore in *both* stages that
  27    considering the undertaking is almost withering.  For now, we brute
  28    force it all */
  29
  30 #include <stdlib.h>
  31 #include <stdio.h>
  32 #include <math.h>
  33 #include <string.h>
  34 #include "vqgen.h"
  35
  36 /* Codebook generation happens in two steps:
  37
  38    1) Train the codebook with data collected from the encoder: We use
  39    one of a few error metrics (which represent the distance between a
  40    given data point and a candidate point in the training set) to
  41    divide the training set up into cells representing roughly equal
  42    probability of occurring.
  43
  44    2) Generate the codebook and auxiliary data from the trained data set
  45 */
  46
  47 /* Codebook training ****************************************************
  48  *
  49  * The basic idea here is that a VQ codebook is like an m-dimensional
  50  * foam with n bubbles.  The bubbles compete for space/volume and are
  51  * 'pressurized' [biased] according to some metric.  The basic alg
  52  * iterates through allowing the bubbles to compete for space until
  53  * they converge (if the damping is dome properly) on a steady-state
  54  * solution. Individual input points, collected from libvorbis, are
  55  * used to train the algorithm monte-carlo style.  */
  56
  57 /* internal helpers *****************************************************/
  58 #define vN(data,i) (data+v->elements*i)
  59
  60 /* default metric; squared 'distance' from desired value. */
  61 double _dist_sq(vqgen *v,double *a, double *b){
  62   int i;
  63   int el=v->elements;
  64   double acc=0.;
  65   for(i=0;i<el;i++){
  66     double val=(a[i]-b[i]);
  67     acc+=val*val;
  68   }
  69   return acc;
  70 }
  71
  72 /* *must* be beefed up. */
  73 void _vqgen_seed(vqgen *v){
  74   memcpy(v->entrylist,v->pointlist,sizeof(double)*v->entries*v->elements);
  75 }
  76
  77 /* External calls *******************************************************/
  78
  79 void vqgen_init(vqgen *v,int elements,int entries,
  80                 double (*metric)(vqgen *,double *, double *)){
  81   memset(v,0,sizeof(vqgen));
  82
  83   v->elements=elements;
  84   v->allocated=32768;
  85   v->pointlist=malloc(v->allocated*v->elements*sizeof(double));
  86
  87   v->entries=entries;
  88   v->entrylist=malloc(v->entries*v->elements*sizeof(double));
  89   v->assigned=malloc(v->entries*sizeof(long));
  90   v->bias=calloc(v->entries,sizeof(double));
  91   if(metric)
  92     v->metric_func=metric;
  93   else
  94     v->metric_func=_dist_sq;
  95 }
  96
  97 void vqgen_addpoint(vqgen *v, double *p){
  98   if(v->points>=v->allocated){
  99     v->allocated*=2;
 100     v->pointlist=realloc(v->pointlist,v->allocated*v->elements*sizeof(double));
 101   }
 102
 103   memcpy(_point(v,v->points),p,sizeof(double)*v->elements);
 104   v->points++;
 105   if(v->points==v->entries)_vqgen_seed(v);
 106 }
 107
 108 double vqgen_iterate(vqgen *v){
 109   long   i,j,k;
 110   double fdesired=(double)v->points/v->entries;
 111   long  desired=fdesired;
 112   double asserror=0.;
 113   double meterror=0.;
 114   double *new=malloc(sizeof(double)*v->entries*v->elements);
 115   long   *nearcount=malloc(v->entries*sizeof(long));
 116   double *nearbias=malloc(v->entries*desired*sizeof(double));
 117
 118 #ifdef NOISY
 119   char buff[80];
 120   FILE *assig;
 121   FILE *bias;
 122   FILE *cells;
 123   sprintf(buff,"cells%d.m",v->it);
 124   cells=fopen(buff,"w");
 125   sprintf(buff,"assig%d.m",v->it);
 126   assig=fopen(buff,"w");
 127   sprintf(buff,"bias%d.m",v->it);
 128   bias=fopen(buff,"w");
 129 #endif
 130
 131   fprintf(stderr,"Pass #%d... ",v->it);
 132
 133   if(v->entries<2){
 134     fprintf(stderr,"generation requires at least two entries\n");
 135     exit(1);
 136   }
 137
 138   /* fill in nearest points for entries */
 139   /*memset(v->bias,0,sizeof(double)*v->entries);*/
 140   memset(nearcount,0,sizeof(long)*v->entries);
 141   memset(v->assigned,0,sizeof(long)*v->entries);
 142   for(i=0;i<v->points;i++){
 143     double *ppt=_point(v,i);
 144     double firstmetric=v->metric_func(v,_now(v,0),ppt)+v->bias[0];
 145     double secondmetric=v->metric_func(v,_now(v,1),ppt)+v->bias[1];
 146     long   firstentry=0;
 147     long   secondentry=1;
 148     if(firstmetric>secondmetric){
 149       double temp=firstmetric;
 150       firstmetric=secondmetric;
 151       secondmetric=temp;
 152       firstentry=1;
 153       secondentry=0;
 154     }
 155
 156     for(j=2;j<v->entries;j++){
 157       double thismetric=v->metric_func(v,_now(v,j),_point(v,i))+v->bias[j];
 158       if(thismetric<secondmetric){
 159         if(thismetric<firstmetric){
 160           secondmetric=firstmetric;
 161           secondentry=firstentry;
 162           firstmetric=thismetric;
 163           firstentry=j;
 164         }else{
 165           secondmetric=thismetric;
 166           secondentry=j;
 167         }
 168       }
 169     }
 170
 171     j=firstentry;
 172     meterror+=sqrt(_dist_sq(v,_now(v,j),_point(v,i)));
 173     /* set up midpoints for next iter */
 174     if(v->assigned[j]++)
 175       for(k=0;k<v->elements;k++)
 176         vN(new,j)[k]+=_point(v,i)[k];
 177     else
 178       for(k=0;k<v->elements;k++)
 179         vN(new,j)[k]=_point(v,i)[k];
 180
 181
 182 #ifdef NOISY
 183     fprintf(cells,"%g %g\n%g %g\n\n",
 184             _now(v,j)[0],_now(v,j)[1],
 185             _point(v,i)[0],_point(v,i)[1]);
 186 #endif
 187
 188     for(j=0;j<v->entries;j++){
 189
 190       double thismetric;
 191       double *nearbiasptr=nearbias+desired*j;
 192       long k=nearcount[j]-1;
 193
 194       /* 'thismetric' is to be the bias value necessary in the current
 195          arrangement for entry j to capture point i */
 196       if(firstentry==j){
 197         /* use the secondary entry as the threshhold */
 198         thismetric=secondmetric-v->metric_func(v,_now(v,j),_point(v,i));
 199       }else{
 200         /* use the primary entry as the threshhold */
 201         thismetric=firstmetric-v->metric_func(v,_now(v,j),_point(v,i));
 202       }
 203
 204       if(k>=0 && thismetric>nearbiasptr[k]){
 205
 206         /* start at the end and search backward for where this entry
 207            belongs */
 208
 209         for(;k>0;k--) if(nearbiasptr[k-1]>=thismetric)break;
 210
 211         /* insert at k.  Shift and inject. */
 212         memmove(nearbiasptr+k+1,nearbiasptr+k,(desired-k-1)*sizeof(double));
 213         nearbiasptr[k]=thismetric;
 214
 215         if(nearcount[j]<desired)nearcount[j]++;
 216
 217       }else{
 218         if(nearcount[j]<desired){
 219           /* we checked the thresh earlier.  We know this is the
 220              last entry */
 221           nearbiasptr[nearcount[j]++]=thismetric;
 222         }
 223       }
 224     }
 225   }
 226
 227   /* inflate/deflate */
 228   for(i=0;i<v->entries;i++)
 229     v->bias[i]=nearbias[(i+1)*desired-1];
 230
 231   /* assign midpoints */
 232
 233   for(j=0;j<v->entries;j++){
 234     asserror+=fabs(v->assigned[j]-fdesired);
 235     if(v->assigned[j])
 236       for(k=0;k<v->elements;k++)
 237         _now(v,j)[k]=vN(new,j)[k]/v->assigned[j];
 238 #ifdef NOISY
 239     fprintf(assig,"%ld\n",v->assigned[j]);
 240     fprintf(bias,"%g\n",v->bias[j]);
 241 #endif
 242   }
 243
 244   asserror/=(v->entries*fdesired);
 245   fprintf(stderr,": dist %g(%g) metric error=%g \n",
 246           asserror,fdesired,meterror/v->points/v->elements);
 247   v->it++;
 248
 249   free(new);
 250   free(nearcount);
 251   free(nearbias);
 252 #ifdef NOISY
 253   fclose(assig);
 254   fclose(bias);
 255   fclose(cells);
 256 #endif
 257   return(asserror);
 258 }
 259