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