From: Monty <xiphmont@xiph.org>
Date: Thu, 30 Dec 1999 01:23:36 +0000 (+0000)
Subject: Incremental update
X-Git-Tag: v1.3.3~1306
X-Git-Url: http://review.tizen.org/git/?a=commitdiff_plain;h=f637dd7b71ce97b2e637b717b4d309e8c55f6b93;p=platform%2Fupstream%2Flibvorbis.git

Incremental update

svn path=/trunk/vorbis/; revision=216
---

diff --git a/vq/build.c b/vq/build.c
index 63b74cd..d23c9aa 100644
--- a/vq/build.c
+++ b/vq/build.c
@@ -209,54 +209,68 @@ int main(int argc,char *argv[]){
   for(j=0;j<entries;j++){
     fprintf(out,"\t");
     for(k=0;k<dim;k++)
-      fprintf(out,"%ld, ",b.quantlist[i++]);
+      fprintf(out,"%5ld, ",b.quantlist[i++]);
     fprintf(out,"\n");
   }
   fprintf(out,"};\n\n");
 
   /* codelist */
   fprintf(out,"static long _vq_codelist_%s[] = {\n",name);
-  for(j=0;j<entries;j++){
-    fprintf(out,"\t%ld,\n",b.codelist[j]);
+  for(j=0;j<entries;){
+    fprintf(out,"\t");
+    for(k=0;k<8 && j<entries;k++,j++)
+      fprintf(out,"%ld,",b.codelist[j]);
+    fprintf(out,"\n");
   }
   fprintf(out,"};\n\n");
 
   /* lengthlist */
   fprintf(out,"static long _vq_lengthlist_%s[] = {\n",name);
-  for(j=0;j<entries;j++){
-    fprintf(out,"\t%ld,\n",b.lengthlist[j]);
+  for(j=0;j<entries;){
+    fprintf(out,"\t");
+    for(k=0;k<16 && j<entries;k++,j++)
+      fprintf(out,"%2ld,",b.lengthlist[j]);
+    fprintf(out,"\n");
   }
   fprintf(out,"};\n\n");
 
   /* ptr0 */
   fprintf(out,"static long _vq_ptr0_%s[] = {\n",name);
-  for(j=0;j<b.aux;j++){
-    fprintf(out,"\t%ld,\n",b.ptr0[j]);
+  for(j=0;j<b.aux;){
+    fprintf(out,"\t");
+    for(k=0;k<8 && j<b.aux;k++,j++)
+      fprintf(out,"%5ld,",b.ptr0[j]);
+    fprintf(out,"\n");
   }
   fprintf(out,"};\n\n");
 
   /* ptr1 */
   fprintf(out,"static long _vq_ptr1_%s[] = {\n",name);
-  for(j=0;j<b.aux;j++){
-    fprintf(out,"\t%ld,\n",b.ptr1[j]);
+  for(j=0;j<b.aux;){
+    fprintf(out,"\t");
+    for(k=0;k<8 && j<b.aux;k++,j++)
+      fprintf(out,"%6ld,",b.ptr1[j]);
+    fprintf(out,"\n");
   }
   fprintf(out,"};\n\n");
 
-  /* n */
-  fprintf(out,"static double _vq_n_%s[] = {\n",name);
-  i=0;
-  for(j=0;j<b.aux;j++){
+  /* p */
+  fprintf(out,"static long _vq_p_%s[] = {\n",name);
+  for(j=0;j<b.aux;){
     fprintf(out,"\t");
-    for(k=0;k<dim;k++)
-      fprintf(out,"%g,",b.n[i++]);
+    for(k=0;k<8 && j<b.aux;k++,j++)
+      fprintf(out,"%6ld,",b.p[j]);
     fprintf(out,"\n");
   }
   fprintf(out,"};\n\n");
 
-  /* c */
-  fprintf(out,"static double _vq_c_%s[] = {\n",name);
-  for(j=0;j<b.aux;j++){
-    fprintf(out,"\t%g,\n",b.c[j]);
+  /* q */
+  fprintf(out,"static long _vq_q_%s[] = {\n",name);
+  for(j=0;j<b.aux;){
+    fprintf(out,"\t");
+    for(k=0;k<8 && j<b.aux;k++,j++)
+      fprintf(out,"%6ld,",b.q[j]);
+    fprintf(out,"\n");
   }
   fprintf(out,"};\n\n");
 
@@ -270,8 +284,10 @@ int main(int argc,char *argv[]){
   fprintf(out,"\t_vq_lengthlist_%s,\n",name);
   fprintf(out,"\t_vq_ptr0_%s,\n",name);
   fprintf(out,"\t_vq_ptr1_%s,\n",name);
-  fprintf(out,"\t_vq_n_%s,\n",name);
-  fprintf(out,"\t_vq_c_%s,\n",name);
+  fprintf(out,"\t0,\n",name);
+  fprintf(out,"\t0,\n",name);
+  fprintf(out,"\t_vq_p_%s,\n",name);
+  fprintf(out,"\t_vq_q_%s,\n",name);
   fprintf(out,"\t%ld, %ld };\n\n",b.aux,b.aux);
 
   fprintf(out,"\n#endif\n");
diff --git a/vq/run.c b/vq/run.c
new file mode 100644
index 0000000..2c33edb
--- /dev/null
+++ b/vq/run.c
@@ -0,0 +1,266 @@
+#include "lsp16.vqh"
+#define CODEBOOK _vq_book_lsp16
+
+/********************************************************************
+ *                                                                  *
+ * 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: utility main for loading/testing/running finished codebooks
+ author: Monty <xiphmont@mit.edu>
+ modifications by: Monty
+ last modification date: Dec 28 1999
+
+ ********************************************************************/
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <string.h>
+#include <errno.h>
+
+/* this is a bit silly; it's a C stub used by a Perl script to build a
+   quick executable against the chosen codebook and run tests */
+
+static char *linebuffer=NULL;
+static int  lbufsize=0;
+static char *rline(FILE *in){
+  long sofar=0;
+  if(feof(in))return NULL;
+
+  while(1){
+    int gotline=0;
+
+    while(!gotline){
+      if(sofar>=lbufsize){
+	if(!lbufsize){	
+	  lbufsize=1024;
+	  linebuffer=malloc(lbufsize);
+	}else{
+	  lbufsize*=2;
+	  linebuffer=realloc(linebuffer,lbufsize);
+	}
+      }
+      {
+	long c=fgetc(in);
+	switch(c){
+	case '\n':
+	case EOF:
+	  gotline=1;
+	  break;
+	default:
+	  linebuffer[sofar++]=c;
+	  linebuffer[sofar]='\0';
+	  break;
+	}
+      }
+    }
+    
+    if(linebuffer[0]=='#'){
+      sofar=0;
+    }else{
+      return(linebuffer);
+    }
+  }
+}
+
+void vqbook_unquantize(vqbook *b){
+  long j,k;
+  double mindel=float24_unpack(b->min);
+  double delta=float24_unpack(b->delta);
+  if(!b->valuelist)b->valuelist=malloc(sizeof(double)*b->entries*b->dim);
+  
+  for(j=0;j<b->entries;j++){
+    double last=0.;
+    for(k=0;k<b->dim;k++){
+      double val=b->quantlist[j*b->dim+k]*delta+last+mindel;
+      b->valuelist[j*b->dim+k]=val;
+      if(b->sequencep)last=val;
+
+    }
+  }
+}
+
+
+/* command line:
+   run outbase [-m] [-s <start>,<n>] datafile [-s <start>,<n>] [datafile...]
+
+   produces: outbase-residue.m (error between input data and chosen codewords;
+                                can be used to cascade)
+             outbase-cells.m   (2d gnuplot file of cells if -m)
+
+   currently assumes a 'sequenced' codebook wants its test data to be
+   normalized to begin at zero... */
+
+int main(int argc,char *argv[]){
+  vqbook *b=&CODEBOOK;
+  FILE *cells=NULL;
+  FILE *residue=NULL;
+  FILE *in=NULL;
+  char *line,*name;
+  long i,j,k;
+  int start=0,num=-1;
+
+  double mean=0.,meansquare=0.,mean_count=0.;
+  
+  argv++;
+
+  if(*argv==NULL){
+    fprintf(stderr,"Need a basename.\n");
+    exit(1);
+  }
+
+  name=strdup(*argv);
+  vqbook_unquantize(b);
+
+  {
+    char buf[80];
+    snprintf(buf,80,"%s-residue.vqd",name);
+    residue=fopen(buf,"w");
+    if(!residue){
+      fprintf(stderr,"Unable to open output file %s\n",buf);
+      exit(1);
+    }
+  }   
+
+  /* parse the command line; handle things as the come */
+  argv=argv++;
+  while(*argv){
+    if(argv[0][0]=='-'){
+      /* it's an option */
+      switch(argv[0][1]){
+      case 'm':
+	{
+	  char buf[80];
+	  snprintf(buf,80,"%s-cells.m",name);
+	  cells=fopen(buf,"w");
+	  if(!cells){
+	    fprintf(stderr,"Unable to open output file %s\n",buf);
+	    exit(1);
+	  }
+	}   
+	argv++;
+	break;
+      case 's':
+	if(!argv[1]){
+	  fprintf(stderr,"Option %s missing argument.\n",argv[0]);
+	  exit(1);
+	}      
+	if(sscanf(argv[1],"%d,%d",&start,&num)!=2){
+          num= -1;
+          if(sscanf(argv[1],"%d",&start)!=1){
+	    fprintf(stderr,"Option %s syntax error.\n",argv[0]);
+            exit(1);
+	  }
+        }
+	argv+=2;
+        break;
+      }
+    }else{
+      /* it's an input file */
+      char *file=strdup(*argv++);
+      FILE *in;
+      int cols=-1;
+
+      in=fopen(file,"r");
+      if(in==NULL){
+        fprintf(stderr,"Could not open input file %s\n",file);
+        exit(1);
+      }
+
+      while((line=rline(in))){
+        if(cols==-1){
+          char *temp=line;
+          while(*temp==' ')temp++;
+          for(cols=0;*temp;cols++){
+            while(*temp>32)temp++;
+            while(*temp==' ')temp++;
+          }
+        }
+        {
+          int i;
+          double *p=malloc(cols*sizeof(double));
+          if(start*num+b->dim>cols){
+            fprintf(stderr,"ran out of columns reading %s\n",file);
+            exit(1);
+          }
+          while(*line==' ')line++;
+          for(i=0;i<cols;i++){
+
+            /* static length buffer bug workaround */
+            char *temp=line;
+            char old;
+            while(*temp>32)temp++;
+
+            old=temp[0];
+            temp[0]='\0';
+            p[i]=atof(line);
+            temp[0]=old;
+            
+            while(*line>32)line++;
+            while(*line==' ')line++;
+          }
+          if(num<=0)num=(cols-start)/b->dim;
+          for(i=num-1;i>=0;i--){
+	    int entry;
+	    double *base=p+start+i*b->dim;
+
+	    /* if this is a sequenced book and i||start,
+	       normalize the beginning to zero */
+	    if(b->sequencep && (i>0 || start>0)){
+	      for(k=0;k<b->dim;k++)
+		base[k]-= *(base-1);
+	    }
+
+	    /* assign the point */
+	    entry=vqenc_entry(b,base);
+
+	    /* accumulate metrics */
+	    for(k=0;k<b->dim;k++){
+	      double err=base[k]-b->valuelist[k+entry*b->dim];
+	      mean+=fabs(err);
+	      meansquare+=err*err;
+	      mean_count++;
+	    }
+
+	    /* brute force it... did that work better? */
+
+
+	    /* paint the cell if -m */
+	    if(cells){
+	      fprintf(cells,"%g %g\n%g %g\n\n",
+		      base[0],base[1],
+		      b->valuelist[0+entry*b->dim],
+		      b->valuelist[1+entry*b->dim]);
+	    }
+
+	    /* write cascading data */
+
+
+
+	  }
+          free(b);
+        }
+      }
+      fclose(in);
+    }
+  }
+  
+  
+  fclose(residue);
+  if(cells)fclose(cells);
+
+  /* print accumulated error statistics */
+  fprintf(stderr,"results:\n\tmean squared error:%g\n\tmean error:%g\n\n",
+	  sqrt(meansquare/mean_count),mean/mean_count);
+
+  return 0;
+}
diff --git a/vq/vqgen.h b/vq/vqgen.h
index f351c24..bbc942e 100644
--- a/vq/vqgen.h
+++ b/vq/vqgen.h
@@ -57,6 +57,8 @@ typedef struct vqbook{
   /* auxiliary encoding information. Not used in decode */
   double *n;         /* decision hyperplanes: sum(x_i*n_i)[0<=i<dim]=c */ 
   double *c;
+  long   *p;
+  long   *q;
   long   aux;
   long   alloc;
 
diff --git a/vq/vqsplit.c b/vq/vqsplit.c
index 3a83aec..99a891d 100644
--- a/vq/vqsplit.c
+++ b/vq/vqsplit.c
@@ -101,6 +101,8 @@ void vqsp_count(vqgen *v,long *membership,
 
   /* The entry splitting isn't total, so that storage has to be
      allocated for recursion.  Reuse the entryA/entryB vectors */
+  /* keep the entries in ascending order (relative to the original
+     list); we rely on that stability when ordering p/q choice */
   for(j=0;j<entries;j++){
     if(entryA[j] && entryB[j])C++;
     if(entryA[j])entryA[A++]=entryindex[j];
@@ -121,15 +123,6 @@ void pq_in_out(vqgen *v,double *n,double *c,double *p,double *q){
   }
 }
 
-void pq_center_out(vqgen *v,double *n,double *c,double *center,double *q){
-  int k;
-  *c=0.;
-  for(k=0;k<v->elements;k++){
-    n[k]=(center[k]-q[k])*2.;
-    *c+=center[k]*n[k];
-  }
-}
-
 static void spinnit(void){
   static int p=0;
   static long lasttime=0;
@@ -172,9 +165,7 @@ int lp_split(vqgen *v,vqbook *b,
      even a midpoint division won't disturb the basic properties) */
 
   long ret;
-  double *p;
-  double *q;
-  double *n;
+  double *n=alloca(sizeof(double)*v->elements);
   double c;
   long *entryA=calloc(entries,sizeof(long));
   long *entryB=calloc(entries,sizeof(long));
@@ -185,6 +176,9 @@ int lp_split(vqgen *v,vqbook *b,
   long i,j,k;
 
   long *membership=malloc(sizeof(long)*points);
+
+  long   besti=-1;
+  long   bestj=-1;
   
   /* which cells do points belong to?  Do this before n^2 best pair chooser. */
 
@@ -206,11 +200,6 @@ int lp_split(vqgen *v,vqbook *b,
     membership[i]=firstentry;
   }
 
-  p=alloca(sizeof(double)*v->elements);
-  q=alloca(sizeof(double)*v->elements);
-  n=alloca(sizeof(double)*v->elements);
-  memset(p,0,sizeof(double)*v->elements);
-
   /* We need to find the dividing hyperplane. find the median of each
      axis as the centerpoint and the normal facing farthest point */
 
@@ -220,8 +209,6 @@ int lp_split(vqgen *v,vqbook *b,
   if(entries<8 || points*entries*entries<128*1024*1024){
     /* try every pair possibility */
     double best=0;
-    long   besti=0;
-    long   bestj=0;
     double this;
     for(i=0;i<entries-1;i++){
       for(j=i+1;j<entries;j++){
@@ -235,20 +222,30 @@ int lp_split(vqgen *v,vqbook *b,
 		   &entriesA,&entriesB,&entriesC);
 	this=(entriesA-entriesC)*(entriesB-entriesC);
 
-	if(this>best){
+	/* when choosing best, we also want some form of stability to
+           make sure more branches are pared later; secondary
+           weighting isn;t needed as the entry lists are in ascending
+           order, and we always try p/q in the same sequence */
+	
+	if( (besti==-1) ||
+	    (this>best) ){
+	  
 	  best=this;
-	  besti=i;
-	  bestj=j;
+	  besti=entryindex[i];
+	  bestj=entryindex[j];
+
 	}
       }
     }
-    pq_in_out(v,n,&c,_now(v,entryindex[besti]),_now(v,entryindex[bestj]));
   }else{
+    double *p=alloca(v->elements*sizeof(double));
+    double *q=alloca(v->elements*sizeof(double));
     double best=0.;
-    long   bestj=0;
     
     /* try COG/normal and furthest pairs */
-    /* medianpoint */
+    /* meanpoint */
+    /* eventually, we want to select the closest entry and figure n/c
+       from p/q (because storing n/c is too large */
     for(k=0;k<v->elements;k++){
       spinnit();
       
@@ -258,14 +255,32 @@ int lp_split(vqgen *v,vqbook *b,
       p[k]/=entries;
 
     }
-    
-    /* try every normal */
-    for(j=0;j<entries;j++){
-      double *ppj=_now(v,entryindex[j]);
+
+    /* we go through the entries one by one, looking for the entry on
+       the other side closest to the point of reflection through the
+       center */
+
+    for(i=0;i<entries;i++){
+      double *ppi=_now(v,entryindex[i]);
+      double ref_best=0.;
+      double ref_j=-1;
       double this;
       spinnit();
+      
+      for(k=0;k<v->elements;k++)
+	q[k]=2*p[k]-ppi[k];
+
+      for(j=0;j<entries;j++){
+	if(j!=i){
+	  double this=_dist_sq(v,q,_now(v,entryindex[j]));
+	  if(ref_j==-1 || this<ref_best){
+	    ref_best=this;
+	    ref_j=entryindex[j];
+	  }
+	}
+      }
 
-      pq_center_out(v,n,&c,p,_now(v,entryindex[j]));
+      pq_in_out(v,n,&c,ppi,_now(v,ref_j));
       vqsp_count(v,membership,
 		 entryindex,entries, 
 		 pointindex,points,
@@ -273,19 +288,33 @@ int lp_split(vqgen *v,vqbook *b,
 		 n, c, 
 		 &entriesA,&entriesB,&entriesC);
       this=(entriesA-entriesC)*(entriesB-entriesC);
-      
-      if(this>best){
+
+	/* when choosing best, we also want some form of stability to
+           make sure more branches are pared later; secondary
+           weighting isn;t needed as the entry lists are in ascending
+           order, and we always try p/q in the same sequence */
+	
+      if( (besti==-1) ||
+	  (this>best) ){
+	
 	best=this;
-	bestj=j;
+	besti=entryindex[i];
+	bestj=ref_j;
+
       }
     }
-    
-    pq_center_out(v,n,&c,p,_now(v,entryindex[bestj]));
+    if(besti>bestj){
+      long temp=besti;
+      besti=bestj;
+      bestj=temp;
+    }
+
   }
   
   /* find cells enclosing points */
   /* count A/B points */
 
+  pq_in_out(v,n,&c,_now(v,besti),_now(v,bestj));
   vqsp_count(v,membership,
 	     entryindex,entries, 
 	     pointindex,points,
@@ -332,13 +361,13 @@ int lp_split(vqgen *v,vqbook *b,
       b->alloc*=2;
       b->ptr0=realloc(b->ptr0,sizeof(long)*b->alloc);
       b->ptr1=realloc(b->ptr1,sizeof(long)*b->alloc);
-      b->n=realloc(b->n,sizeof(double)*b->dim*b->alloc);
-      b->c=realloc(b->c,sizeof(double)*b->alloc);
+      b->p=realloc(b->p,sizeof(long)*b->alloc);
+      b->q=realloc(b->q,sizeof(long)*b->alloc);
     }
     
-    memcpy(b->n+b->dim*thisaux,n,sizeof(double)*v->elements);
-    b->c[thisaux]=c;
-
+    b->p[thisaux]=besti;
+    b->q[thisaux]=bestj;
+    
     if(entriesA==1){
       ret=1;
       b->ptr0[thisaux]=entryA[0];
@@ -360,6 +389,27 @@ int lp_split(vqgen *v,vqbook *b,
   return(ret);
 }
 
+static int _node_eq(vqbook *v, long a, long b){
+  long    Aptr0=v->ptr0[a];
+  long    Aptr1=v->ptr1[a];
+  long    Ap   =v->p[a];
+  long    Aq   =v->q[a];
+  long    Bptr0=v->ptr0[b];
+  long    Bptr1=v->ptr1[b];
+  long    Bp   =v->p[b];
+  long    Bq   =v->q[b];
+  int i;
+
+  /* the possibility of choosing the same p and q, but switched, can;t
+     happen because we always look for the best p/q in the same search
+     order and the search is stable */
+
+  if(Aptr0==Bptr0 && Aptr1==Bptr1 && Ap==Bp && Aq==Bq)
+    return(1);
+
+  return(0);
+}
+
 void vqsp_book(vqgen *v, vqbook *b, long *quantlist){
   long *entryindex=malloc(sizeof(long)*v->entries);
   long *pointindex=malloc(sizeof(long)*v->points);
@@ -376,8 +426,8 @@ void vqsp_book(vqgen *v, vqbook *b, long *quantlist){
 
   b->ptr0=malloc(sizeof(long)*b->alloc);
   b->ptr1=malloc(sizeof(long)*b->alloc);
-  b->n=malloc(sizeof(double)*b->dim*b->alloc);
-  b->c=malloc(sizeof(double)*b->alloc);
+  b->p=malloc(sizeof(long)*b->alloc);
+  b->q=malloc(sizeof(long)*b->alloc);
   b->lengthlist=calloc(b->entries,sizeof(long));
   
   /* which cells do points belong to?  Only do this once. */
@@ -407,12 +457,69 @@ void vqsp_book(vqgen *v, vqbook *b, long *quantlist){
   free(entryindex);
   free(pointindex);
 
+  fprintf(stderr,"Paring and rerouting redundant branches:\n");
+  /* The tree is likely big and redundant.  Pare and reroute branches */
+  {
+    int changedflag=1;
+    long *unique_entries=alloca(b->aux*sizeof(long));
+
+    while(changedflag){
+      int nodes=0;
+      changedflag=0;
+
+      fprintf(stderr,"\t...");
+      /* span the tree node by node; list unique decision nodes and
+	 short circuit redundant branches */
+
+      for(i=0;i<b->aux;){
+	int k;
+
+	/* check list of unique decisions */
+	for(j=0;j<nodes;j++)
+	  if(_node_eq(b,i,unique_entries[j]))break;
+
+	if(j==nodes){
+	  /* a new entry */
+	  unique_entries[nodes++]=i++;
+	}else{
+	  /* a redundant entry; find all higher nodes referencing it and
+             short circuit them to the previously noted unique entry */
+	  changedflag=1;
+	  for(k=0;k<b->aux;k++){
+	    if(b->ptr0[k]==-i)b->ptr0[k]=-unique_entries[j];
+	    if(b->ptr1[k]==-i)b->ptr1[k]=-unique_entries[j];
+	  }
+
+	  /* Now, we need to fill in the hole from this redundant
+             entry in the listing.  Insert the last entry in the list.
+             Fix the forward pointers to that last entry */
+	  b->aux--;
+	  b->ptr0[i]=b->ptr0[b->aux];
+	  b->ptr1[i]=b->ptr1[b->aux];
+	  b->p[i]=b->p[b->aux];
+	  b->q[i]=b->q[b->aux];
+	  for(k=0;k<b->aux;k++){
+	    if(b->ptr0[k]==-b->aux)b->ptr0[k]=-i;
+	    if(b->ptr1[k]==-b->aux)b->ptr1[k]=-i;
+	  }
+	  /* hole plugged */
+
+	}
+      }
+
+      fprintf(stderr," %ld remaining\n",b->aux);
+    }
+  }
+
   /* run all training points through the decision tree to get a final
      probability count */
   {
-    long *probability=calloc(b->entries,sizeof(long));
+    long *probability=malloc(b->entries*sizeof(long));
     long *membership=malloc(b->entries*sizeof(long));
 
+    for(i=0;i<b->entries;i++)probability[i]=1; /* trivial guard */
+
+    b->valuelist=v->entrylist; /* temporary for vqenc_entry; replaced later */
     for(i=0;i<v->points;i++){
       int ret=vqenc_entry(b,v->pointlist+i*v->elements);
       probability[ret]++;
@@ -464,10 +571,6 @@ void vqsp_book(vqgen *v, vqbook *b, long *quantlist){
       probability[ret]++;
     }
     
-    /* print the results */
-    for(i=0;i<v->entries;i++)
-      fprintf(stderr,"%ld: count=%ld codelength=%ld\n",i,probability[i],b->lengthlist[i]);
-
     free(probability);
     free(membership);
   }
@@ -489,6 +592,8 @@ void vqsp_book(vqgen *v, vqbook *b, long *quantlist){
     for(i=0;i<b->aux;i++){
       if(b->ptr0[i]>=0)b->ptr0[i]=revindex[b->ptr0[i]];
       if(b->ptr1[i]>=0)b->ptr1[i]=revindex[b->ptr1[i]];
+      b->p[i]=revindex[b->p[i]];
+      b->q[i]=revindex[b->q[i]];
     }
     free(revindex);
 
@@ -519,8 +624,6 @@ void vqsp_book(vqgen *v, vqbook *b, long *quantlist){
 	length=b->lengthlist[i];
       }
       b->codelist[i]=current;
-      fprintf(stderr,"codeword %ld: %ld (length %d)\n",
-	      i, current, b->lengthlist[i]);
       current++;
     }
   }
@@ -535,16 +638,27 @@ void vqsp_book(vqgen *v, vqbook *b, long *quantlist){
       }
     }
   }
-
+  fprintf(stderr,"Done.\n\n");
 }
 
+/* slow version for use here that does not use a preexpanded n/c. */
 int vqenc_entry(vqbook *b,double *val){
   int ptr=0,k;
+  double *n=alloca(b->dim*sizeof(double));
+
   while(1){
-    double c= -b->c[ptr];
-    double *nptr=b->n+b->dim*ptr;
+    double c=0.;
+    double *p=b->valuelist+b->p[ptr]*b->dim;
+    double *q=b->valuelist+b->q[ptr]*b->dim;
+    
+    for(k=0;k<b->dim;k++){
+      n[k]=p[k]-q[k];
+      c-=(p[k]+q[k])*n[k];
+    }
+    c/=2.;
+
     for(k=0;k<b->dim;k++)
-      c+=nptr[k]*val[k];
+      c+=n[k]*val[k];
     if(c>0.) /* in A */
       ptr= -b->ptr0[ptr];
     else     /* in B */