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 
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 
SourceForge Logo

Generated on Thu Apr 23 20:04:05 2009 for IT++ by Doxygen 1.5.8