18e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels/******************************************************************** 28e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels * * 38e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE. * 48e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS * 58e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE * 68e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING. * 78e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels * * 88e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2001 * 98e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels * by the Xiph.Org Foundation http://www.xiph.org/ * 108e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels * * 118e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels ******************************************************************** 128e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 138e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels function: train a VQ codebook 148e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels last mod: $Id: vqgen.c 16037 2009-05-26 21:10:58Z xiphmont $ 158e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 168e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels ********************************************************************/ 178e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 188e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels/* This code is *not* part of libvorbis. It is used to generate 198e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels trained codebooks offline and then spit the results into a 208e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels pregenerated codebook that is compiled into libvorbis. It is an 218e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels expensive (but good) algorithm. Run it on big iron. */ 228e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 238e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels/* There are so many optimizations to explore in *both* stages that 248e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels considering the undertaking is almost withering. For now, we brute 258e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels force it all */ 268e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 278e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels#include <stdlib.h> 288e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels#include <stdio.h> 298e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels#include <math.h> 308e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels#include <string.h> 318e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 328e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels#include "vqgen.h" 338e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels#include "bookutil.h" 348e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 358e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels/* Codebook generation happens in two steps: 368e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 378e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 1) Train the codebook with data collected from the encoder: We use 388e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels one of a few error metrics (which represent the distance between a 398e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels given data point and a candidate point in the training set) to 408e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels divide the training set up into cells representing roughly equal 418e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels probability of occurring. 428e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 438e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 2) Generate the codebook and auxiliary data from the trained data set 448e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels*/ 458e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 468e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels/* Codebook training **************************************************** 478e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels * 488e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels * The basic idea here is that a VQ codebook is like an m-dimensional 498e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels * foam with n bubbles. The bubbles compete for space/volume and are 508e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels * 'pressurized' [biased] according to some metric. The basic alg 518e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels * iterates through allowing the bubbles to compete for space until 528e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels * they converge (if the damping is dome properly) on a steady-state 538e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels * solution. Individual input points, collected from libvorbis, are 548e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels * used to train the algorithm monte-carlo style. */ 558e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 568e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels/* internal helpers *****************************************************/ 578e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels#define vN(data,i) (data+v->elements*i) 588e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 598e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels/* default metric; squared 'distance' from desired value. */ 608e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckelsfloat _dist(vqgen *v,float *a, float *b){ 618e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels int i; 628e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels int el=v->elements; 638e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels float acc=0.f; 648e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels for(i=0;i<el;i++){ 658e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels float val=(a[i]-b[i]); 668e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels acc+=val*val; 678e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels } 688e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels return sqrt(acc); 698e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels} 708e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 718e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckelsfloat *_weight_null(vqgen *v,float *a){ 728e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels return a; 738e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels} 748e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 758e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels/* *must* be beefed up. */ 768e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckelsvoid _vqgen_seed(vqgen *v){ 778e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels long i; 788e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels for(i=0;i<v->entries;i++) 798e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels memcpy(_now(v,i),_point(v,i),sizeof(float)*v->elements); 808e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels v->seeded=1; 818e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels} 828e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 838e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckelsint directdsort(const void *a, const void *b){ 848e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels float av=*((float *)a); 858e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels float bv=*((float *)b); 868e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels return (av<bv)-(av>bv); 878e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels} 888e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 898e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckelsvoid vqgen_cellmetric(vqgen *v){ 908e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels int j,k; 918e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels float min=-1.f,max=-1.f,mean=0.f,acc=0.f; 928e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels long dup=0,unused=0; 938e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels #ifdef NOISY 948e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels int i; 958e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels char buff[80]; 968e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels float spacings[v->entries]; 978e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels int count=0; 988e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels FILE *cells; 998e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels sprintf(buff,"cellspace%d.m",v->it); 1008e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels cells=fopen(buff,"w"); 1018e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels#endif 1028e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 1038e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels /* minimum, maximum, cell spacing */ 1048e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels for(j=0;j<v->entries;j++){ 1058e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels float localmin=-1.; 1068e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 1078e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels for(k=0;k<v->entries;k++){ 1088e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels if(j!=k){ 1098e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels float this=_dist(v,_now(v,j),_now(v,k)); 1108e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels if(this>0){ 1118e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels if(v->assigned[k] && (localmin==-1 || this<localmin)) 1128e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels localmin=this; 1138e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels }else{ 1148e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels if(k<j){ 1158e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels dup++; 1168e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels break; 1178e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels } 1188e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels } 1198e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels } 1208e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels } 1218e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels if(k<v->entries)continue; 1228e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 1238e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels if(v->assigned[j]==0){ 1248e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels unused++; 1258e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels continue; 1268e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels } 1278e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 1288e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels localmin=v->max[j]+localmin/2; /* this gives us rough diameter */ 1298e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels if(min==-1 || localmin<min)min=localmin; 1308e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels if(max==-1 || localmin>max)max=localmin; 1318e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels mean+=localmin; 1328e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels acc++; 1338e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels#ifdef NOISY 1348e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels spacings[count++]=localmin; 1358e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels#endif 1368e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels } 1378e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 1388e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels fprintf(stderr,"cell diameter: %.03g::%.03g::%.03g (%ld unused/%ld dup)\n", 1398e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels min,mean/acc,max,unused,dup); 1408e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 1418e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels#ifdef NOISY 1428e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels qsort(spacings,count,sizeof(float),directdsort); 1438e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels for(i=0;i<count;i++) 1448e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels fprintf(cells,"%g\n",spacings[i]); 1458e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels fclose(cells); 1468e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels#endif 1478e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 1488e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels} 1498e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 1508e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels/* External calls *******************************************************/ 1518e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 1528e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels/* We have two forms of quantization; in the first, each vector 1538e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels element in the codebook entry is orthogonal. Residues would use this 1548e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels quantization for example. 1558e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 1568e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels In the second, we have a sequence of monotonically increasing 1578e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels values that we wish to quantize as deltas (to save space). We 1588e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels still need to quantize so that absolute values are accurate. For 1598e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels example, LSP quantizes all absolute values, but the book encodes 1608e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels distance between values because each successive value is larger 1618e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels than the preceeding value. Thus the desired quantibits apply to 1628e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels the encoded (delta) values, not abs positions. This requires minor 1638e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels additional encode-side trickery. */ 1648e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 1658e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckelsvoid vqgen_quantize(vqgen *v,quant_meta *q){ 1668e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 1678e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels float maxdel; 1688e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels float mindel; 1698e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 1708e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels float delta; 1718e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels float maxquant=((1<<q->quant)-1); 1728e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 1738e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels int j,k; 1748e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 1758e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels mindel=maxdel=_now(v,0)[0]; 1768e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 1778e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels for(j=0;j<v->entries;j++){ 1788e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels float last=0.f; 1798e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels for(k=0;k<v->elements;k++){ 1808e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels if(mindel>_now(v,j)[k]-last)mindel=_now(v,j)[k]-last; 1818e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels if(maxdel<_now(v,j)[k]-last)maxdel=_now(v,j)[k]-last; 1828e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels if(q->sequencep)last=_now(v,j)[k]; 1838e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels } 1848e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels } 1858e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 1868e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 1878e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels /* first find the basic delta amount from the maximum span to be 1888e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels encoded. Loosen the delta slightly to allow for additional error 1898e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels during sequence quantization */ 1908e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 1918e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels delta=(maxdel-mindel)/((1<<q->quant)-1.5f); 1928e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 1938e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels q->min=_float32_pack(mindel); 1948e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels q->delta=_float32_pack(delta); 1958e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 1968e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels mindel=_float32_unpack(q->min); 1978e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels delta=_float32_unpack(q->delta); 1988e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 1998e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels for(j=0;j<v->entries;j++){ 2008e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels float last=0; 2018e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels for(k=0;k<v->elements;k++){ 2028e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels float val=_now(v,j)[k]; 2038e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels float now=rint((val-last-mindel)/delta); 2048e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 2058e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels _now(v,j)[k]=now; 2068e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels if(now<0){ 2078e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels /* be paranoid; this should be impossible */ 2088e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels fprintf(stderr,"fault; quantized value<0\n"); 2098e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels exit(1); 2108e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels } 2118e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 2128e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels if(now>maxquant){ 2138e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels /* be paranoid; this should be impossible */ 2148e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels fprintf(stderr,"fault; quantized value>max\n"); 2158e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels exit(1); 2168e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels } 2178e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels if(q->sequencep)last=(now*delta)+mindel+last; 2188e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels } 2198e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels } 2208e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels} 2218e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 2228e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels/* much easier :-). Unlike in the codebook, we don't un-log log 2238e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels scales; we just make sure they're properly offset. */ 2248e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckelsvoid vqgen_unquantize(vqgen *v,quant_meta *q){ 2258e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels long j,k; 2268e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels float mindel=_float32_unpack(q->min); 2278e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels float delta=_float32_unpack(q->delta); 2288e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 2298e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels for(j=0;j<v->entries;j++){ 2308e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels float last=0.f; 2318e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels for(k=0;k<v->elements;k++){ 2328e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels float now=_now(v,j)[k]; 2338e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels now=fabs(now)*delta+last+mindel; 2348e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels if(q->sequencep)last=now; 2358e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels _now(v,j)[k]=now; 2368e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels } 2378e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels } 2388e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels} 2398e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 2408e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckelsvoid vqgen_init(vqgen *v,int elements,int aux,int entries,float mindist, 2418e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels float (*metric)(vqgen *,float *, float *), 2428e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels float *(*weight)(vqgen *,float *),int centroid){ 2438e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels memset(v,0,sizeof(vqgen)); 2448e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 2458e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels v->centroid=centroid; 2468e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels v->elements=elements; 2478e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels v->aux=aux; 2488e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels v->mindist=mindist; 2498e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels v->allocated=32768; 2508e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels v->pointlist=_ogg_malloc(v->allocated*(v->elements+v->aux)*sizeof(float)); 2518e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 2528e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels v->entries=entries; 2538e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels v->entrylist=_ogg_malloc(v->entries*v->elements*sizeof(float)); 2548e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels v->assigned=_ogg_malloc(v->entries*sizeof(long)); 2558e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels v->bias=_ogg_calloc(v->entries,sizeof(float)); 2568e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels v->max=_ogg_calloc(v->entries,sizeof(float)); 2578e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels if(metric) 2588e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels v->metric_func=metric; 2598e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels else 2608e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels v->metric_func=_dist; 2618e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels if(weight) 2628e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels v->weight_func=weight; 2638e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels else 2648e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels v->weight_func=_weight_null; 2658e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 2668e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels v->asciipoints=tmpfile(); 2678e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 2688e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels} 2698e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 2708e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckelsvoid vqgen_addpoint(vqgen *v, float *p,float *a){ 2718e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels int k; 2728e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels for(k=0;k<v->elements;k++) 2738e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels fprintf(v->asciipoints,"%.12g\n",p[k]); 2748e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels for(k=0;k<v->aux;k++) 2758e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels fprintf(v->asciipoints,"%.12g\n",a[k]); 2768e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 2778e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels if(v->points>=v->allocated){ 2788e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels v->allocated*=2; 2798e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels v->pointlist=_ogg_realloc(v->pointlist,v->allocated*(v->elements+v->aux)* 2808e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels sizeof(float)); 2818e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels } 2828e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 2838e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels memcpy(_point(v,v->points),p,sizeof(float)*v->elements); 2848e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels if(v->aux)memcpy(_point(v,v->points)+v->elements,a,sizeof(float)*v->aux); 2858e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 2868e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels /* quantize to the density mesh if it's selected */ 2878e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels if(v->mindist>0.f){ 2888e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels /* quantize to the mesh */ 2898e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels for(k=0;k<v->elements+v->aux;k++) 2908e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels _point(v,v->points)[k]= 2918e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels rint(_point(v,v->points)[k]/v->mindist)*v->mindist; 2928e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels } 2938e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels v->points++; 2948e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels if(!(v->points&0xff))spinnit("loading... ",v->points); 2958e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels} 2968e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 2978e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels/* yes, not threadsafe. These utils aren't */ 2988e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckelsstatic int sortit=0; 2998e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckelsstatic int sortsize=0; 3008e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckelsstatic int meshcomp(const void *a,const void *b){ 3018e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels if(((sortit++)&0xfff)==0)spinnit("sorting mesh...",sortit); 3028e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels return(memcmp(a,b,sortsize)); 3038e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels} 3048e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 3058e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckelsvoid vqgen_sortmesh(vqgen *v){ 3068e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels sortit=0; 3078e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels if(v->mindist>0.f){ 3088e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels long i,march=1; 3098e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 3108e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels /* sort to make uniqueness detection trivial */ 3118e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels sortsize=(v->elements+v->aux)*sizeof(float); 3128e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels qsort(v->pointlist,v->points,sortsize,meshcomp); 3138e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 3148e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels /* now march through and eliminate dupes */ 3158e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels for(i=1;i<v->points;i++){ 3168e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels if(memcmp(_point(v,i),_point(v,i-1),sortsize)){ 3178e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels /* a new, unique entry. march it down */ 3188e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels if(i>march)memcpy(_point(v,march),_point(v,i),sortsize); 3198e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels march++; 3208e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels } 3218e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels spinnit("eliminating density... ",v->points-i); 3228e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels } 3238e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 3248e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels /* we're done */ 3258e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels fprintf(stderr,"\r%ld training points remining out of %ld" 3268e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels " after density mesh (%ld%%)\n",march,v->points,march*100/v->points); 3278e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels v->points=march; 3288e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 3298e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels } 3308e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels v->sorted=1; 3318e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels} 3328e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 3338e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckelsfloat vqgen_iterate(vqgen *v,int biasp){ 3348e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels long i,j,k; 3358e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 3368e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels float fdesired; 3378e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels long desired; 3388e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels long desired2; 3398e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 3408e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels float asserror=0.f; 3418e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels float meterror=0.f; 3428e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels float *new; 3438e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels float *new2; 3448e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels long *nearcount; 3458e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels float *nearbias; 3468e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels #ifdef NOISY 3478e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels char buff[80]; 3488e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels FILE *assig; 3498e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels FILE *bias; 3508e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels FILE *cells; 3518e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels sprintf(buff,"cells%d.m",v->it); 3528e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels cells=fopen(buff,"w"); 3538e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels sprintf(buff,"assig%d.m",v->it); 3548e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels assig=fopen(buff,"w"); 3558e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels sprintf(buff,"bias%d.m",v->it); 3568e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels bias=fopen(buff,"w"); 3578e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels #endif 3588e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 3598e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 3608e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels if(v->entries<2){ 3618e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels fprintf(stderr,"generation requires at least two entries\n"); 3628e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels exit(1); 3638e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels } 3648e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 3658e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels if(!v->sorted)vqgen_sortmesh(v); 3668e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels if(!v->seeded)_vqgen_seed(v); 3678e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 3688e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels fdesired=(float)v->points/v->entries; 3698e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels desired=fdesired; 3708e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels desired2=desired*2; 3718e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels new=_ogg_malloc(sizeof(float)*v->entries*v->elements); 3728e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels new2=_ogg_malloc(sizeof(float)*v->entries*v->elements); 3738e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels nearcount=_ogg_malloc(v->entries*sizeof(long)); 3748e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels nearbias=_ogg_malloc(v->entries*desired2*sizeof(float)); 3758e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 3768e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels /* fill in nearest points for entry biasing */ 3778e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels /*memset(v->bias,0,sizeof(float)*v->entries);*/ 3788e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels memset(nearcount,0,sizeof(long)*v->entries); 3798e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels memset(v->assigned,0,sizeof(long)*v->entries); 3808e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels if(biasp){ 3818e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels for(i=0;i<v->points;i++){ 3828e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels float *ppt=v->weight_func(v,_point(v,i)); 3838e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels float firstmetric=v->metric_func(v,_now(v,0),ppt)+v->bias[0]; 3848e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels float secondmetric=v->metric_func(v,_now(v,1),ppt)+v->bias[1]; 3858e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels long firstentry=0; 3868e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels long secondentry=1; 3878e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 3888e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels if(!(i&0xff))spinnit("biasing... ",v->points+v->points+v->entries-i); 3898e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 3908e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels if(firstmetric>secondmetric){ 3918e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels float temp=firstmetric; 3928e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels firstmetric=secondmetric; 3938e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels secondmetric=temp; 3948e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels firstentry=1; 3958e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels secondentry=0; 3968e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels } 3978e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 3988e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels for(j=2;j<v->entries;j++){ 3998e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels float thismetric=v->metric_func(v,_now(v,j),ppt)+v->bias[j]; 4008e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels if(thismetric<secondmetric){ 4018e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels if(thismetric<firstmetric){ 4028e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels secondmetric=firstmetric; 4038e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels secondentry=firstentry; 4048e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels firstmetric=thismetric; 4058e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels firstentry=j; 4068e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels }else{ 4078e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels secondmetric=thismetric; 4088e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels secondentry=j; 4098e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels } 4108e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels } 4118e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels } 4128e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 4138e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels j=firstentry; 4148e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels for(j=0;j<v->entries;j++){ 4158e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 4168e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels float thismetric,localmetric; 4178e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels float *nearbiasptr=nearbias+desired2*j; 4188e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels long k=nearcount[j]; 4198e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 4208e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels localmetric=v->metric_func(v,_now(v,j),ppt); 4218e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels /* 'thismetric' is to be the bias value necessary in the current 4228e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels arrangement for entry j to capture point i */ 4238e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels if(firstentry==j){ 4248e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels /* use the secondary entry as the threshhold */ 4258e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels thismetric=secondmetric-localmetric; 4268e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels }else{ 4278e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels /* use the primary entry as the threshhold */ 4288e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels thismetric=firstmetric-localmetric; 4298e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels } 4308e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 4318e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels /* support the idea of 'minimum distance'... if we want the 4328e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels cells in a codebook to be roughly some minimum size (as with 4338e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels the low resolution residue books) */ 4348e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 4358e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels /* a cute two-stage delayed sorting hack */ 4368e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels if(k<desired){ 4378e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels nearbiasptr[k]=thismetric; 4388e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels k++; 4398e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels if(k==desired){ 4408e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels spinnit("biasing... ",v->points+v->points+v->entries-i); 4418e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels qsort(nearbiasptr,desired,sizeof(float),directdsort); 4428e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels } 4438e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 4448e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels }else if(thismetric>nearbiasptr[desired-1]){ 4458e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels nearbiasptr[k]=thismetric; 4468e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels k++; 4478e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels if(k==desired2){ 4488e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels spinnit("biasing... ",v->points+v->points+v->entries-i); 4498e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels qsort(nearbiasptr,desired2,sizeof(float),directdsort); 4508e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels k=desired; 4518e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels } 4528e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels } 4538e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels nearcount[j]=k; 4548e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels } 4558e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels } 4568e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 4578e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels /* inflate/deflate */ 4588e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 4598e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels for(i=0;i<v->entries;i++){ 4608e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels float *nearbiasptr=nearbias+desired2*i; 4618e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 4628e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels spinnit("biasing... ",v->points+v->entries-i); 4638e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 4648e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels /* due to the delayed sorting, we likely need to finish it off....*/ 4658e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels if(nearcount[i]>desired) 4668e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels qsort(nearbiasptr,nearcount[i],sizeof(float),directdsort); 4678e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 4688e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels v->bias[i]=nearbiasptr[desired-1]; 4698e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 4708e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels } 4718e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels }else{ 4728e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels memset(v->bias,0,v->entries*sizeof(float)); 4738e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels } 4748e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 4758e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels /* Now assign with new bias and find new midpoints */ 4768e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels for(i=0;i<v->points;i++){ 4778e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels float *ppt=v->weight_func(v,_point(v,i)); 4788e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels float firstmetric=v->metric_func(v,_now(v,0),ppt)+v->bias[0]; 4798e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels long firstentry=0; 4808e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 4818e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels if(!(i&0xff))spinnit("centering... ",v->points-i); 4828e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 4838e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels for(j=0;j<v->entries;j++){ 4848e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels float thismetric=v->metric_func(v,_now(v,j),ppt)+v->bias[j]; 4858e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels if(thismetric<firstmetric){ 4868e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels firstmetric=thismetric; 4878e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels firstentry=j; 4888e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels } 4898e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels } 4908e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 4918e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels j=firstentry; 4928e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 4938e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels#ifdef NOISY 4948e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels fprintf(cells,"%g %g\n%g %g\n\n", 4958e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels _now(v,j)[0],_now(v,j)[1], 4968e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels ppt[0],ppt[1]); 4978e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels#endif 4988e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 4998e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels firstmetric-=v->bias[j]; 5008e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels meterror+=firstmetric; 5018e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 5028e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels if(v->centroid==0){ 5038e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels /* set up midpoints for next iter */ 5048e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels if(v->assigned[j]++){ 5058e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels for(k=0;k<v->elements;k++) 5068e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels vN(new,j)[k]+=ppt[k]; 5078e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels if(firstmetric>v->max[j])v->max[j]=firstmetric; 5088e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels }else{ 5098e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels for(k=0;k<v->elements;k++) 5108e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels vN(new,j)[k]=ppt[k]; 5118e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels v->max[j]=firstmetric; 5128e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels } 5138e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels }else{ 5148e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels /* centroid */ 5158e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels if(v->assigned[j]++){ 5168e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels for(k=0;k<v->elements;k++){ 5178e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels if(vN(new,j)[k]>ppt[k])vN(new,j)[k]=ppt[k]; 5188e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels if(vN(new2,j)[k]<ppt[k])vN(new2,j)[k]=ppt[k]; 5198e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels } 5208e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels if(firstmetric>v->max[firstentry])v->max[j]=firstmetric; 5218e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels }else{ 5228e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels for(k=0;k<v->elements;k++){ 5238e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels vN(new,j)[k]=ppt[k]; 5248e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels vN(new2,j)[k]=ppt[k]; 5258e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels } 5268e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels v->max[firstentry]=firstmetric; 5278e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels } 5288e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels } 5298e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels } 5308e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 5318e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels /* assign midpoints */ 5328e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 5338e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels for(j=0;j<v->entries;j++){ 5348e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels#ifdef NOISY 5358e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels fprintf(assig,"%ld\n",v->assigned[j]); 5368e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels fprintf(bias,"%g\n",v->bias[j]); 5378e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels#endif 5388e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels asserror+=fabs(v->assigned[j]-fdesired); 5398e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels if(v->assigned[j]){ 5408e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels if(v->centroid==0){ 5418e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels for(k=0;k<v->elements;k++) 5428e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels _now(v,j)[k]=vN(new,j)[k]/v->assigned[j]; 5438e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels }else{ 5448e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels for(k=0;k<v->elements;k++) 5458e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels _now(v,j)[k]=(vN(new,j)[k]+vN(new2,j)[k])/2.f; 5468e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels } 5478e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels } 5488e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels } 5498e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 5508e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels asserror/=(v->entries*fdesired); 5518e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 5528e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels fprintf(stderr,"Pass #%d... ",v->it); 5538e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels fprintf(stderr,": dist %g(%g) metric error=%g \n", 5548e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels asserror,fdesired,meterror/v->points); 5558e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels v->it++; 5568e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 5578e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels free(new); 5588e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels free(nearcount); 5598e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels free(nearbias); 5608e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels#ifdef NOISY 5618e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels fclose(assig); 5628e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels fclose(bias); 5638e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels fclose(cells); 5648e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels#endif 5658e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels return(asserror); 5668e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels} 5678e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels 568