Incremental update toward first cut
[platform/upstream/libvorbis.git] / vq / vqsplit.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: build a VQ codebook and the encoding decision 'tree'
15  last mod: $Id: vqsplit.c,v 1.18 2000/02/23 09:10:13 xiphmont Exp $
16
17  ********************************************************************/
18
19 /* This code is *not* part of libvorbis.  It is used to generate
20    trained codebooks offline and then spit the results into a
21    pregenerated codebook that is compiled into libvorbis.  It is an
22    expensive (but good) algorithm.  Run it on big iron. */
23
24 /* There are so many optimizations to explore in *both* stages that
25    considering the undertaking is almost withering.  For now, we brute
26    force it all */
27
28 #include <stdlib.h>
29 #include <stdio.h>
30 #include <math.h>
31 #include <string.h>
32 #include <sys/time.h>
33
34 #include "vqgen.h"
35 #include "vqsplit.h"
36 #include "bookutil.h"
37
38 /* Codebook generation happens in two steps: 
39
40    1) Train the codebook with data collected from the encoder: We use
41    one of a few error metrics (which represent the distance between a
42    given data point and a candidate point in the training set) to
43    divide the training set up into cells representing roughly equal
44    probability of occurring. 
45
46    2) Generate the codebook and auxiliary data from the trained data set
47 */
48
49 /* Building a codebook from trained set **********************************
50
51    The codebook in raw form is technically finished once it's trained.
52    However, we want to finalize the representative codebook values for
53    each entry and generate auxiliary information to optimize encoding.
54    We generate the auxiliary coding tree using collected data,
55    probably the same data as in the original training */
56
57 /* At each recursion, the data set is split in half.  Cells with data
58    points on side A go into set A, same with set B.  The sets may
59    overlap.  If the cell overlaps the deviding line only very slightly
60    (provided parameter), we may choose to ignore the overlap in order
61    to pare the tree down */
62
63 long *isortvals;
64 int iascsort(const void *a,const void *b){
65   long av=isortvals[*((long *)a)];
66   long bv=isortvals[*((long *)b)];
67   return(av-bv);
68 }
69
70 /* grab it from vqgen.c */
71 extern double _dist(vqgen *v,double *a, double *b);
72
73 /* goes through the split, but just counts it and returns a metric*/
74 void vqsp_count(vqgen *v,long *membership,
75                 long *entryindex,long entries, 
76                 long *pointindex,long points,
77                 long *entryA,long *entryB,
78                 double *n, double c,
79                 long *entriesA,long *entriesB,long *entriesC){
80   long i,j,k;
81   long A=0,B=0,C=0;
82
83   memset(entryA,0,sizeof(long)*entries);
84   memset(entryB,0,sizeof(long)*entries);
85
86   /* Do the points belonging to this cell occur on sideA, sideB or
87      both? */
88
89   for(i=0;i<points;i++){
90     double *ppt=_point(v,pointindex[i]);
91     long   firstentry=membership[i];
92     double position=-c;
93
94     if(!entryA[firstentry] || !entryB[firstentry]){
95       for(k=0;k<v->elements;k++)
96       position+=ppt[k]*n[k];
97       if(position>0.)
98         entryA[firstentry]=1;
99       else
100         entryB[firstentry]=1;
101     }
102   }
103
104   /* The entry splitting isn't total, so that storage has to be
105      allocated for recursion.  Reuse the entryA/entryB vectors */
106   /* keep the entries in ascending order (relative to the original
107      list); we rely on that stability when ordering p/q choice */
108   for(j=0;j<entries;j++){
109     if(entryA[j] && entryB[j])C++;
110     if(entryA[j])entryA[A++]=entryindex[j];
111     if(entryB[j])entryB[B++]=entryindex[j];
112   }
113   *entriesA=A;
114   *entriesB=B;
115   *entriesC=C;
116 }
117
118 void pq_in_out(vqgen *v,double *n,double *c,double *p,double *q){
119   int k;
120   *c=0.;
121   for(k=0;k<v->elements;k++){
122     double center=(p[k]+q[k])/2.;
123     n[k]=(center-q[k])*2.;
124     *c+=center*n[k];
125   }
126 }
127
128 int lp_split(vqgen *v,codebook *b,
129              long *entryindex,long entries, 
130              long *pointindex,long points,
131              long depth, long *pointsofar){
132
133   encode_aux *t=b->c->encode_tree;
134
135   /* The encoder, regardless of book, will be using a straight
136      euclidian distance-to-point metric to determine closest point.
137      Thus we split the cells using the same (we've already trained the
138      codebook set spacing and distribution using special metrics and
139      even a midpoint division won't disturb the basic properties) */
140
141   long ret;
142   double *n=alloca(sizeof(double)*v->elements);
143   double c;
144   long *entryA=calloc(entries,sizeof(long));
145   long *entryB=calloc(entries,sizeof(long));
146   long entriesA=0;
147   long entriesB=0;
148   long entriesC=0;
149   long pointsA=0;
150   long i,j,k;
151
152   long *membership=malloc(sizeof(long)*points);
153
154   long   besti=-1;
155   long   bestj=-1;
156
157   char spinbuf[80];
158   sprintf(spinbuf,"splitting [%ld left]... ",v->points-*pointsofar);
159
160   if(depth==22 && points==9 && entries==2 && *pointsofar==252935){
161     fprintf(stderr,"HERE\n");
162
163   }
164
165   
166   /* which cells do points belong to?  Do this before n^2 best pair chooser. */
167
168   for(i=0;i<points;i++){
169     double *ppt=_point(v,pointindex[i]);
170     long   firstentry=0;
171     double firstmetric=_dist(v,_now(v,entryindex[0]),ppt);
172     
173     if(points*entries>64*1024)spinnit(spinbuf,entries);
174
175     for(j=1;j<entries;j++){
176       double thismetric=_dist(v,_now(v,entryindex[j]),ppt);
177       if(thismetric<=firstmetric){ /* Not <; on the line goes to higher number */
178         firstmetric=thismetric;
179         firstentry=j;
180       }
181     }
182     
183     membership[i]=firstentry;
184   }
185
186   /* We need to find the dividing hyperplane. find the median of each
187      axis as the centerpoint and the normal facing farthest point */
188
189   /* more than one way to do this part.  For small sets, we can brute
190      force it. */
191
192   if(entries<8 || (double)points*entries*entries<16.*1024*1024){
193     /* try every pair possibility */
194     double best=0;
195     double this;
196     for(i=0;i<entries-1;i++){
197       for(j=i+1;j<entries;j++){
198         spinnit(spinbuf,entries-i);
199         pq_in_out(v,n,&c,_now(v,entryindex[i]),_now(v,entryindex[j]));
200         vqsp_count(v,membership,
201                    entryindex,entries, 
202                    pointindex,points,
203                    entryA,entryB,
204                    n, c, 
205                    &entriesA,&entriesB,&entriesC);
206         this=(entriesA-entriesC)*(entriesB-entriesC);
207
208         /* when choosing best, we also want some form of stability to
209            make sure more branches are pared later; secondary
210            weighting isn;t needed as the entry lists are in ascending
211            order, and we always try p/q in the same sequence */
212         
213         if( (besti==-1) ||
214             (this>best) ){
215           
216           best=this;
217           besti=entryindex[i];
218           bestj=entryindex[j];
219
220         }
221       }
222     }
223   }else{
224     double *p=alloca(v->elements*sizeof(double));
225     double *q=alloca(v->elements*sizeof(double));
226     double best=0.;
227     
228     /* try COG/normal and furthest pairs */
229     /* meanpoint */
230     /* eventually, we want to select the closest entry and figure n/c
231        from p/q (because storing n/c is too large */
232     for(k=0;k<v->elements;k++){
233       spinnit(spinbuf,entries);
234       
235       p[k]=0.;
236       for(j=0;j<entries;j++)
237         p[k]+=v->entrylist[entryindex[j]*v->elements+k];
238       p[k]/=entries;
239
240     }
241
242     /* we go through the entries one by one, looking for the entry on
243        the other side closest to the point of reflection through the
244        center */
245
246     for(i=0;i<entries;i++){
247       double *ppi=_now(v,entryindex[i]);
248       double ref_best=0.;
249       double ref_j=-1;
250       double this;
251       spinnit(spinbuf,entries-i);
252       
253       for(k=0;k<v->elements;k++)
254         q[k]=2*p[k]-ppi[k];
255
256       for(j=0;j<entries;j++){
257         if(j!=i){
258           double this=_dist(v,q,_now(v,entryindex[j]));
259           if(ref_j==-1 || this<=ref_best){ /* <=, not <; very important */
260             ref_best=this;
261             ref_j=entryindex[j];
262           }
263         }
264       }
265
266       pq_in_out(v,n,&c,ppi,_now(v,ref_j));
267       vqsp_count(v,membership,
268                  entryindex,entries, 
269                  pointindex,points,
270                  entryA,entryB,
271                  n, c, 
272                  &entriesA,&entriesB,&entriesC);
273       this=(entriesA-entriesC)*(entriesB-entriesC);
274
275         /* when choosing best, we also want some form of stability to
276            make sure more branches are pared later; secondary
277            weighting isn;t needed as the entry lists are in ascending
278            order, and we always try p/q in the same sequence */
279         
280       if( (besti==-1) ||
281           (this>best) ){
282         
283         best=this;
284         besti=entryindex[i];
285         bestj=ref_j;
286
287       }
288     }
289     if(besti>bestj){
290       long temp=besti;
291       besti=bestj;
292       bestj=temp;
293     }
294
295   }
296   
297   /* find cells enclosing points */
298   /* count A/B points */
299
300   pq_in_out(v,n,&c,_now(v,besti),_now(v,bestj));
301   vqsp_count(v,membership,
302              entryindex,entries, 
303              pointindex,points,
304              entryA,entryB,
305              n, c,
306              &entriesA,&entriesB,&entriesC);
307
308   free(membership);
309
310   /* the point index is split, so we do an Order n rearrangement into
311      A first/B last and just pass it on */
312   {
313     long Aptr=0;
314     long Bptr=points-1;
315     while(Aptr<=Bptr){
316       while(Aptr<=Bptr){
317         double position=-c;
318         for(k=0;k<v->elements;k++)
319           position+=_point(v,pointindex[Aptr])[k]*n[k];
320         if(position<=0.)break; /* not in A */
321         Aptr++;
322       }
323       while(Aptr<=Bptr){
324         double position=-c;
325         for(k=0;k<v->elements;k++)
326           position+=_point(v,pointindex[Bptr])[k]*n[k];
327         if(position>0.)break; /* not in B */
328         Bptr--;
329       }
330       if(Aptr<Bptr){
331         long temp=pointindex[Aptr];
332         pointindex[Aptr]=pointindex[Bptr];
333         pointindex[Bptr]=temp;
334       }
335       pointsA=Aptr;
336     }
337   }
338
339   /*  fprintf(stderr,"split: total=%ld depth=%ld set A=%ld:%ld:%ld=B\n",
340       entries,depth,entriesA-entriesC,entriesC,entriesB-entriesC);*/
341   {
342     long thisaux=t->aux++;
343     if(t->aux>=t->alloc){
344       t->alloc*=2;
345       t->ptr0=realloc(t->ptr0,sizeof(long)*t->alloc);
346       t->ptr1=realloc(t->ptr1,sizeof(long)*t->alloc);
347       t->p=realloc(t->p,sizeof(long)*t->alloc);
348       t->q=realloc(t->q,sizeof(long)*t->alloc);
349     }
350     
351     t->p[thisaux]=besti;
352     t->q[thisaux]=bestj;
353     
354     if(entriesA==1){
355       ret=1;
356       t->ptr0[thisaux]=entryA[0];
357       *pointsofar+=pointsA;
358     }else{
359       t->ptr0[thisaux]= -t->aux;
360       ret=lp_split(v,b,entryA,entriesA,pointindex,pointsA,depth+1,pointsofar); 
361     }
362     if(entriesB==1){
363       ret++;
364       t->ptr1[thisaux]=entryB[0];
365       *pointsofar+=points-pointsA;
366     }else{
367       t->ptr1[thisaux]= -t->aux;
368       ret+=lp_split(v,b,entryB,entriesB,pointindex+pointsA,points-pointsA,
369                     depth+1,pointsofar); 
370     }
371   }
372   free(entryA);
373   free(entryB);
374   return(ret);
375 }
376
377 static int _node_eq(encode_aux *v, long a, long b){
378   long    Aptr0=v->ptr0[a];
379   long    Aptr1=v->ptr1[a];
380   long    Bptr0=v->ptr0[b];
381   long    Bptr1=v->ptr1[b];
382
383   /* the possibility of choosing the same p and q, but switched, can;t
384      happen because we always look for the best p/q in the same search
385      order and the search is stable */
386
387   if(Aptr0==Bptr0 && Aptr1==Bptr1)
388     return(1);
389
390   return(0);
391 }
392
393 void vqsp_book(vqgen *v, codebook *b, long *quantlist){
394   long i,j;
395   static_codebook *c=(static_codebook *)b->c;
396   encode_aux *t;
397
398   memset(b,0,sizeof(codebook));
399   memset(c,0,sizeof(static_codebook));
400   b->c=c;
401   t=c->encode_tree=calloc(1,sizeof(encode_aux));
402
403   /* make sure there are no duplicate entries and that every 
404      entry has points */
405
406   for(i=0;i<v->entries;){
407     /* duplicate? if so, eliminate */
408     for(j=0;j<i;j++){
409       if(_dist(v,_now(v,i),_now(v,j))==0.){
410         fprintf(stderr,"found a duplicate entry!  removing...\n");
411         v->entries--;
412         memcpy(_now(v,i),_now(v,v->entries),sizeof(double)*v->elements);
413         memcpy(quantlist+i*v->elements,quantlist+v->entries*v->elements,
414                sizeof(long)*v->elements);
415         break;
416       }
417     }
418     if(j==i)i++;
419   }
420
421   {
422     v->assigned=calloc(v->entries,sizeof(long));
423     for(i=0;i<v->points;i++){
424       double *ppt=_point(v,i);
425       double firstmetric=_dist(v,_now(v,0),ppt);
426       long   firstentry=0;
427
428       if(!(i&0xff))spinnit("checking... ",v->points-i);
429
430       for(j=0;j<v->entries;j++){
431         double thismetric=_dist(v,_now(v,j),ppt);
432         if(thismetric<firstmetric){
433           firstmetric=thismetric;
434           firstentry=j;
435         }
436       }
437       
438       v->assigned[firstentry]++;
439     }
440
441     for(j=0;j<v->entries;){
442       if(v->assigned[j]==0){
443         fprintf(stderr,"found an unused entry!  removing...\n");
444         v->entries--;
445         memcpy(_now(v,j),_now(v,v->entries),sizeof(double)*v->elements);
446         v->assigned[j]=v->assigned[v->elements];
447         memcpy(quantlist+j*v->elements,quantlist+v->entries*v->elements,
448                sizeof(long)*v->elements);
449         continue;
450       }
451       j++;
452     }
453   }
454
455   fprintf(stderr,"Building a book with %ld unique entries...\n",v->entries);
456
457   {
458     long *entryindex=malloc(v->entries*sizeof(long *));
459     long *pointindex=malloc(v->points*sizeof(long));
460     long pointssofar=0;
461       
462     for(i=0;i<v->entries;i++)entryindex[i]=i;
463     for(i=0;i<v->points;i++)pointindex[i]=i;
464
465     t->alloc=4096;
466     t->ptr0=malloc(sizeof(long)*t->alloc);
467     t->ptr1=malloc(sizeof(long)*t->alloc);
468     t->p=malloc(sizeof(long)*t->alloc);
469     t->q=malloc(sizeof(long)*t->alloc);
470     t->aux=0;
471     c->dim=v->elements;
472     c->entries=v->entries;
473     c->lengthlist=calloc(c->entries,sizeof(long));
474     
475     fprintf(stderr,"Leaves added: %d              \n",
476             lp_split(v,b,entryindex,v->entries,
477                      pointindex,v->points,0,&pointssofar));
478       
479     free(pointindex);
480     
481     fprintf(stderr,"Paring/rerouting redundant branches... ");
482     
483     /* The tree is likely big and redundant.  Pare and reroute branches */
484     {
485       int changedflag=1;
486       
487       while(changedflag){
488         changedflag=0;
489         
490         /* span the tree node by node; list unique decision nodes and
491            short circuit redundant branches */
492         
493         for(i=0;i<t->aux;){
494           int k;
495           
496           /* check list of unique decisions */
497           for(j=0;j<i;j++)
498             if(_node_eq(t,i,j))break;
499           
500           if(j<i){
501             /* a redundant entry; find all higher nodes referencing it and
502                short circuit them to the previously noted unique entry */
503             changedflag=1;
504             for(k=0;k<t->aux;k++){
505               if(t->ptr0[k]==-i)t->ptr0[k]=-j;
506               if(t->ptr1[k]==-i)t->ptr1[k]=-j;
507             }
508             
509             /* Now, we need to fill in the hole from this redundant
510                entry in the listing.  Insert the last entry in the list.
511                Fix the forward pointers to that last entry */
512             t->aux--;
513             t->ptr0[i]=t->ptr0[t->aux];
514             t->ptr1[i]=t->ptr1[t->aux];
515             t->p[i]=t->p[t->aux];
516             t->q[i]=t->q[t->aux];
517             for(k=0;k<t->aux;k++){
518               if(t->ptr0[k]==-t->aux)t->ptr0[k]=-i;
519               if(t->ptr1[k]==-t->aux)t->ptr1[k]=-i;
520             }
521             /* hole plugged */
522             
523           }else
524             i++;
525         }
526         
527         fprintf(stderr,"\rParing/rerouting redundant branches... "
528                 "%ld remaining   ",t->aux);
529       }
530       fprintf(stderr,"\n");
531     }
532   }
533   
534   /* run all training points through the decision tree to get a final
535      probability count */
536   {
537     long *probability=malloc(c->entries*sizeof(long));
538     for(i=0;i<c->entries;i++)probability[i]=1; /* trivial guard */
539     b->valuelist=v->entrylist; /* temporary for vqenc_entry; replaced later */
540
541     /* sigh.  A necessary hack */
542     for(i=0;i<t->aux;i++)t->p[i]*=c->dim;
543     for(i=0;i<t->aux;i++)t->q[i]*=c->dim;
544     
545     for(i=0;i<v->points;i++){
546       int ret=codebook_entry(b,v->pointlist+i*v->elements);
547       probability[ret]++;
548       if(!(i&0xff))spinnit("counting hits... ",v->points-i);
549     }
550     for(i=0;i<t->aux;i++)t->p[i]/=c->dim;
551     for(i=0;i<t->aux;i++)t->q[i]/=c->dim;
552
553     build_tree_from_lengths(c->entries,probability,c->lengthlist);
554     
555     free(probability);
556   }
557
558   /* Sort the entries by codeword length, short to long (eases
559      assignment and packing to do it now) */
560   {
561     long *wordlen=c->lengthlist;
562     long *index=malloc(c->entries*sizeof(long));
563     long *revindex=malloc(c->entries*sizeof(long));
564     int k;
565     for(i=0;i<c->entries;i++)index[i]=i;
566     isortvals=c->lengthlist;
567     qsort(index,c->entries,sizeof(long),iascsort);
568
569     /* rearrange storage; ptr0/1 first as it needs a reverse index */
570     /* n and c stay unchanged */
571     for(i=0;i<c->entries;i++)revindex[index[i]]=i;
572     for(i=0;i<t->aux;i++){
573       if(!(i&0x3f))spinnit("sorting... ",t->aux-i);
574
575       if(t->ptr0[i]>=0)t->ptr0[i]=revindex[t->ptr0[i]];
576       if(t->ptr1[i]>=0)t->ptr1[i]=revindex[t->ptr1[i]];
577       t->p[i]=revindex[t->p[i]];
578       t->q[i]=revindex[t->q[i]];
579     }
580     free(revindex);
581
582     /* map lengthlist and vallist with index */
583     c->lengthlist=calloc(c->entries,sizeof(long));
584     b->valuelist=malloc(sizeof(double)*c->entries*c->dim);
585     c->quantlist=malloc(sizeof(long)*c->entries*c->dim);
586     for(i=0;i<c->entries;i++){
587       long e=index[i];
588       for(k=0;k<c->dim;k++){
589         b->valuelist[i*c->dim+k]=v->entrylist[e*c->dim+k];
590         c->quantlist[i*c->dim+k]=quantlist[e*c->dim+k];
591       }
592       c->lengthlist[i]=wordlen[e];
593     }
594
595     free(wordlen);
596   }
597
598   fprintf(stderr,"Done.                            \n\n");
599 }
600