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
Generated on Fri Jun 8 00:27:18 2007 for IT++ by Doxygen 1.5.2