cannam@86: /********************************************************************
cannam@86:  *                                                                  *
cannam@86:  * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE.   *
cannam@86:  * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS     *
cannam@86:  * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE *
cannam@86:  * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING.       *
cannam@86:  *                                                                  *
cannam@86:  * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2001             *
cannam@86:  * by the Xiph.Org Foundation http://www.xiph.org/                  *
cannam@86:  *                                                                  *
cannam@86:  ********************************************************************
cannam@86: 
cannam@86:  function: train a VQ codebook 
cannam@86:  last mod: $Id: vqgen.c 16037 2009-05-26 21:10:58Z xiphmont $
cannam@86: 
cannam@86:  ********************************************************************/
cannam@86: 
cannam@86: /* This code is *not* part of libvorbis.  It is used to generate
cannam@86:    trained codebooks offline and then spit the results into a
cannam@86:    pregenerated codebook that is compiled into libvorbis.  It is an
cannam@86:    expensive (but good) algorithm.  Run it on big iron. */
cannam@86: 
cannam@86: /* There are so many optimizations to explore in *both* stages that
cannam@86:    considering the undertaking is almost withering.  For now, we brute
cannam@86:    force it all */
cannam@86: 
cannam@86: #include <stdlib.h>
cannam@86: #include <stdio.h>
cannam@86: #include <math.h>
cannam@86: #include <string.h>
cannam@86: 
cannam@86: #include "vqgen.h"
cannam@86: #include "bookutil.h"
cannam@86: 
cannam@86: /* Codebook generation happens in two steps: 
cannam@86: 
cannam@86:    1) Train the codebook with data collected from the encoder: We use
cannam@86:    one of a few error metrics (which represent the distance between a
cannam@86:    given data point and a candidate point in the training set) to
cannam@86:    divide the training set up into cells representing roughly equal
cannam@86:    probability of occurring. 
cannam@86: 
cannam@86:    2) Generate the codebook and auxiliary data from the trained data set
cannam@86: */
cannam@86: 
cannam@86: /* Codebook training ****************************************************
cannam@86:  *
cannam@86:  * The basic idea here is that a VQ codebook is like an m-dimensional
cannam@86:  * foam with n bubbles.  The bubbles compete for space/volume and are
cannam@86:  * 'pressurized' [biased] according to some metric.  The basic alg
cannam@86:  * iterates through allowing the bubbles to compete for space until
cannam@86:  * they converge (if the damping is dome properly) on a steady-state
cannam@86:  * solution. Individual input points, collected from libvorbis, are
cannam@86:  * used to train the algorithm monte-carlo style.  */
cannam@86: 
cannam@86: /* internal helpers *****************************************************/
cannam@86: #define vN(data,i) (data+v->elements*i)
cannam@86: 
cannam@86: /* default metric; squared 'distance' from desired value. */
cannam@86: float _dist(vqgen *v,float *a, float *b){
cannam@86:   int i;
cannam@86:   int el=v->elements;
cannam@86:   float acc=0.f;
cannam@86:   for(i=0;i<el;i++){
cannam@86:     float val=(a[i]-b[i]);
cannam@86:     acc+=val*val;
cannam@86:   }
cannam@86:   return sqrt(acc);
cannam@86: }
cannam@86: 
cannam@86: float *_weight_null(vqgen *v,float *a){
cannam@86:   return a;
cannam@86: }
cannam@86: 
cannam@86: /* *must* be beefed up. */
cannam@86: void _vqgen_seed(vqgen *v){
cannam@86:   long i;
cannam@86:   for(i=0;i<v->entries;i++)
cannam@86:     memcpy(_now(v,i),_point(v,i),sizeof(float)*v->elements);
cannam@86:   v->seeded=1;
cannam@86: }
cannam@86: 
cannam@86: int directdsort(const void *a, const void *b){
cannam@86:   float av=*((float *)a);
cannam@86:   float bv=*((float *)b);
cannam@86:   return (av<bv)-(av>bv);
cannam@86: }
cannam@86: 
cannam@86: void vqgen_cellmetric(vqgen *v){
cannam@86:   int j,k;
cannam@86:   float min=-1.f,max=-1.f,mean=0.f,acc=0.f;
cannam@86:   long dup=0,unused=0;
cannam@86:  #ifdef NOISY
cannam@86:   int i;
cannam@86:    char buff[80];
cannam@86:    float spacings[v->entries];
cannam@86:    int count=0;
cannam@86:    FILE *cells;
cannam@86:    sprintf(buff,"cellspace%d.m",v->it);
cannam@86:    cells=fopen(buff,"w");
cannam@86: #endif
cannam@86: 
cannam@86:   /* minimum, maximum, cell spacing */
cannam@86:   for(j=0;j<v->entries;j++){
cannam@86:     float localmin=-1.;
cannam@86: 
cannam@86:     for(k=0;k<v->entries;k++){
cannam@86:       if(j!=k){
cannam@86:         float this=_dist(v,_now(v,j),_now(v,k));
cannam@86:         if(this>0){
cannam@86:           if(v->assigned[k] && (localmin==-1 || this<localmin))
cannam@86:             localmin=this;
cannam@86:         }else{        
cannam@86:           if(k<j){
cannam@86:             dup++;
cannam@86:             break;
cannam@86:           }
cannam@86:         }
cannam@86:       }
cannam@86:     }
cannam@86:     if(k<v->entries)continue;
cannam@86: 
cannam@86:     if(v->assigned[j]==0){
cannam@86:       unused++;
cannam@86:       continue;
cannam@86:     }
cannam@86:     
cannam@86:     localmin=v->max[j]+localmin/2; /* this gives us rough diameter */
cannam@86:     if(min==-1 || localmin<min)min=localmin;
cannam@86:     if(max==-1 || localmin>max)max=localmin;
cannam@86:     mean+=localmin;
cannam@86:     acc++;
cannam@86: #ifdef NOISY
cannam@86:     spacings[count++]=localmin;
cannam@86: #endif
cannam@86:   }
cannam@86: 
cannam@86:   fprintf(stderr,"cell diameter: %.03g::%.03g::%.03g (%ld unused/%ld dup)\n",
cannam@86:           min,mean/acc,max,unused,dup);
cannam@86: 
cannam@86: #ifdef NOISY
cannam@86:   qsort(spacings,count,sizeof(float),directdsort);
cannam@86:   for(i=0;i<count;i++)
cannam@86:     fprintf(cells,"%g\n",spacings[i]);
cannam@86:   fclose(cells);
cannam@86: #endif            
cannam@86: 
cannam@86: }
cannam@86: 
cannam@86: /* External calls *******************************************************/
cannam@86: 
cannam@86: /* We have two forms of quantization; in the first, each vector
cannam@86:    element in the codebook entry is orthogonal.  Residues would use this
cannam@86:    quantization for example.
cannam@86: 
cannam@86:    In the second, we have a sequence of monotonically increasing
cannam@86:    values that we wish to quantize as deltas (to save space).  We
cannam@86:    still need to quantize so that absolute values are accurate. For
cannam@86:    example, LSP quantizes all absolute values, but the book encodes
cannam@86:    distance between values because each successive value is larger
cannam@86:    than the preceeding value.  Thus the desired quantibits apply to
cannam@86:    the encoded (delta) values, not abs positions. This requires minor
cannam@86:    additional encode-side trickery. */
cannam@86: 
cannam@86: void vqgen_quantize(vqgen *v,quant_meta *q){
cannam@86: 
cannam@86:   float maxdel;
cannam@86:   float mindel;
cannam@86: 
cannam@86:   float delta;
cannam@86:   float maxquant=((1<<q->quant)-1);
cannam@86: 
cannam@86:   int j,k;
cannam@86: 
cannam@86:   mindel=maxdel=_now(v,0)[0];
cannam@86:   
cannam@86:   for(j=0;j<v->entries;j++){
cannam@86:     float last=0.f;
cannam@86:     for(k=0;k<v->elements;k++){
cannam@86:       if(mindel>_now(v,j)[k]-last)mindel=_now(v,j)[k]-last;
cannam@86:       if(maxdel<_now(v,j)[k]-last)maxdel=_now(v,j)[k]-last;
cannam@86:       if(q->sequencep)last=_now(v,j)[k];
cannam@86:     }
cannam@86:   }
cannam@86: 
cannam@86: 
cannam@86:   /* first find the basic delta amount from the maximum span to be
cannam@86:      encoded.  Loosen the delta slightly to allow for additional error
cannam@86:      during sequence quantization */
cannam@86: 
cannam@86:   delta=(maxdel-mindel)/((1<<q->quant)-1.5f);
cannam@86: 
cannam@86:   q->min=_float32_pack(mindel);
cannam@86:   q->delta=_float32_pack(delta);
cannam@86: 
cannam@86:   mindel=_float32_unpack(q->min);
cannam@86:   delta=_float32_unpack(q->delta);
cannam@86: 
cannam@86:   for(j=0;j<v->entries;j++){
cannam@86:     float last=0;
cannam@86:     for(k=0;k<v->elements;k++){
cannam@86:       float val=_now(v,j)[k];
cannam@86:       float now=rint((val-last-mindel)/delta);
cannam@86:       
cannam@86:       _now(v,j)[k]=now;
cannam@86:       if(now<0){
cannam@86:         /* be paranoid; this should be impossible */
cannam@86:         fprintf(stderr,"fault; quantized value<0\n");
cannam@86:         exit(1);
cannam@86:       }
cannam@86: 
cannam@86:       if(now>maxquant){
cannam@86:         /* be paranoid; this should be impossible */
cannam@86:         fprintf(stderr,"fault; quantized value>max\n");
cannam@86:         exit(1);
cannam@86:       }
cannam@86:       if(q->sequencep)last=(now*delta)+mindel+last;
cannam@86:     }
cannam@86:   }
cannam@86: }
cannam@86: 
cannam@86: /* much easier :-).  Unlike in the codebook, we don't un-log log
cannam@86:    scales; we just make sure they're properly offset. */
cannam@86: void vqgen_unquantize(vqgen *v,quant_meta *q){
cannam@86:   long j,k;
cannam@86:   float mindel=_float32_unpack(q->min);
cannam@86:   float delta=_float32_unpack(q->delta);
cannam@86: 
cannam@86:   for(j=0;j<v->entries;j++){
cannam@86:     float last=0.f;
cannam@86:     for(k=0;k<v->elements;k++){
cannam@86:       float now=_now(v,j)[k];
cannam@86:       now=fabs(now)*delta+last+mindel;
cannam@86:       if(q->sequencep)last=now;
cannam@86:       _now(v,j)[k]=now;
cannam@86:     }
cannam@86:   }
cannam@86: }
cannam@86: 
cannam@86: void vqgen_init(vqgen *v,int elements,int aux,int entries,float mindist,
cannam@86:                 float  (*metric)(vqgen *,float *, float *),
cannam@86:                 float *(*weight)(vqgen *,float *),int centroid){
cannam@86:   memset(v,0,sizeof(vqgen));
cannam@86: 
cannam@86:   v->centroid=centroid;
cannam@86:   v->elements=elements;
cannam@86:   v->aux=aux;
cannam@86:   v->mindist=mindist;
cannam@86:   v->allocated=32768;
cannam@86:   v->pointlist=_ogg_malloc(v->allocated*(v->elements+v->aux)*sizeof(float));
cannam@86: 
cannam@86:   v->entries=entries;
cannam@86:   v->entrylist=_ogg_malloc(v->entries*v->elements*sizeof(float));
cannam@86:   v->assigned=_ogg_malloc(v->entries*sizeof(long));
cannam@86:   v->bias=_ogg_calloc(v->entries,sizeof(float));
cannam@86:   v->max=_ogg_calloc(v->entries,sizeof(float));
cannam@86:   if(metric)
cannam@86:     v->metric_func=metric;
cannam@86:   else
cannam@86:     v->metric_func=_dist;
cannam@86:   if(weight)
cannam@86:     v->weight_func=weight;
cannam@86:   else
cannam@86:     v->weight_func=_weight_null;
cannam@86: 
cannam@86:   v->asciipoints=tmpfile();
cannam@86: 
cannam@86: }
cannam@86: 
cannam@86: void vqgen_addpoint(vqgen *v, float *p,float *a){
cannam@86:   int k;
cannam@86:   for(k=0;k<v->elements;k++)
cannam@86:     fprintf(v->asciipoints,"%.12g\n",p[k]);
cannam@86:   for(k=0;k<v->aux;k++)
cannam@86:     fprintf(v->asciipoints,"%.12g\n",a[k]);
cannam@86: 
cannam@86:   if(v->points>=v->allocated){
cannam@86:     v->allocated*=2;
cannam@86:     v->pointlist=_ogg_realloc(v->pointlist,v->allocated*(v->elements+v->aux)*
cannam@86:                          sizeof(float));
cannam@86:   }
cannam@86: 
cannam@86:   memcpy(_point(v,v->points),p,sizeof(float)*v->elements);
cannam@86:   if(v->aux)memcpy(_point(v,v->points)+v->elements,a,sizeof(float)*v->aux);
cannam@86:  
cannam@86:   /* quantize to the density mesh if it's selected */
cannam@86:   if(v->mindist>0.f){
cannam@86:     /* quantize to the mesh */
cannam@86:     for(k=0;k<v->elements+v->aux;k++)
cannam@86:       _point(v,v->points)[k]=
cannam@86:         rint(_point(v,v->points)[k]/v->mindist)*v->mindist;
cannam@86:   }
cannam@86:   v->points++;
cannam@86:   if(!(v->points&0xff))spinnit("loading... ",v->points);
cannam@86: }
cannam@86: 
cannam@86: /* yes, not threadsafe.  These utils aren't */
cannam@86: static int sortit=0;
cannam@86: static int sortsize=0;
cannam@86: static int meshcomp(const void *a,const void *b){
cannam@86:   if(((sortit++)&0xfff)==0)spinnit("sorting mesh...",sortit);
cannam@86:   return(memcmp(a,b,sortsize));
cannam@86: }
cannam@86: 
cannam@86: void vqgen_sortmesh(vqgen *v){
cannam@86:   sortit=0;
cannam@86:   if(v->mindist>0.f){
cannam@86:     long i,march=1;
cannam@86: 
cannam@86:     /* sort to make uniqueness detection trivial */
cannam@86:     sortsize=(v->elements+v->aux)*sizeof(float);
cannam@86:     qsort(v->pointlist,v->points,sortsize,meshcomp);
cannam@86: 
cannam@86:     /* now march through and eliminate dupes */
cannam@86:     for(i=1;i<v->points;i++){
cannam@86:       if(memcmp(_point(v,i),_point(v,i-1),sortsize)){
cannam@86:         /* a new, unique entry.  march it down */
cannam@86:         if(i>march)memcpy(_point(v,march),_point(v,i),sortsize);
cannam@86:         march++;
cannam@86:       }
cannam@86:       spinnit("eliminating density... ",v->points-i);
cannam@86:     }
cannam@86: 
cannam@86:     /* we're done */
cannam@86:     fprintf(stderr,"\r%ld training points remining out of %ld"
cannam@86:             " after density mesh (%ld%%)\n",march,v->points,march*100/v->points);
cannam@86:     v->points=march;
cannam@86: 
cannam@86:   }
cannam@86:   v->sorted=1;
cannam@86: }
cannam@86: 
cannam@86: float vqgen_iterate(vqgen *v,int biasp){
cannam@86:   long   i,j,k;
cannam@86: 
cannam@86:   float fdesired;
cannam@86:   long  desired;
cannam@86:   long  desired2;
cannam@86: 
cannam@86:   float asserror=0.f;
cannam@86:   float meterror=0.f;
cannam@86:   float *new;
cannam@86:   float *new2;
cannam@86:   long   *nearcount;
cannam@86:   float *nearbias;
cannam@86:  #ifdef NOISY
cannam@86:    char buff[80];
cannam@86:    FILE *assig;
cannam@86:    FILE *bias;
cannam@86:    FILE *cells;
cannam@86:    sprintf(buff,"cells%d.m",v->it);
cannam@86:    cells=fopen(buff,"w");
cannam@86:    sprintf(buff,"assig%d.m",v->it);
cannam@86:    assig=fopen(buff,"w");
cannam@86:    sprintf(buff,"bias%d.m",v->it);
cannam@86:    bias=fopen(buff,"w");
cannam@86:  #endif
cannam@86:  
cannam@86: 
cannam@86:   if(v->entries<2){
cannam@86:     fprintf(stderr,"generation requires at least two entries\n");
cannam@86:     exit(1);
cannam@86:   }
cannam@86: 
cannam@86:   if(!v->sorted)vqgen_sortmesh(v);
cannam@86:   if(!v->seeded)_vqgen_seed(v);
cannam@86: 
cannam@86:   fdesired=(float)v->points/v->entries;
cannam@86:   desired=fdesired;
cannam@86:   desired2=desired*2;
cannam@86:   new=_ogg_malloc(sizeof(float)*v->entries*v->elements);
cannam@86:   new2=_ogg_malloc(sizeof(float)*v->entries*v->elements);
cannam@86:   nearcount=_ogg_malloc(v->entries*sizeof(long));
cannam@86:   nearbias=_ogg_malloc(v->entries*desired2*sizeof(float));
cannam@86: 
cannam@86:   /* fill in nearest points for entry biasing */
cannam@86:   /*memset(v->bias,0,sizeof(float)*v->entries);*/
cannam@86:   memset(nearcount,0,sizeof(long)*v->entries);
cannam@86:   memset(v->assigned,0,sizeof(long)*v->entries);
cannam@86:   if(biasp){
cannam@86:     for(i=0;i<v->points;i++){
cannam@86:       float *ppt=v->weight_func(v,_point(v,i));
cannam@86:       float firstmetric=v->metric_func(v,_now(v,0),ppt)+v->bias[0];
cannam@86:       float secondmetric=v->metric_func(v,_now(v,1),ppt)+v->bias[1];
cannam@86:       long   firstentry=0;
cannam@86:       long   secondentry=1;
cannam@86:       
cannam@86:       if(!(i&0xff))spinnit("biasing... ",v->points+v->points+v->entries-i);
cannam@86:       
cannam@86:       if(firstmetric>secondmetric){
cannam@86:         float temp=firstmetric;
cannam@86:         firstmetric=secondmetric;
cannam@86:         secondmetric=temp;
cannam@86:         firstentry=1;
cannam@86:         secondentry=0;
cannam@86:       }
cannam@86:       
cannam@86:       for(j=2;j<v->entries;j++){
cannam@86:         float thismetric=v->metric_func(v,_now(v,j),ppt)+v->bias[j];
cannam@86:         if(thismetric<secondmetric){
cannam@86:           if(thismetric<firstmetric){
cannam@86:             secondmetric=firstmetric;
cannam@86:             secondentry=firstentry;
cannam@86:             firstmetric=thismetric;
cannam@86:             firstentry=j;
cannam@86:           }else{
cannam@86:             secondmetric=thismetric;
cannam@86:             secondentry=j;
cannam@86:           }
cannam@86:         }
cannam@86:       }
cannam@86:       
cannam@86:       j=firstentry;
cannam@86:       for(j=0;j<v->entries;j++){
cannam@86:         
cannam@86:         float thismetric,localmetric;
cannam@86:         float *nearbiasptr=nearbias+desired2*j;
cannam@86:         long k=nearcount[j];
cannam@86:         
cannam@86:         localmetric=v->metric_func(v,_now(v,j),ppt);
cannam@86:         /* 'thismetric' is to be the bias value necessary in the current
cannam@86:            arrangement for entry j to capture point i */
cannam@86:         if(firstentry==j){
cannam@86:           /* use the secondary entry as the threshhold */
cannam@86:           thismetric=secondmetric-localmetric;
cannam@86:         }else{
cannam@86:           /* use the primary entry as the threshhold */
cannam@86:           thismetric=firstmetric-localmetric;
cannam@86:         }
cannam@86:         
cannam@86:         /* support the idea of 'minimum distance'... if we want the
cannam@86:            cells in a codebook to be roughly some minimum size (as with
cannam@86:            the low resolution residue books) */
cannam@86:         
cannam@86:         /* a cute two-stage delayed sorting hack */
cannam@86:         if(k<desired){
cannam@86:           nearbiasptr[k]=thismetric;
cannam@86:           k++;
cannam@86:           if(k==desired){
cannam@86:             spinnit("biasing... ",v->points+v->points+v->entries-i);
cannam@86:             qsort(nearbiasptr,desired,sizeof(float),directdsort);
cannam@86:           }
cannam@86:           
cannam@86:         }else if(thismetric>nearbiasptr[desired-1]){
cannam@86:           nearbiasptr[k]=thismetric;
cannam@86:           k++;
cannam@86:           if(k==desired2){
cannam@86:             spinnit("biasing... ",v->points+v->points+v->entries-i);
cannam@86:             qsort(nearbiasptr,desired2,sizeof(float),directdsort);
cannam@86:             k=desired;
cannam@86:           }
cannam@86:         }
cannam@86:         nearcount[j]=k;
cannam@86:       }
cannam@86:     }
cannam@86:     
cannam@86:     /* inflate/deflate */
cannam@86:     
cannam@86:     for(i=0;i<v->entries;i++){
cannam@86:       float *nearbiasptr=nearbias+desired2*i;
cannam@86:       
cannam@86:       spinnit("biasing... ",v->points+v->entries-i);
cannam@86:       
cannam@86:       /* due to the delayed sorting, we likely need to finish it off....*/
cannam@86:       if(nearcount[i]>desired)
cannam@86:         qsort(nearbiasptr,nearcount[i],sizeof(float),directdsort);
cannam@86: 
cannam@86:       v->bias[i]=nearbiasptr[desired-1];
cannam@86: 
cannam@86:     }
cannam@86:   }else{ 
cannam@86:     memset(v->bias,0,v->entries*sizeof(float));
cannam@86:   }
cannam@86: 
cannam@86:   /* Now assign with new bias and find new midpoints */
cannam@86:   for(i=0;i<v->points;i++){
cannam@86:     float *ppt=v->weight_func(v,_point(v,i));
cannam@86:     float firstmetric=v->metric_func(v,_now(v,0),ppt)+v->bias[0];
cannam@86:     long   firstentry=0;
cannam@86: 
cannam@86:     if(!(i&0xff))spinnit("centering... ",v->points-i);
cannam@86: 
cannam@86:     for(j=0;j<v->entries;j++){
cannam@86:       float thismetric=v->metric_func(v,_now(v,j),ppt)+v->bias[j];
cannam@86:       if(thismetric<firstmetric){
cannam@86:         firstmetric=thismetric;
cannam@86:         firstentry=j;
cannam@86:       }
cannam@86:     }
cannam@86: 
cannam@86:     j=firstentry;
cannam@86:       
cannam@86: #ifdef NOISY
cannam@86:     fprintf(cells,"%g %g\n%g %g\n\n",
cannam@86:           _now(v,j)[0],_now(v,j)[1],
cannam@86:           ppt[0],ppt[1]);
cannam@86: #endif
cannam@86: 
cannam@86:     firstmetric-=v->bias[j];
cannam@86:     meterror+=firstmetric;
cannam@86: 
cannam@86:     if(v->centroid==0){
cannam@86:       /* set up midpoints for next iter */
cannam@86:       if(v->assigned[j]++){
cannam@86:         for(k=0;k<v->elements;k++)
cannam@86:           vN(new,j)[k]+=ppt[k];
cannam@86:         if(firstmetric>v->max[j])v->max[j]=firstmetric;
cannam@86:       }else{
cannam@86:         for(k=0;k<v->elements;k++)
cannam@86:           vN(new,j)[k]=ppt[k];
cannam@86:         v->max[j]=firstmetric;
cannam@86:       }
cannam@86:     }else{
cannam@86:       /* centroid */
cannam@86:       if(v->assigned[j]++){
cannam@86:         for(k=0;k<v->elements;k++){
cannam@86:           if(vN(new,j)[k]>ppt[k])vN(new,j)[k]=ppt[k];
cannam@86:           if(vN(new2,j)[k]<ppt[k])vN(new2,j)[k]=ppt[k];
cannam@86:         }
cannam@86:         if(firstmetric>v->max[firstentry])v->max[j]=firstmetric;
cannam@86:       }else{
cannam@86:         for(k=0;k<v->elements;k++){
cannam@86:           vN(new,j)[k]=ppt[k];
cannam@86:           vN(new2,j)[k]=ppt[k];
cannam@86:         }
cannam@86:         v->max[firstentry]=firstmetric;
cannam@86:       }
cannam@86:     }
cannam@86:   }
cannam@86: 
cannam@86:   /* assign midpoints */
cannam@86: 
cannam@86:   for(j=0;j<v->entries;j++){
cannam@86: #ifdef NOISY
cannam@86:     fprintf(assig,"%ld\n",v->assigned[j]);
cannam@86:     fprintf(bias,"%g\n",v->bias[j]);
cannam@86: #endif
cannam@86:     asserror+=fabs(v->assigned[j]-fdesired);
cannam@86:     if(v->assigned[j]){
cannam@86:       if(v->centroid==0){
cannam@86:         for(k=0;k<v->elements;k++)
cannam@86:           _now(v,j)[k]=vN(new,j)[k]/v->assigned[j];
cannam@86:       }else{
cannam@86:         for(k=0;k<v->elements;k++)
cannam@86:           _now(v,j)[k]=(vN(new,j)[k]+vN(new2,j)[k])/2.f;
cannam@86:       }
cannam@86:     }
cannam@86:   }
cannam@86: 
cannam@86:   asserror/=(v->entries*fdesired);
cannam@86: 
cannam@86:   fprintf(stderr,"Pass #%d... ",v->it);
cannam@86:   fprintf(stderr,": dist %g(%g) metric error=%g \n",
cannam@86:           asserror,fdesired,meterror/v->points);
cannam@86:   v->it++;
cannam@86:   
cannam@86:   free(new);
cannam@86:   free(nearcount);
cannam@86:   free(nearbias);
cannam@86: #ifdef NOISY
cannam@86:   fclose(assig);
cannam@86:   fclose(bias);
cannam@86:   fclose(cells);
cannam@86: #endif
cannam@86:   return(asserror);
cannam@86: }
cannam@86: