IT++ Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
vqtrain.cpp
Go to the documentation of this file.
1 
29 #include <itpp/srccode/vqtrain.h>
30 #include <itpp/base/matfunc.h>
31 #include <itpp/base/random.h>
32 #include <itpp/base/sort.h>
33 #include <itpp/base/specmat.h>
34 #include <itpp/base/math/min_max.h>
35 #include <itpp/stat/misc_stat.h>
36 #include <iostream>
37 
39 
40 namespace itpp
41 {
42 
43 // the cols contains codevectors
44 double kmeansiter(Array<vec> &DB, mat &codebook)
45 {
46  int DIM = DB(0).length(), SIZE = codebook.cols(), T = DB.length();
47  vec x, xnum(SIZE);
48  mat xsum(DIM, SIZE);
49  int n, MinIndex, i, j, k;
50  double MinS, S, D, Dold, *xp, *cp;
51 
52  xsum.clear();
53  xnum.clear();
54 
55  n = 0;
56  D = 1E20;
57  Dold = D;
58  D = 0;
59  for (k = 0;k < T;k++) {
60  x = DB(k);
61  xp = x._data();
62  MinS = 1E20;
63  MinIndex = 0;
64  for (i = 0;i < SIZE;i++) {
65  cp = &codebook(0, i);
66  S = sqr(xp[0] - cp[0]);
67  for (j = 1;j < DIM;j++) {
68  S += sqr(xp[j] - cp[j]);
69  if (S >= MinS) goto sune;
70  }
71  MinS = S;
72  MinIndex = i;
73  sune:
74  i = i;
75  }
76  D += MinS;
77  cp = &xsum(0, MinIndex);
78  for (j = 0;j < DIM;j++) {
79  cp[j] += xp[j];
80  }
81  xnum(MinIndex)++;
82  }
83  for (i = 0;i < SIZE;i++) {
84  for (j = 0;j < DIM;j++) {
85  codebook(j, i) = xsum(j, i) / xnum(i);
86  }
87  }
88  return D;
89 }
90 
91 mat kmeans(Array<vec> &DB, int SIZE, int NOITER, bool VERBOSE)
92 {
93  int DIM = DB(0).length(), T = DB.length();
94  mat codebook(DIM, SIZE);
95  int n, i, j;
96  double D, Dold;
97  ivec ind(SIZE);
98 
99  for (i = 0;i < SIZE;i++) {
100  ind(i) = randi(0, T - 1);
101  j = 0;
102  while (j < i) {
103  if (ind(j) == ind(i)) {
104  ind(i) = randi(0, T - 1);
105  j = 0;
106  }
107  j++;
108  }
109  codebook.set_col(i, DB(ind(i)));
110  }
111 
112 
113  if (VERBOSE) std::cout << "Training VQ..." << std::endl ;
114 
115  D = 1E20;
116  for (n = 0;n < NOITER;n++) {
117  Dold = D;
118  D = kmeansiter(DB, codebook);
119  if (VERBOSE) std::cout << n << ": " << D / T << " ";
120  if (std::abs((D - Dold) / D) < 1e-4) break;
121  }
122  return codebook;
123 }
124 
125 mat lbg(Array<vec> &DB, int SIZE, int NOITER, bool VERBOSE)
126 {
127  int S = 1, DIM = DB(0).length(), T = DB.length(), i, n;
128  mat cb;
129  vec delta = 0.001 * randn(DIM), x;
130  double D, Dold;
131 
132  x = zeros(DIM);
133  for (i = 0;i < T;i++) {
134  x += DB(i);
135  }
136  x /= T;
137  cb.set_size(DIM, 1);
138  cb.set_col(0, x);
139  while (cb.cols() < SIZE) {
140  S = cb.cols();
141  if (2*S <= SIZE) cb.set_size(DIM, 2*S, true);
142  else cb.set_size(DIM, 2*S, true);
143  for (i = S;i < cb.cols();i++) {
144  cb.set_col(i, cb.get_col(i - S) + delta);
145  }
146  D = 1E20;
147  for (n = 0;n < NOITER;n++) {
148  Dold = D;
149  D = kmeansiter(DB, cb);
150  if (VERBOSE) std::cout << n << ": " << D / T << " ";
151  if (std::abs((D - Dold) / D) < 1e-4) break;
152  }
153  }
154 
155  return cb;
156 }
157 
158 mat vqtrain(Array<vec> &DB, int SIZE, int NOITER, double STARTSTEP, bool VERBOSE)
159 {
160  int DIM = DB(0).length();
161  vec x;
162  vec codebook(DIM*SIZE);
163  int n, MinIndex, i, j;
164  double MinS, S, D, step, *xp, *cp;
165 
166  for (i = 0;i < SIZE;i++) {
167  codebook.replace_mid(i*DIM, DB(randi(0, DB.length() - 1)));
168  }
169  if (VERBOSE) std::cout << "Training VQ..." << std::endl ;
170 
171 res:
172  D = 0;
173  for (n = 0;n < NOITER;n++) {
174  step = STARTSTEP * (1.0 - double(n) / NOITER);
175  if (step < 0) step = 0;
176  x = DB(randi(0, DB.length() - 1)); // seems unnecessary! Check it up.
177  xp = x._data();
178 
179  MinS = 1E20;
180  MinIndex = 0;
181  for (i = 0;i < SIZE;i++) {
182  cp = &codebook(i * DIM);
183  S = sqr(xp[0] - cp[0]);
184  for (j = 1;j < DIM;j++) {
185  S += sqr(xp[j] - cp[j]);
186  if (S >= MinS) goto sune;
187  }
188  MinS = S;
189  MinIndex = i;
190  sune:
191  i = i;
192  }
193  D += MinS;
194  cp = &codebook(MinIndex * DIM);
195  for (j = 0;j < DIM;j++) {
196  cp[j] += step * (xp[j] - cp[j]);
197  }
198  if ((n % 20000 == 0) && (n > 1)) {
199  if (VERBOSE) std::cout << n << ": " << D / 20000 << " ";
200  D = 0;
201  }
202  }
203 
204  // checking training result
205  vec dist(SIZE), num(SIZE);
206 
207  dist.clear();
208  num.clear();
209  for (n = 0;n < DB.length();n++) {
210  x = DB(n);
211  xp = x._data();
212  MinS = 1E20;
213  MinIndex = 0;
214  for (i = 0;i < SIZE;i++) {
215  cp = &codebook(i * DIM);
216  S = sqr(xp[0] - cp[0]);
217  for (j = 1;j < DIM;j++) {
218  S += sqr(xp[j] - cp[j]);
219  if (S >= MinS) goto sune2;
220  }
221  MinS = S;
222  MinIndex = i;
223  sune2:
224  i = i;
225  }
226  dist(MinIndex) += MinS;
227  num(MinIndex) += 1;
228  }
229  dist = 10 * log10(dist * dist.length() / sum(dist));
230  if (VERBOSE) std::cout << std::endl << "Distortion contribution: " << dist << std::endl ;
231  if (VERBOSE) std::cout << "Num spread: " << num / DB.length()*100 << " %" << std::endl << std::endl ;
232  if (min(dist) < -30) {
233  std::cout << "Points without entries! Retraining" << std::endl ;
234  j = min_index(dist);
235  i = max_index(num);
236  codebook.replace_mid(j*DIM, codebook.mid(i*DIM, DIM));
237  goto res;
238  }
239 
240  mat cb(DIM, SIZE);
241  for (i = 0;i < SIZE;i++) {
242  cb.set_col(i, codebook.mid(i*DIM, DIM));
243  }
244  return cb;
245 }
246 
247 vec sqtrain(const vec &inDB, int SIZE)
248 {
249  vec DB(inDB);
250  vec Levels, Levels_old;
251  ivec indexlist(SIZE + 1);
252  int il, im, ih, i;
253  int SIZEDB = inDB.length();
254  double x;
255 
256  sort(DB);
257  Levels = DB(round_i(linspace(0.01 * SIZEDB, 0.99 * SIZEDB, SIZE)));
258  Levels_old = zeros(SIZE);
259 
260  while (energy(Levels - Levels_old) > 0.0001) {
261  Levels_old = Levels;
262  for (i = 0;i < SIZE - 1;i++) {
263  x = (Levels(i) + Levels(i + 1)) / 2;
264  il = 0;
265  ih = SIZEDB - 1;
266  while (il < ih - 1) {
267  im = (il + ih) / 2;
268  if (x < DB(im)) ih = im;
269  else il = im;
270  }
271  indexlist(i + 1) = il;
272  }
273  indexlist(0) = -1;
274  indexlist(SIZE) = SIZEDB - 1;
275  for (i = 0;i < SIZE;i++) Levels(i) = mean(DB(indexlist(i) + 1, indexlist(i + 1)));
276  }
277  return Levels;
278 }
279 
280 ivec bitalloc(const vec &variances, int nobits)
281 {
282  ivec bitvec(variances.length());
283  bitvec.clear();
284  int i, bits = nobits;
285  vec var = variances;
286 
287  while (bits > 0) {
288  i = max_index(var);
289  var(i) /= 4;
290  bitvec(i)++;
291  bits--;
292  }
293  return bitvec;
294 }
295 
296 } // namespace itpp
297 
SourceForge Logo

Generated on Fri Mar 21 2014 17:14:14 for IT++ by Doxygen 1.8.1.2