IT++ Logo
vqtrain.cpp
Go to the documentation of this file.
00001 
00029 #include <itpp/srccode/vqtrain.h>
00030 #include <itpp/base/matfunc.h>
00031 #include <itpp/base/random.h>
00032 #include <itpp/base/sort.h>
00033 #include <itpp/base/specmat.h>
00034 #include <itpp/base/math/min_max.h>
00035 #include <itpp/stat/misc_stat.h>
00036 #include <iostream>
00037 
00039 
00040 namespace itpp
00041 {
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:
00074       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:
00172   D = 0;
00173   for (n = 0;n < NOITER;n++) {
00174     step = STARTSTEP * (1.0 - double(n) / NOITER);
00175     if (step < 0) step = 0;
00176     x = DB(randi(0, DB.length() - 1)); // seems unnecessary! Check it up.
00177     xp = x._data();
00178 
00179     MinS = 1E20;
00180     MinIndex = 0;
00181     for (i = 0;i < SIZE;i++) {
00182       cp = &codebook(i * DIM);
00183       S = sqr(xp[0] - cp[0]);
00184       for (j = 1;j < DIM;j++) {
00185         S += sqr(xp[j] - cp[j]);
00186         if (S >= MinS) goto sune;
00187       }
00188       MinS = S;
00189       MinIndex = i;
00190     sune:
00191       i = i;
00192     }
00193     D += MinS;
00194     cp = &codebook(MinIndex * DIM);
00195     for (j = 0;j < DIM;j++) {
00196       cp[j] += step * (xp[j] - cp[j]);
00197     }
00198     if ((n % 20000 == 0) && (n > 1)) {
00199       if (VERBOSE) std::cout << n << ": " << D / 20000 << " ";
00200       D = 0;
00201     }
00202   }
00203 
00204   // checking training result
00205   vec dist(SIZE), num(SIZE);
00206 
00207   dist.clear();
00208   num.clear();
00209   for (n = 0;n < DB.length();n++) {
00210     x = DB(n);
00211     xp = x._data();
00212     MinS = 1E20;
00213     MinIndex = 0;
00214     for (i = 0;i < SIZE;i++) {
00215       cp = &codebook(i * DIM);
00216       S = sqr(xp[0] - cp[0]);
00217       for (j = 1;j < DIM;j++) {
00218         S += sqr(xp[j] - cp[j]);
00219         if (S >= MinS) goto sune2;
00220       }
00221       MinS = S;
00222       MinIndex = i;
00223     sune2:
00224       i = i;
00225     }
00226     dist(MinIndex) += MinS;
00227     num(MinIndex) += 1;
00228   }
00229   dist = 10 * log10(dist * dist.length() / sum(dist));
00230   if (VERBOSE) std::cout << std::endl << "Distortion contribution: " << dist << std::endl ;
00231   if (VERBOSE) std::cout << "Num spread: " << num / DB.length()*100 << " %" << std::endl << std::endl ;
00232   if (min(dist) < -30) {
00233     std::cout << "Points without entries! Retraining" << std::endl ;
00234     j = min_index(dist);
00235     i = max_index(num);
00236     codebook.replace_mid(j*DIM, codebook.mid(i*DIM, DIM));
00237     goto res;
00238   }
00239 
00240   mat cb(DIM, SIZE);
00241   for (i = 0;i < SIZE;i++) {
00242     cb.set_col(i, codebook.mid(i*DIM, DIM));
00243   }
00244   return cb;
00245 }
00246 
00247 vec sqtrain(const vec &inDB, int SIZE)
00248 {
00249   vec  DB(inDB);
00250   vec  Levels, Levels_old;
00251   ivec indexlist(SIZE + 1);
00252   int  il, im, ih, i;
00253   int  SIZEDB = inDB.length();
00254   double x;
00255 
00256   sort(DB);
00257   Levels = DB(round_i(linspace(0.01 * SIZEDB, 0.99 * SIZEDB, SIZE)));
00258   Levels_old = zeros(SIZE);
00259 
00260   while (energy(Levels - Levels_old) > 0.0001) {
00261     Levels_old = Levels;
00262     for (i = 0;i < SIZE - 1;i++) {
00263       x = (Levels(i) + Levels(i + 1)) / 2;
00264       il = 0;
00265       ih = SIZEDB - 1;
00266       while (il < ih - 1) {
00267         im = (il + ih) / 2;
00268         if (x < DB(im)) ih = im;
00269         else il = im;
00270       }
00271       indexlist(i + 1) = il;
00272     }
00273     indexlist(0) = -1;
00274     indexlist(SIZE) = SIZEDB - 1;
00275     for (i = 0;i < SIZE;i++) Levels(i) = mean(DB(indexlist(i) + 1, indexlist(i + 1)));
00276   }
00277   return Levels;
00278 }
00279 
00280 ivec bitalloc(const vec &variances, int nobits)
00281 {
00282   ivec bitvec(variances.length());
00283   bitvec.clear();
00284   int  i, bits = nobits;
00285   vec  var = variances;
00286 
00287   while (bits > 0) {
00288     i = max_index(var);
00289     var(i) /= 4;
00290     bitvec(i)++;
00291     bits--;
00292   }
00293   return bitvec;
00294 }
00295 
00296 } // namespace itpp
00297 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
SourceForge Logo

Generated on Sat Jul 9 2011 15:21:33 for IT++ by Doxygen 1.7.4