vq/latticehint.c

   1 /********************************************************************
   2  *                                                                  *
   3  * THIS FILE IS PART OF THE Ogg Vorbis SOFTWARE CODEC SOURCE CODE.  *
   4  * USE, DISTRIBUTION AND REPRODUCTION OF THIS SOURCE IS GOVERNED BY *
   5  * THE GNU PUBLIC LICENSE 2, WHICH IS INCLUDED WITH THIS SOURCE.    *
   6  * PLEASE READ THESE TERMS DISTRIBUTING.                            *
   7  *                                                                  *
   8  * THE OggSQUISH SOURCE CODE IS (C) COPYRIGHT 1994-2000             *
   9  * by Monty <monty@xiph.org> and The XIPHOPHORUS Company            *
  10  * http://www.xiph.org/                                             *
  11  *                                                                  *
  12  ********************************************************************
  13
  14  function: utility main for building thresh/pigeonhole encode hints
  15  last mod: $Id: latticehint.c,v 1.1 2000/07/19 18:10:02 xiphmont Exp $
  16
  17  ********************************************************************/
  18
  19 #include <stdlib.h>
  20 #include <stdio.h>
  21 #include <math.h>
  22 #include <string.h>
  23 #include <errno.h>
  24 #include "vorbis/codebook.h"
  25 #include "../lib/sharedbook.h"
  26 #include "bookutil.h"
  27
  28 /* The purpose of this util is to build encode hints for lattice
  29    codebooks so that brute forcing each codebook entry isn't needed.
  30    Threshhold hints are for books in which each scalar in the vector
  31    is independant (eg, residue) and pigeonhole lookups provide a
  32    minimum error fit for words where the scalars are interdependant
  33    (each affecting the fit of the next in sequence) as in an LSP
  34    sequential book (or can be used along with a sparse threshhold map,
  35    like a splitting tree that need not be trained)
  36
  37    If the input book is non-sequential, a threshhold hint is built.
  38    If the input book is sequential, a pigeonholing hist is built.
  39    If the book is sparse, a pigeonholing hint is built, possibly in addition
  40      to the threshhold hint
  41
  42    command line:
  43    latticehint book.vqh
  44
  45    latticehint produces book.vqh on stdout */
  46
  47 static int longsort(const void *a, const void *b){
  48   return(**((long **)a)-**((long **)b));
  49 }
  50
  51 static int addtosearch(int entry,long **tempstack,long *tempcount,int add){
  52   long *ptr=tempstack[entry];
  53   long i=tempcount[entry];
  54
  55   if(ptr){
  56     while(i--)
  57       if(*ptr++==add)return(0);
  58     tempstack[entry]=realloc(tempstack[entry],
  59                              (tempcount[entry]+1)*sizeof(long));
  60   }else{
  61     tempstack[entry]=malloc(sizeof(long));
  62   }
  63
  64   tempstack[entry][tempcount[entry]++]=add;
  65   return(1);
  66 }
  67
  68 static void setvals(int dim,encode_aux_pigeonhole *p,
  69                     long *temptrack,double *tempmin,double *tempmax,
  70                     int seqp){
  71   int i;
  72   double last=0.;
  73   for(i=0;i<dim;i++){
  74     tempmin[i]=(temptrack[i])*p->del+p->min+last;
  75     tempmax[i]=tempmin[i]+p->del;
  76     if(seqp)last=tempmin[i];
  77   }
  78 }
  79
  80 /* note that things are currently set up such that input fits that
  81    quantize outside the pigeonmap are dropped and brute-forced.  So we
  82    can ignore the <0 and >=n boundary cases in min/max error */
  83
  84 static double minerror(int dim,double *a,encode_aux_pigeonhole *p,
  85                        long *temptrack,double *tempmin,double *tempmax){
  86   int i;
  87   double err=0.;
  88   for(i=0;i<dim;i++){
  89     double eval=0.;
  90     if(a[i]<tempmin[i]){
  91       eval=tempmin[i]-a[i];
  92     }else if(a[i]>tempmax[i]){
  93       eval=a[i]-tempmax[i];
  94     }
  95     err+=eval*eval;
  96   }
  97   return(err);
  98 }
  99
 100 static double maxerror(int dim,double *a,encode_aux_pigeonhole *p,
 101                        long *temptrack,double *tempmin,double *tempmax){
 102   int i;
 103   double err=0.,eval;
 104   for(i=0;i<dim;i++){
 105     if(a[i]<tempmin[i]){
 106       eval=tempmax[i]-a[i];
 107     }else if(a[i]>tempmax[i]){
 108       eval=a[i]-tempmin[i];
 109     }else{
 110       double t1=a[i]-tempmin[i];
 111       eval=tempmax[i]-a[i];
 112       if(t1>eval)eval=t1;
 113     }
 114     err+=eval*eval;
 115   }
 116   return(err);
 117 }
 118
 119 int main(int argc,char *argv[]){
 120   codebook *b;
 121   static_codebook *c;
 122   int entries=-1,dim=-1;
 123   double min,del,cutoff=1.;
 124   char *name;
 125   long i,j;
 126
 127   if(argv[1]==NULL){
 128     fprintf(stderr,"Need a lattice book on the command line.\n");
 129     exit(1);
 130   }
 131
 132   if(argv[2])cutoff=atof(argv[2]);
 133
 134   {
 135     char *ptr;
 136     char *filename=strdup(argv[1]);
 137
 138     b=codebook_load(filename);
 139     c=(static_codebook *)(b->c);
 140
 141     ptr=strrchr(filename,'.');
 142     if(ptr){
 143       *ptr='\0';
 144       name=strdup(filename);
 145     }else{
 146       name=strdup(filename);
 147     }
 148   }
 149
 150   if(c->maptype!=1){
 151     fprintf(stderr,"Provided book is not a latticebook.\n");
 152     exit(1);
 153   }
 154
 155   entries=b->entries;
 156   dim=b->dim;
 157   min=_float32_unpack(c->q_min);
 158   del=_float32_unpack(c->q_delta);
 159
 160   /* Do we want to gen a threshold hint? */
 161   if(c->q_sequencep==0){
 162     /* yes. Discard any preexisting threshhold hint */
 163     long quantvals=_book_maptype1_quantvals(c);
 164     long **quantsort=alloca(quantvals*sizeof(long *));
 165     encode_aux_threshmatch *t=calloc(1,sizeof(encode_aux_threshmatch));
 166     c->thresh_tree=t;
 167
 168     fprintf(stderr,"Adding threshold hint to %s...\n",name);
 169
 170     /* simplest possible threshold hint only */
 171     t->quantthresh=calloc(quantvals-1,sizeof(double));
 172     t->quantmap=calloc(quantvals,sizeof(int));
 173     t->threshvals=quantvals;
 174     t->quantvals=quantvals;
 175
 176     /* the quantvals may not be in order; sort em first */
 177     for(i=0;i<quantvals;i++)quantsort[i]=c->quantlist+i;
 178     qsort(quantsort,quantvals,sizeof(long *),longsort);
 179
 180     /* ok, gen the map and thresholds */
 181     for(i=0;i<quantvals;i++)t->quantmap[i]=quantsort[i]-c->quantlist;
 182     for(i=0;i<quantvals-1;i++)
 183       t->quantthresh[i]=(*(quantsort[i])+*(quantsort[i+1]))*.5*del+min;
 184   }
 185
 186   /* Do we want to gen a pigeonhole hint? */
 187   for(i=0;i<entries;i++)if(c->lengthlist[i]==0)break;
 188   if(c->q_sequencep || i<entries){
 189     long **tempstack;
 190     long *tempcount;
 191     long *temptrack;
 192     double *tempmin;
 193     double *tempmax;
 194     long totalstack=0;
 195     long subpigeons;
 196     long quantvals=_book_maptype1_quantvals(c);
 197     int changep=1;
 198
 199     encode_aux_pigeonhole *p=calloc(1,sizeof(encode_aux_pigeonhole));
 200     c->pigeon_tree=p;
 201
 202     fprintf(stderr,"Adding pigeonhole hint to %s...\n",name);
 203
 204     /* the idea is that we quantize uniformly, even in a nonuniform
 205        lattice, so that quantization of one scalar has a predictable
 206        result on the next sequential scalar in a greedy matching
 207        algorithm.  We generate a lookup based on the quantization of
 208        the vector (pigeonmap groups quantized entries together) and
 209        list the entries that could possible be the best fit for any
 210        given member of that pigeonhole.  The encode process then has a
 211        much smaller list to brute force */
 212
 213     /* find our pigeonhole-specific quantization values, fill in the
 214        quant value->pigeonhole map */
 215     p->del=del;
 216     p->min=min-del*.5;
 217     p->quantvals=quantvals;
 218     {
 219       int max=0;
 220       for(i=0;i<quantvals;i++)if(max<c->quantlist[i])max=c->quantlist[i];
 221       p->mapentries=max;
 222     }
 223     p->pigeonmap=malloc(p->mapentries*sizeof(long));
 224
 225     /* pigeonhole roughly on the boundaries of the quantvals; the
 226        exact pigeonhole grouping is an optimization issue, not a
 227        correctness issue */
 228     for(i=0;i<p->mapentries;i++){
 229       double thisval=del*(i+.5)+min; /* middle of the quant zone */
 230       int quant=0;
 231       double err=fabs(c->quantlist[0]*del+min-thisval);
 232       for(j=1;j<quantvals;j++){
 233         double thiserr=fabs(c->quantlist[j]*del+min-thisval);
 234         if(thiserr<err){
 235           quant=j;
 236           err=thiserr;
 237         }
 238       }
 239       p->pigeonmap[i]=quant;
 240     }
 241
 242     /* pigeonmap complete.  Now do the grungy business of finding the
 243     entries that could possibly be the best fit for a value appearing
 244     in the pigeonhole. The trick that allows the below to work is the
 245     uniform quantization; even though the scalars may be 'sequential'
 246     (each a delta from the last), the uniform quantization means that
 247     the error variance is *not* dependant.  Given a pigeonhole and an
 248     entry, we can find the minimum and maximum possible errors
 249     (relative to the entry) for any point that could appear in the
 250     pigeonhole */
 251
 252     /* must iterate over both pigeonholes and entries */
 253     /* temporarily (in order to avoid thinking hard), we grow each
 254        pigeonhole seperately, the build a stack of 'em later */
 255     subpigeons=1;
 256     for(i=0;i<dim;i++)subpigeons*=p->mapentries;
 257     temptrack=calloc(dim,sizeof(long));
 258     tempmin=calloc(dim,sizeof(double));
 259     tempmax=calloc(dim,sizeof(double));
 260     tempstack=calloc(entries,sizeof(long *));
 261     tempcount=calloc(entries,sizeof(long));
 262
 263     while(1){
 264       double errorpost=-1;
 265       char buffer[80];
 266
 267       /* map our current pigeonhole to a 'big pigeonhole' so we know
 268          what list we're after */
 269       int entry=0;
 270       for(i=dim-1;i>=0;i--)entry=entry*quantvals+p->pigeonmap[temptrack[i]];
 271       setvals(dim,p,temptrack,tempmin,tempmax,c->q_sequencep);
 272       sprintf(buffer,"Building pigeonhole search list [%ld]...",totalstack);
 273
 274
 275       /* Search all entries to find the one with the minimum possible
 276          maximum error.  Record that error */
 277       for(i=0;i<entries;i++){
 278         if(c->lengthlist[i]>0){
 279           double this=maxerror(dim,b->valuelist+i*dim,p,
 280                                temptrack,tempmin,tempmax);
 281           if(errorpost==-1 || this<errorpost)errorpost=this;
 282           spinnit(buffer,subpigeons);
 283         }
 284       }
 285
 286       /* Our search list will contain all entries with a minimum
 287          possible error <= our errorpost */
 288       for(i=0;i<entries;i++)
 289         if(c->lengthlist[i]>0){
 290           spinnit(buffer,subpigeons);
 291           if(minerror(dim,b->valuelist+i*dim,p,
 292                       temptrack,tempmin,tempmax)<errorpost)
 293             totalstack+=addtosearch(entry,tempstack,tempcount,i);
 294         }
 295
 296       for(i=0;i<dim;i++){
 297         temptrack[i]++;
 298         if(temptrack[i]<p->mapentries)break;
 299         temptrack[i]=0;
 300       }
 301       if(i==dim)break;
 302       subpigeons--;
 303     }
 304
 305     fprintf(stderr,"\r                                                     "
 306             "\rTotal search list size (all entries): %ld\n",totalstack);
 307
 308     /* pare the index of lists for improbable quantizations (where
 309        improbable is determined by c->lengthlist; we assume that
 310        pigeonholing is in sync with the codeword cells, which it is */
 311     for(i=0;i<entries;i++){
 312       double probability= 1./(1<<c->lengthlist[i]);
 313       if(c->lengthlist[i]==0 || probability*entries<cutoff){
 314         totalstack-=tempcount[i];
 315         tempcount[i]=0;
 316       }
 317     }
 318
 319     /* pare the list of shortlists; merge contained and similar lists
 320        together */
 321     p->fitmap=malloc(entries*sizeof(long));
 322     for(i=0;i<entries;i++)p->fitmap[i]=-1;
 323     while(changep){
 324       char buffer[80];
 325       changep=0;
 326
 327       for(i=0;i<entries;i++){
 328         if(p->fitmap[i]<0 && tempcount[i]){
 329           for(j=i+1;j<entries;j++){
 330             if(p->fitmap[j]<0 && tempcount[j]){
 331               /* is one list a superset, or are they sufficiently similar? */
 332               int amiss=0,bmiss=0,ii,jj;
 333               for(ii=0;ii<tempcount[i];ii++){
 334                 for(jj=0;jj<tempcount[j];jj++)
 335                   if(tempstack[i][ii]==tempstack[j][jj])break;
 336                 if(jj==tempcount[j])amiss++;
 337               }
 338               for(jj=0;jj<tempcount[j];jj++){
 339                 for(ii=0;ii<tempcount[i];ii++)
 340                   if(tempstack[i][ii]==tempstack[j][jj])break;
 341                 if(ii==tempcount[i])bmiss++;
 342               }
 343               if(amiss==0 ||
 344                  bmiss==0 ||
 345                  (amiss*4<tempcount[i] && bmiss*4<tempcount[j] &&
 346                   tempcount[i]+bmiss<entries/30)){
 347
 348                 /*superset/similar  Add all of one to the other. */
 349                 for(jj=0;jj<tempcount[j];jj++)
 350                   totalstack+=addtosearch(i,tempstack,tempcount,
 351                                           tempstack[j][jj]);
 352                 totalstack-=tempcount[j];
 353                 p->fitmap[j]=i;
 354                 changep=1;
 355               }
 356             }
 357           }
 358           sprintf(buffer,"Consolidating [%ld total, %s]... ",totalstack,
 359                   changep?"reit":"nochange");
 360           spinnit(buffer,entries-i);
 361         }
 362       }
 363     }
 364
 365     /* repack the temp stack in final form */
 366
 367     p->fittotal=totalstack;
 368     p->fitlist=malloc((totalstack+1)*sizeof(long));
 369     p->fitlength=malloc(entries*sizeof(long));
 370     {
 371       long usage=0;
 372       for(i=0;i<entries;i++){
 373         if(p->fitmap[i]==-1){
 374           if(tempcount[i])
 375             memcpy(p->fitlist+usage,tempstack[i],tempcount[i]*sizeof(long));
 376           p->fitmap[i]=usage;
 377           p->fitlength[i]=tempcount[i];
 378           usage+=tempcount[i];
 379           if(usage>totalstack){
 380             fprintf(stderr,"Internal error; usage>totalstack\n");
 381             exit(1);
 382           }
 383         }else{
 384           p->fitlength[i]=p->fitlength[p->fitmap[i]];
 385           p->fitmap[i]=p->fitmap[p->fitmap[i]];
 386         }
 387       }
 388     }
 389   }
 390
 391   write_codebook(stdout,name,c);
 392   fprintf(stderr,"\r                                                     "
 393           "\nDone.\n");
 394   exit(0);
 395 }