first cut of complete VQ utilities
authorMonty <xiphmont@xiph.org>
Thu, 6 Jan 2000 13:57:15 +0000 (13:57 +0000)
committerMonty <xiphmont@xiph.org>
Thu, 6 Jan 2000 13:57:15 +0000 (13:57 +0000)
Monty

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

vq/Makefile.in
vq/bookutil.c
vq/bookutil.h
vq/cascade.c
vq/metrics.c [new file with mode: 0644]
vq/partition.c [new file with mode: 0644]
vq/vqgen.c
vq/vqgen.h
vq/vqsplit.c

index 5e4ba0e..11da4a8 100644 (file)
@@ -1,4 +1,4 @@
-# $Id: Makefile.in,v 1.6 2000/01/05 10:14:53 xiphmont Exp $
+# $Id: Makefile.in,v 1.7 2000/01/06 13:57:09 xiphmont Exp $
 
 ###############################################################################
 #                                                                             #
@@ -29,7 +29,7 @@ HFILES =      ../include/vorbis/codebook.h vqgen.h vqext.h bookutil.h
 
 OFILES =       vqgen.o vqsplit.o bookutil.o
 ALLOFILES =    $(OFILES) lspdata.o genericdata.o train.o build.o run.o\
-               cascade.o
+               cascade.o partition.o metrics.o
 
 all:
        $(MAKE) target CFLAGS="$(OPT)"
@@ -40,7 +40,7 @@ debug:
 profile: 
        $(MAKE) target CFLAGS="$(PROFILE)"
 
-target:        lspvqtrain genericvqtrain vqbuild vqcascade
+target:        lspvqtrain genericvqtrain vqbuild vqcascade vqpartition vqmetrics
 
 lspvqtrain:    $(OFILES) lspdata.o train.o
        $(CC) $(CFLAGS) $(LDFLAGS) $^ -o $@ $(LIBS)
@@ -54,6 +54,12 @@ vqbuild:     $(OFILES) build.o
 vqcascade:     $(OFILES) run.o cascade.o
        $(CC) $(CFLAGS) $(LDFLAGS) $^ -o $@ $(LIBS)
 
+vqpartition:   $(OFILES) run.o partition.o
+       $(CC) $(CFLAGS) $(LDFLAGS) $^ -o $@ $(LIBS)
+
+vqmetrics:     $(OFILES) run.o metrics.o
+       $(CC) $(CFLAGS) $(LDFLAGS) $^ -o $@ $(LIBS)
+
 $(ALLOFILES):  $(HFILES)
 
 .c.o:
index d4e25af..d69be48 100644 (file)
@@ -12,7 +12,7 @@
  ********************************************************************
 
  function: utility functions for loading .vqh and .vqd files
- last mod: $Id: bookutil.c,v 1.2 2000/01/05 15:04:55 xiphmont Exp $
+ last mod: $Id: bookutil.c,v 1.3 2000/01/06 13:57:10 xiphmont Exp $
 
  ********************************************************************/
 
@@ -67,6 +67,7 @@ char *get_line(FILE *in){
         switch(c){
         case '\n':
         case EOF:
+          linebuffer[sofar]='\0';
           gotline=1;
           break;
         default:
@@ -132,16 +133,20 @@ int get_vector(codebook *b,FILE *in,double *a){
 
   while(1){
 
-    if(get_next_value(in,a))
-      break;
+    if(get_line_value(in,a)){
+      sequence_base=0.;
+      if(get_next_value(in,a))
+       break;
+    }
 
     for(i=1;i<b->dim;i++)
       if(get_line_value(in,a+i))
        break;
     
     if(i==b->dim){
+      double temp=a[b->dim-1];
       for(i=0;i<b->dim;i++)a[i]-=sequence_base;
-      if(b->q_sequencep)sequence_base=a[b->dim-1];
+      if(b->q_sequencep)sequence_base=temp;
       return(0);
     }
     sequence_base=0.;
@@ -345,3 +350,35 @@ double float24_unpack(long val){
   return(ldexp(mant,exp-17-VQ_FEXP_BIAS));
 }
 
+void spinnit(char *s,int n){
+  static int p=0;
+  static long lasttime=0;
+  long test;
+  struct timeval thistime;
+
+  gettimeofday(&thistime,NULL);
+  test=thistime.tv_sec*10+thistime.tv_usec/100000;
+  if(lasttime!=test){
+    lasttime=test;
+
+    fprintf(stderr,"%s%d ",s,n);
+
+    p++;if(p>3)p=0;
+    switch(p){
+    case 0:
+      fprintf(stderr,"|    \r");
+      break;
+    case 1:
+      fprintf(stderr,"/    \r");
+      break;
+    case 2:
+      fprintf(stderr,"-    \r");
+      break;
+    case 3:
+      fprintf(stderr,"\\    \r");
+      break;
+    }
+    fflush(stderr);
+  }
+}
+
index 174fa48..ebddf25 100644 (file)
@@ -12,7 +12,7 @@
  ********************************************************************
 
  function: utility functions for loading .vqh and .vqd files
- last mod: $Id: bookutil.h,v 1.1 2000/01/05 10:14:54 xiphmont Exp $
+ last mod: $Id: bookutil.h,v 1.2 2000/01/06 13:57:11 xiphmont Exp $
 
  ********************************************************************/
 
@@ -20,6 +20,8 @@
 #define _V_BOOKUTIL_H_
 
 #include <stdio.h>
+#include <sys/time.h>
+
 #include "vorbis/codebook.h"
 
 extern void      codebook_unquantize(codebook *b);
@@ -34,7 +36,7 @@ extern char     *find_seek_to(FILE *in,char *s);
 extern codebook *codebook_load(char *filename);
 extern int       codebook_entry(codebook *b,double *val);
 
-
+extern void spinnit(char *s,int n);
 extern long float24_pack(double val);
 extern double float24_unpack(long val);
 
index c7dcc5a..51819d0 100644 (file)
  ********************************************************************
 
  function: function call to do simple data cascading
- last mod: $Id: cascade.c,v 1.2 2000/01/05 15:04:56 xiphmont Exp $
+ last mod: $Id: cascade.c,v 1.3 2000/01/06 13:57:12 xiphmont Exp $
 
  ********************************************************************/
 
-/* this one just outputs to stdout */
+/* this one outputs residue to stdout. */
 
+#include <stdlib.h>
+#include <unistd.h>
+#include <math.h>
 #include "bookutil.h"
 
+/* set up metrics */
+
+double count=0.;
+
 void process_preprocess(codebook *b,char *basename){
 }
 void process_postprocess(codebook *b,char *basename){
+  fprintf(stderr,"Done.                      \n");
 }
 
 void process_vector(codebook *b,double *a){
@@ -33,6 +41,8 @@ void process_vector(codebook *b,double *a){
   for(i=0;i<b->dim;i++)
     fprintf(stdout,"%f, ",a[i]-e[i]);
   fprintf(stdout,"\n");
+
+  if((long)(count++)%100)spinnit("working.... lines: ",count);
 }
 
 void process_usage(void){
diff --git a/vq/metrics.c b/vq/metrics.c
new file mode 100644 (file)
index 0000000..ee632d2
--- /dev/null
@@ -0,0 +1,182 @@
+/********************************************************************
+ *                                                                  *
+ * 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-2000             *
+ * by Monty <monty@xiph.org> and The XIPHOPHORUS Company            *
+ * http://www.xiph.org/                                             *
+ *                                                                  *
+ ********************************************************************
+
+ function: function calls to collect codebook metrics
+ last mod: $Id: metrics.c,v 1.1 2000/01/06 13:57:13 xiphmont Exp $
+
+ ********************************************************************/
+
+
+#include <stdlib.h>
+#include <unistd.h>
+#include <math.h>
+#include "bookutil.h"
+
+/* set up metrics */
+
+double meanamplitude_acc=0.;
+double meanamplitudesq_acc=0.;
+double meanerror_acc=0.;
+double meanerrorsq_acc=0.;
+double meandev_acc=0.;
+
+double *histogram=NULL;
+double *histogram_errorsq=NULL;
+double *histogram_distance=NULL;
+double *histogram_hi=NULL;
+double *histogram_lo=NULL;
+
+double count=0.;
+
+int histerrsort(const void *a, const void *b){
+  double av=histogram_distance[*((long *)a)];
+  double bv=histogram_distance[*((long *)b)];
+  if(av<bv)return(-1);
+  return(1);
+}
+
+void process_preprocess(codebook *b,char *basename){
+  histogram=calloc(b->entries,sizeof(double));
+  histogram_distance=calloc(b->entries,sizeof(double));
+  histogram_errorsq=calloc(b->entries*b->dim,sizeof(double));
+  histogram_hi=calloc(b->entries*b->dim,sizeof(double));
+  histogram_lo=calloc(b->entries*b->dim,sizeof(double));
+}
+
+void process_postprocess(codebook *b,char *basename){
+  int i,j,k;
+  char *buffer=alloca(strlen(basename)+80);
+
+  fprintf(stderr,"Done.  Processed %ld data points for %ld entries:\n",
+         (long)count,b->entries);
+  fprintf(stderr,"\tglobal mean amplitude: %g\n",
+         meanamplitude_acc/(count*b->dim));
+  fprintf(stderr,"\tglobal mean squared amplitude: %g\n",
+         sqrt(meanamplitudesq_acc/(count*b->dim)));
+
+  fprintf(stderr,"\tglobal mean error: %g\n",
+         meanerror_acc/(count*b->dim));
+  fprintf(stderr,"\tglobal mean squared error: %g\n",
+         sqrt(meanerrorsq_acc/(count*b->dim)));
+  fprintf(stderr,"\tglobal mean deviation: %g\n",
+         meandev_acc/(count*b->dim));
+  {
+    FILE *out;
+    sprintf(buffer,"%s-fit.m",basename);
+    out=fopen(buffer,"w");
+    if(!out){
+      fprintf(stderr,"Could not open file %s for writing\n",buffer);
+      exit(1);
+    }
+
+    for(i=0;i<b->entries;i++){
+      for(k=0;k<b->dim;k++){
+       fprintf(out,"%d, %g, %g\n",
+               i*b->dim+k,(b->valuelist+i*b->dim)[k],
+               sqrt((histogram_errorsq+i*b->dim)[k]/histogram[i]));
+      }
+    }
+    fclose(out);
+
+    sprintf(buffer,"%s-worst.m",basename);
+    out=fopen(buffer,"w");
+    if(!out){
+      fprintf(stderr,"Could not open file %s for writing\n",buffer);
+      exit(1);
+    }
+
+    for(i=0;i<b->entries;i++){
+      for(k=0;k<b->dim;k++){
+       fprintf(out,"%d, %g, %g, %g\n",
+               i*b->dim+k,(b->valuelist+i*b->dim)[k],
+               (b->valuelist+i*b->dim)[k]+(histogram_lo+i*b->dim)[k],
+               (b->valuelist+i*b->dim)[k]+(histogram_hi+i*b->dim)[k]);
+      }
+    }
+    fclose(out);
+  }
+
+  {
+    FILE *out;
+    long *index=alloca(sizeof(long)*b->entries);
+    sprintf(buffer,"%s-distance.m",basename);
+    out=fopen(buffer,"w");
+    if(!out){
+      fprintf(stderr,"Could not open file %s for writing\n",buffer);
+      exit(1);
+    }
+    for(j=0;j<b->entries;j++){
+      if(histogram[j])histogram_distance[j]/=histogram[j];
+      index[j]=j;
+    }
+
+    qsort(index,b->entries,sizeof(long),histerrsort);
+
+    for(j=0;j<b->entries;j++)
+      for(k=0;k<histogram[index[j]];k++)
+       fprintf(out,"%g,\n",histogram_distance[index[j]]);
+    fclose(out);
+               
+  }
+
+}
+
+void process_vector(codebook *b,double *a){
+  int entry=codebook_entry(b,a);
+  double *e=b->valuelist+b->dim*entry;
+  int i;
+  double amplitude=0.;
+  double distance=0.;
+  double base=0.;
+  for(i=0;i<b->dim;i++){
+    double error=a[i]-e[i];
+    if(b->q_sequencep){
+      amplitude=a[i]-base;
+      base=a[i];
+    }else
+      amplitude=a[i];
+    
+    meanamplitude_acc+=fabs(amplitude);
+    meanamplitudesq_acc+=amplitude*amplitude;
+    meanerror_acc+=fabs(error);
+    meanerrorsq_acc+=error*error;
+    
+    if(amplitude)
+      meandev_acc+=fabs((a[i]-e[i])/amplitude);
+    else
+      meandev_acc+=fabs(a[i]-e[i]); /* yeah, yeah */
+    
+    histogram_errorsq[entry*b->dim+i]+=error*error;
+    if(histogram[entry]==0 || histogram_hi[entry*b->dim+i]<error)
+      histogram_hi[entry*b->dim+i]=error;
+    if(histogram[entry]==0 || histogram_lo[entry*b->dim+i]>error)
+      histogram_lo[entry*b->dim+i]=error;
+    distance+=error*error;
+  }
+
+  histogram[entry]++;  
+  histogram_distance[entry]+=sqrt(distance);
+
+  if((long)(count++)%100)spinnit("working.... lines: ",count);
+}
+
+void process_usage(void){
+  fprintf(stderr,
+         "usage: vqmetrics <codebook>.vqh datafile.vqd [datafile.vqd]...\n\n"
+         "       data can be taken on stdin.  Output goes to output files:\n"
+         "       basename-fit.m:      gnuplot: mean square error by entry value\n"
+         "       basename-worst.m:    gnuplot: worst error by entry value\n"
+         "       basename-distance.m: gnuplot file showing distance probability\n"
+         "\n");
+
+}
diff --git a/vq/partition.c b/vq/partition.c
new file mode 100644 (file)
index 0000000..52cc819
--- /dev/null
@@ -0,0 +1,70 @@
+/********************************************************************
+ *                                                                  *
+ * 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-2000             *
+ * by Monty <monty@xiph.org> and The XIPHOPHORUS Company            *
+ * http://www.xiph.org/                                             *
+ *                                                                  *
+ ********************************************************************
+
+ function: function call to do data partitioning
+ last mod: $Id: partition.c,v 1.1 2000/01/06 13:57:13 xiphmont Exp $
+
+ ********************************************************************/
+
+/* make n files, one for each cell.  partition the input data into the
+   N cells */
+
+#include <stdlib.h>
+#include <stdio.h>
+#include "bookutil.h"
+
+/* open all the partitioning files */
+FILE **out;
+long count;
+
+void process_preprocess(codebook *b,char *basename){
+  int i;
+  char *buffer=alloca(strlen(basename)+80);
+
+  out=malloc(sizeof(FILE *)*b->entries);
+  for(i=0;i<b->entries;i++){
+    sprintf(buffer,"%s-%dp.vqd",basename,i);
+    out[i]=fopen(buffer,"w");
+    if(out[i]==NULL){
+      fprintf(stderr,"Could not open file %s\n",buffer);
+      exit(1);
+    }
+  }
+}
+
+void process_postprocess(codebook *b,char *basename){
+  int i;
+  for(i=0;i<b->entries;i++)
+    fclose(out[i]);
+  fprintf(stderr,"Done.                      \n");
+}
+
+/* just redirect this entry to the appropriate partition file */
+void process_vector(codebook *b,double *a){
+  int entry=codebook_entry(b,a);
+  int i;
+
+  for(i=0;i<b->dim;i++)
+    fprintf(out[entry],"%f, ",a[i]);
+  fprintf(out[entry],"\n");
+  if(count++%100)spinnit("working.... lines: ",count);
+}
+
+void process_usage(void){
+  fprintf(stderr,
+         "usage: vqpartition <codebook>.vqh datafile.vqd [datafile.vqd]...\n\n"
+         "       data can be taken on stdin.  partitioning data is written to\n"
+  
+         "       files named <basename>-<n>p.vqh.\n\n");
+
+}
index 7101531..610df4b 100644 (file)
@@ -12,7 +12,7 @@
  ********************************************************************
 
  function: train a VQ codebook 
- last mod: $Id: vqgen.c,v 1.27 2000/01/05 15:05:01 xiphmont Exp $
+ last mod: $Id: vqgen.c,v 1.28 2000/01/06 13:57:13 xiphmont Exp $
 
  ********************************************************************/
 
@@ -82,38 +82,6 @@ void _vqgen_seed(vqgen *v){
 
 /* External calls *******************************************************/
 
-void spinnit(char *s,int n){
-  static int p=0;
-  static long lasttime=0;
-  long test;
-  struct timeval thistime;
-
-  gettimeofday(&thistime,NULL);
-  test=thistime.tv_sec*10+thistime.tv_usec/100000;
-  if(lasttime!=test){
-    lasttime=test;
-
-    fprintf(stderr,"%s%d ",s,n);
-
-    p++;if(p>3)p=0;
-    switch(p){
-    case 0:
-      fprintf(stderr,"|    \r");
-      break;
-    case 1:
-      fprintf(stderr,"/    \r");
-      break;
-    case 2:
-      fprintf(stderr,"-    \r");
-      break;
-    case 3:
-      fprintf(stderr,"\\    \r");
-      break;
-    }
-    fflush(stderr);
-  }
-}
-
 /* We have two forms of quantization; in the first, each vector
    element in the codebook entry is orthogonal.  Residues would use this
    quantization for example.
@@ -272,7 +240,7 @@ double vqgen_iterate(vqgen *v){
     exit(1);
   }
 
-  /* fill in nearest points for entries */
+  /* fill in nearest points for entry biasing */
   /*memset(v->bias,0,sizeof(double)*v->entries);*/
   memset(nearcount,0,sizeof(long)*v->entries);
   memset(v->assigned,0,sizeof(long)*v->entries);
@@ -283,7 +251,7 @@ double vqgen_iterate(vqgen *v){
     long   firstentry=0;
     long   secondentry=1;
 
-    spinnit("working... ",v->points-i);
+    if(i%100)spinnit("biasing... ",v->points+v->points+v->entries-i);
 
     if(firstmetric>secondmetric){
       double temp=firstmetric;
@@ -309,23 +277,6 @@ double vqgen_iterate(vqgen *v){
     }
 
     j=firstentry;
-      
-#ifdef NOISY
-    fprintf(cells,"%g %g\n%g %g\n\n",
-          _now(v,j)[0],_now(v,j)[1],
-          ppt[0],ppt[1]);
-#endif
-
-    j=firstentry;
-    meterror+=firstmetric-v->bias[firstentry];
-    /* set up midpoints for next iter */
-    if(v->assigned[j]++)
-      for(k=0;k<v->elements;k++)
-       vN(new,j)[k]+=ppt[k];
-    else
-      for(k=0;k<v->elements;k++)
-       vN(new,j)[k]=ppt[k];
-
     for(j=0;j<v->entries;j++){
       
       double thismetric;
@@ -346,13 +297,16 @@ double vqgen_iterate(vqgen *v){
       if(k<desired){
        nearbiasptr[k]=thismetric;
        k++;
-       if(k==desired)
+       if(k==desired){
+         spinnit("biasing... ",v->points+v->points+v->entries-i);
          qsort(nearbiasptr,desired,sizeof(double),directdsort);
+       }
        
       }else if(thismetric>nearbiasptr[desired-1]){
        nearbiasptr[k]=thismetric;
        k++;
        if(k==desired2){
+         spinnit("biasing... ",v->points+v->points+v->entries-i);
          qsort(nearbiasptr,desired2,sizeof(double),directdsort);
          k=desired;
        }
@@ -365,7 +319,7 @@ double vqgen_iterate(vqgen *v){
   for(i=0;i<v->entries;i++){
     double *nearbiasptr=nearbias+desired2*i;
 
-    spinnit("working... ",v->entries-i);
+    spinnit("biasing... ",v->points+v->entries-i);
 
     /* due to the delayed sorting, we likely need to finish it off....*/
     if(nearcount[i]>desired)
@@ -374,6 +328,41 @@ double vqgen_iterate(vqgen *v){
     v->bias[i]=nearbiasptr[desired-1];
   }
 
+  /* Now assign with new bias and find new midpoints */
+  for(i=0;i<v->points;i++){
+    double *ppt=v->weight_func(v,_point(v,i));
+    double firstmetric=v->metric_func(v,_now(v,0),ppt)+v->bias[0];
+    long   firstentry=0;
+
+    if(i%100)spinnit("centering... ",v->points-i);
+
+    for(j=0;j<v->entries;j++){
+      double thismetric=v->metric_func(v,_now(v,j),ppt)+v->bias[j];
+      if(thismetric<firstmetric){
+       firstmetric=thismetric;
+       firstentry=j;
+      }
+    }
+
+    j=firstentry;
+      
+#ifdef NOISY
+    fprintf(cells,"%g %g\n%g %g\n\n",
+          _now(v,j)[0],_now(v,j)[1],
+          ppt[0],ppt[1]);
+#endif
+
+    meterror+=firstmetric-v->bias[firstentry];
+    /* set up midpoints for next iter */
+    if(v->assigned[j]++)
+      for(k=0;k<v->elements;k++)
+       vN(new,j)[k]+=ppt[k];
+    else
+      for(k=0;k<v->elements;k++)
+       vN(new,j)[k]=ppt[k];
+
+  }
+
   /* assign midpoints */
 
   for(j=0;j<v->entries;j++){
index 7cbd86d..8d5a1b4 100644 (file)
  ********************************************************************
 
  function: build a VQ codebook 
- last mod: $Id: vqgen.h,v 1.9 2000/01/05 10:14:58 xiphmont Exp $
+ last mod: $Id: vqgen.h,v 1.10 2000/01/06 13:57:14 xiphmont Exp $
 
  ********************************************************************/
 
 #ifndef _VQGEN_H_
 #define _VQGEN_H_
 
-#include <sys/time.h>
-
 typedef struct vqgen{
   int it;
   int elements;
@@ -69,11 +67,6 @@ extern double vqgen_iterate(vqgen *v);
 extern void vqgen_unquantize(vqgen *v,quant_meta *q);
 extern void vqgen_quantize(vqgen *v,quant_meta *q);
 
-extern void spinnit(char *s,int n);
-
-extern long float24_pack(double val);
-extern double float24_unpack(long val);
-
 #endif
 
 
index 5da1ea5..49bd82e 100644 (file)
@@ -12,7 +12,7 @@
  ********************************************************************
 
  function: build a VQ codebook and the encoding decision 'tree'
- last mod: $Id: vqsplit.c,v 1.10 2000/01/05 15:05:02 xiphmont Exp $
+ last mod: $Id: vqsplit.c,v 1.11 2000/01/06 13:57:15 xiphmont Exp $
 
  ********************************************************************/
 
@@ -183,7 +183,7 @@ int lp_split(vqgen *v,codebook *b,
   /* more than one way to do this part.  For small sets, we can brute
      force it. */
 
-  if(entries<8 || (double)points*entries*entries<128.*1024*1024){
+  if(entries<8 || (double)points*entries*entries<16.*1024*1024){
     /* try every pair possibility */
     double best=0;
     double this;