IT++ Logo Newcom Logo

vqtrain.cpp

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

Generated on Wed Apr 18 11:45:35 2007 for IT++ by Doxygen 1.5.2