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