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
Generated on Sat Apr 19 10:41:57 2008 for IT++ by Doxygen 1.5.5