Enough of the VQ generation works to no longer want to lose it :-)

Monty

svn path=/trunk/vorbis/; revision=166

diff --git a/vq/vqgen.c b/vq/vqgen.c
new file mode 100644 (file)
index 0000000..884603e
--- /dev/null
+++ b/vq/vqgen.c
@@ -0,0 +1,245 @@
+/********************************************************************
+ *                                                                  *
+ * THIS FILE IS PART OF THE Ogg Vorbis SOFTWARE CODEC SOURCE CODE.  *
+ * USE, DISTRIBUTION AND REPRODUCTION OF THIS SOURCE IS GOVERNED BY *
+ * THE GNU PUBLIC LICENSE 2, WHICH IS INCLUDED WITH THIS SOURCE.    *
+ * PLEASE READ THESE TERMS DISTRIBUTING.                            *
+ *                                                                  *
+ * THE OggSQUISH SOURCE CODE IS (C) COPYRIGHT 1994-1999             *
+ * by 1999 Monty <monty@xiph.org> and The XIPHOPHORUS Company       *
+ * http://www.xiph.org/                                             *
+ *                                                                  *
+ ********************************************************************
+
+ function: build a VQ codebook 
+ author: Monty <xiphmont@mit.edu>
+ modifications by: Monty
+ last modification date: Nov 09 1999
+
+ ********************************************************************/
+
+/* This code is *not* part of libvorbis.  It is used to generate
+   trained codebooks offline and then spit the results into a
+   pregenerated codebook that is compiled into libvorbis.  It is an
+   expensive (but good) algorithm.  Run it on big iron. */
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <string.h>
+
+
+typedef struct vqgen{
+  int    elements;
+
+  /* point cache */
+  double *pointlist; 
+  long   points;
+  long   allocated;
+
+  /* entries */
+  double *entry_now;
+  double *entry_next;
+  long   *assigned;
+  double *error;
+  double *bias;
+  long   entries;
+
+  double (*error_func)   (struct vqgen *v,double *a,double *b);
+} vqgen;
+
+
+/*************************************************************************
+ * Vector quantization is a *bit* like an n-body problem in
+ * m-dimensional space implemented as a monte-carlo simulation.  But
+ * not really.  We're first-order only, and the 'forces' are
+ * attractive or repulsive. Think 'xspringies' but without velocity or
+ * mass. */
+
+/* internal helpers *****************************************************/
+double *_point(vqgen *v,long ptr){
+  return v->pointlist+(v->elements*ptr);
+}
+
+double *_now(vqgen *v,long ptr){
+  return v->entry_now+(v->elements*ptr);
+}
+
+double *_next(vqgen *v,long ptr){
+  return v->entry_next+(v->elements*ptr);
+}
+
+double _error_func(vqgen *v,double *a, double *b){
+  int i;
+  int el=v->elements;
+  double acc=0.;
+  for(i=0;i<el;i++)
+    acc+=(a[i]-b[i])*(a[i]-b[i]);
+  return acc;
+}
+
+void vqgen_init(vqgen *v,int elements,int entries){
+  memset(v,0,sizeof(vqgen));
+
+  v->elements=elements;
+  v->allocated=8192;
+  v->pointlist=malloc(v->allocated*v->elements*sizeof(double));
+
+  v->entries=entries;
+  v->entry_now=malloc(v->entries*v->elements*sizeof(double));
+  v->entry_next=malloc(v->entries*v->elements*sizeof(double));
+  v->assigned=malloc(v->entries*sizeof(long));
+  v->error=malloc(v->entries*sizeof(double));
+  v->bias=malloc(v->entries*sizeof(double));
+  {
+    int i;
+    for(i=0;i<v->entries;i++)
+      v->bias[i]=1.;
+  }
+  
+  /*v->lasterror=-1;*/
+
+  v->error_func=_error_func;
+}
+
+void _vqgen_seed(vqgen *v){
+  memcpy(v->entry_now,v->pointlist,sizeof(double)*v->entries*v->elements);
+}
+
+void vqgen_addpoint(vqgen *v, double *p){
+  if(v->points>=v->allocated){
+    v->allocated*=2;
+    v->pointlist=realloc(v->pointlist,v->allocated*v->elements*sizeof(double));
+  }
+  
+  memcpy(_point(v,v->points),p,sizeof(double)*v->elements);
+  v->points++;
+  if(v->points==v->entries)_vqgen_seed(v);
+}
+
+void vqgen_iterate(vqgen *v){
+  static int iteration=0;
+  long i,j;
+  double averror=0.;
+
+  FILE *graph;
+  FILE *err;
+  FILE *bias;
+  char name[80];
+  
+
+  sprintf(name,"space%d.m",iteration);
+  graph=fopen(name,"w");
+
+  sprintf(name,"err%d.m",iteration);
+  err=fopen(name,"w");
+
+  sprintf(name,"bias%d.m",iteration);
+  bias=fopen(name,"w");
+
+
+  /* init */
+  memset(v->entry_next,0,sizeof(double)*v->elements*v->entries);
+  memset(v->assigned,0,sizeof(long)*v->entries);
+  memset(v->error,0,sizeof(double)*v->entries);
+
+  /* assign all the points, accumulate error */
+  for(i=0;i<v->points;i++){
+    double besterror=v->error_func(v,_now(v,0),_point(v,i))*v->bias[0];
+    long   bestentry=0;
+    for(j=1;j<v->entries;j++){
+      double thiserror=v->error_func(v,_now(v,j),_point(v,i))*v->bias[j];
+      if(thiserror<besterror){
+       besterror=thiserror;
+       bestentry=j;
+      }
+    }
+    
+    v->error[bestentry]+=v->error_func(v,_now(v,bestentry),_point(v,i));
+    v->assigned[bestentry]++;
+    {
+      double *n=_next(v,bestentry);
+      double *p=_point(v,i);
+      for(j=0;j<v->elements;j++)
+       n[j]+=p[j];
+    }
+
+    {
+      double *n=_now(v,bestentry);
+      double *p=_point(v,i);
+      fprintf(graph,"%g %g\n%g %g\n\n",p[0],p[1],n[0],n[1]);
+    }
+  }
+
+  /* new midpoints */
+  for(i=0;i<v->entries;i++){
+    double *next=_next(v,i);
+    double *now=_now(v,i);
+    if(v->assigned[i]){
+      for(j=0;j<v->elements;j++)
+       next[j]/=v->assigned[i];
+    }else{
+      for(j=0;j<v->elements;j++)
+       next[j]=now[j];
+    }
+  }
+  
+  /* average global error */
+  for(i=0;i<v->entries;i++)
+    averror+=v->error[i];
+  
+  averror/=v->entries;
+
+  /* positive/negative 'pressure' */
+  
+  if(iteration%10){
+    for(i=0;i<v->entries;i++){
+      double bias=0;
+      if(v->error[i])
+       bias=1.+ (averror-v->error[i])/v->error[i];
+      
+      if(bias>5)bias=5;
+      if(bias<.02)bias=.2;
+      v->bias[i]/=bias;
+    }
+  }else{
+    int i;
+    for(i=0;i<v->entries;i++)
+      v->bias[i]=1.;
+  }
+
+  /* dump state, report error */
+  for(i=0;i<v->entries;i++){
+    fprintf(err,"%g\n",v->error[i]);
+    fprintf(bias,"%g\n",v->bias[i]);
+  }
+
+  fprintf(stderr,"average error: %g\n",averror);
+
+  fclose(err);
+  fclose(bias);
+  fclose(graph);
+  memcpy(v->entry_now,v->entry_next,sizeof(double)*v->entries*v->elements);
+  iteration++;
+}
+
+int main(int argc,char *argv[]){
+  FILE *in=fopen(argv[1],"r");
+  vqgen v;
+  char buffer[80];
+  int i;
+
+  vqgen_init(&v,2,128);
+
+  while(fgets(buffer,80,in)){
+    double a[2];
+    if(sscanf(buffer,"%lf %lf",a,a+1)==2)
+      vqgen_addpoint(&v,a);
+  }
+  fclose(in);
+
+  for(i=0;i<100;i++)
+    vqgen_iterate(&v);
+
+  return(0);
+}