IT++ Logo

vqtrain.cpp

Go to the documentation of this file.
00001 
00030 #include <itpp/srccode/vqtrain.h>
00031 #include <itpp/base/matfunc.h>
00032 #include <itpp/base/random.h>
00033 #include <itpp/base/sort.h>
00034 #include <itpp/base/specmat.h>
00035 #include <itpp/base/math/min_max.h>
00036 #include <itpp/stat/misc_stat.h>
00037 #include <iostream>
00038 
00040 
00041 namespace itpp {
00042 
00043   // the cols contains codevectors
00044   double kmeansiter(Array<vec> &DB, mat &codebook)
00045   {
00046     int                         DIM=DB(0).length(),SIZE=codebook.cols(),T=DB.length();
00047     vec                         x,xnum(SIZE);
00048     mat                         xsum(DIM,SIZE);
00049     int                         n,MinIndex,i,j,k;
00050     double                      MinS,S,D,Dold,*xp,*cp;
00051 
00052     xsum.clear();
00053     xnum.clear();
00054 
00055     n=0;
00056     D=1E20;
00057     Dold=D;
00058     D=0;
00059     for (k=0;k<T;k++) {
00060       x=DB(k);
00061       xp=x._data();
00062       MinS=1E20;
00063       MinIndex=0;
00064       for (i=0;i<SIZE;i++) {
00065         cp=&codebook(0,i);
00066         S=sqr(xp[0]-cp[0]);
00067         for (j=1;j<DIM;j++) {
00068           S+=sqr(xp[j]-cp[j]);
00069           if (S>=MinS) goto sune;
00070         }
00071         MinS=S;
00072         MinIndex=i;
00073       sune:                     i=i;
00074       }
00075       D+=MinS;
00076       cp=&xsum(0,MinIndex);
00077       for (j=0;j<DIM;j++) {
00078         cp[j]+=xp[j];
00079       }
00080       xnum(MinIndex)++;
00081     }
00082     for (i=0;i<SIZE;i++) {
00083       for (j=0;j<DIM;j++) {
00084         codebook(j,i)=xsum(j,i)/xnum(i);
00085       }
00086     }
00087     return D;
00088   }
00089 
00090   mat kmeans(Array<vec> &DB, int SIZE, int NOITER, bool VERBOSE)
00091   {
00092     int                         DIM=DB(0).length(),T=DB.length();
00093     mat                         codebook(DIM,SIZE);
00094     int                         n,i,j;
00095     double                      D,Dold;
00096     ivec                        ind(SIZE);
00097 
00098     for (i=0;i<SIZE;i++) {
00099       ind(i)=randi(0,T-1);
00100       j=0;
00101       while (j<i) {
00102         if (ind(j)==ind(i)) {
00103           ind(i)=randi(0,T-1);
00104           j=0;
00105         }
00106         j++;
00107       }
00108       codebook.set_col(i,DB(ind(i)));
00109     }
00110 
00111 
00112     if (VERBOSE) std::cout << "Training VQ..." << std::endl ;
00113 
00114     D=1E20;
00115     for (n=0;n<NOITER;n++) {
00116       Dold=D;
00117       D=kmeansiter(DB,codebook);
00118       if (VERBOSE) std::cout << n << ": " << D/T << " ";
00119       if (std::abs((D-Dold)/D) < 1e-4) break;
00120     }
00121     return codebook;
00122   }
00123 
00124   mat lbg(Array<vec> &DB, int SIZE, int NOITER, bool VERBOSE)
00125   {
00126     int         S=1,DIM=DB(0).length(),T=DB.length(),i,n;
00127     mat         cb;
00128     vec         delta=0.001*randn(DIM),x;
00129     double      D,Dold;
00130 
00131     x=zeros(DIM);
00132     for (i=0;i<T;i++) {
00133       x+=DB(i);
00134     }
00135     x/=T;
00136     cb.set_size(DIM,1);
00137     cb.set_col(0,x);
00138     while (cb.cols()<SIZE) {
00139       S=cb.cols();
00140       if (2*S<=SIZE) cb.set_size(DIM,2*S,true);
00141       else cb.set_size(DIM,2*S,true);
00142       for (i=S;i<cb.cols();i++) {
00143         cb.set_col(i,cb.get_col(i-S)+delta);
00144       }
00145       D=1E20;
00146       for (n=0;n<NOITER;n++) {
00147         Dold=D;
00148         D=kmeansiter(DB,cb);
00149         if (VERBOSE) std::cout << n << ": " << D/T << " ";
00150         if (std::abs((D-Dold)/D) < 1e-4) break;
00151       }
00152     }
00153 
00154     return cb;
00155   }
00156 
00157   mat vqtrain(Array<vec> &DB, int SIZE, int NOITER, double STARTSTEP, bool VERBOSE)
00158   {
00159     int                         DIM=DB(0).length();
00160     vec                         x;
00161     vec                         codebook(DIM*SIZE);
00162     int                         n,MinIndex,i,j;
00163     double                      MinS,S,D,step,*xp,*cp;
00164 
00165     for (i=0;i<SIZE;i++) {
00166       codebook.replace_mid(i*DIM,DB(randi(0,DB.length()-1)));
00167     }
00168     if (VERBOSE) std::cout << "Training VQ..." << std::endl ;
00169 
00170   res: D=0;
00171     for (n=0;n<NOITER;n++) {
00172       step=STARTSTEP*(1.0-double(n)/NOITER);if (step<0) step=0;
00173       x=DB(randi(0,DB.length()-1)); // seems unnecessary! Check it up.
00174       xp=x._data();
00175 
00176       MinS=1E20;
00177       MinIndex=0;
00178       for (i=0;i<SIZE;i++) {
00179         cp=&codebook(i*DIM);
00180         S=sqr(xp[0]-cp[0]);
00181         for (j=1;j<DIM;j++) {
00182           S+=sqr(xp[j]-cp[j]);
00183           if (S>=MinS) goto sune;
00184         }
00185         MinS=S;
00186         MinIndex=i;
00187       sune:             i=i;
00188       }
00189       D+=MinS;
00190       cp=&codebook(MinIndex*DIM);
00191       for (j=0;j<DIM;j++) {
00192         cp[j]+=step*(xp[j]-cp[j]);
00193       }
00194       if ((n%20000==0) && (n>1)) {
00195         if (VERBOSE) std::cout << n << ": " << D/20000 << " ";
00196         D=0;
00197       }
00198     }
00199 
00200     // checking training result
00201     vec dist(SIZE),num(SIZE);
00202 
00203     dist.clear();num.clear();
00204     for (n=0;n<DB.length();n++) {
00205       x=DB(n);
00206       xp=x._data();
00207       MinS=1E20;
00208       MinIndex=0;
00209       for (i=0;i<SIZE;i++) {
00210         cp=&codebook(i*DIM);
00211         S=sqr(xp[0]-cp[0]);
00212         for (j=1;j<DIM;j++) {
00213           S+=sqr(xp[j]-cp[j]);
00214           if (S>=MinS) goto sune2;
00215         }
00216         MinS=S;
00217         MinIndex=i;
00218       sune2:            i=i;
00219       }
00220       dist(MinIndex)+=MinS;
00221       num(MinIndex)+=1;
00222     }
00223     dist=10*log10(dist*dist.length()/sum(dist));
00224     if (VERBOSE) std::cout << std::endl << "Distortion contribution: " << dist << std::endl ;
00225     if (VERBOSE) std::cout << "Num spread: " << num/DB.length()*100 << " %" << std::endl << std::endl ;
00226     if (min(dist)<-30) {
00227       std::cout << "Points without entries! Retraining" << std::endl ;
00228       j=min_index(dist);
00229       i=max_index(num);
00230       codebook.replace_mid(j*DIM,codebook.mid(i*DIM,DIM));
00231       goto res;
00232     }
00233 
00234     mat cb(DIM,SIZE);
00235     for (i=0;i<SIZE;i++) {
00236       cb.set_col(i,codebook.mid(i*DIM,DIM));
00237     }
00238     return cb;
00239   }
00240 
00241   vec sqtrain(const vec &inDB, int SIZE)
00242   {
00243     vec         DB(inDB);
00244     vec         Levels,Levels_old;
00245     ivec        indexlist(SIZE+1);
00246     int         il,im,ih,i;
00247     int         SIZEDB=inDB.length();
00248     double      x;
00249 
00250     sort(DB);
00251     Levels=DB(round_i(linspace(0.01*SIZEDB,0.99*SIZEDB,SIZE)));
00252     Levels_old=zeros(SIZE);
00253 
00254     while (energy(Levels-Levels_old)>0.0001) {
00255       Levels_old=Levels;
00256       for (i=0;i<SIZE-1;i++) {
00257         x=(Levels(i)+Levels(i+1))/2;
00258         il=0;
00259         ih=SIZEDB-1;
00260         while (il < ih-1) {
00261           im = (il + ih)/2;
00262           if (x < DB(im)) ih = im;
00263           else il = im;
00264         }
00265         indexlist(i+1)=il;
00266       }
00267       indexlist(0)=-1;
00268       indexlist(SIZE)=SIZEDB-1;
00269       for (i=0;i<SIZE;i++) Levels(i)= mean( DB( indexlist(i)+1, indexlist(i+1)));
00270     }
00271     return Levels;
00272   }
00273 
00274   ivec bitalloc(const vec &variances, int nobits)
00275   {
00276     ivec        bitvec(variances.length());bitvec.clear();
00277     int         i,bits=nobits;
00278     vec         var=variances;
00279 
00280     while (bits>0) {
00281       i=max_index(var);
00282       var(i)/=4;
00283       bitvec(i)++;
00284       bits--;
00285     }
00286     return bitvec;
00287   }
00288 
00289 } // namespace itpp
00290 
SourceForge Logo

Generated on Sat Apr 19 10:59:23 2008 for IT++ by Doxygen 1.5.5