IT++ Logo

mog_generic.cpp

Go to the documentation of this file.
00001 
00030 #include <itpp/stat/mog_generic.h>
00031 #include <itpp/base/algebra/inv.h>
00032 #include <itpp/base/algebra/det.h>
00033 #include <itpp/base/matfunc.h>
00034 #include <itpp/base/itfile.h>
00035 #include <itpp/base/math/misc.h>
00036 #include <itpp/base/math/log_exp.h>
00037 
00038 
00039 namespace itpp
00040 {
00041 
00042 
00043 void MOG_generic::init() { cleanup(); }
00044 
00045 
00046 void MOG_generic::init(const int &K_in, const int &D_in, bool full_in)
00047 {
00048   valid = false;
00049 
00050   it_assert(K_in >= 0, "MOG_generic::init(): number of Gaussians must be greater than zero");
00051   it_assert(D_in >= 0, "MOG_generic::init(): dimensionality must be greater than zero");
00052 
00053   K = K_in;
00054   D = D_in;
00055   full = full_in;
00056 
00057   set_means_zero_internal();
00058   full ? set_full_covs_unity_internal() : set_diag_covs_unity_internal();
00059   set_weights_uniform_internal();
00060   setup_misc();
00061 
00062   valid = true;
00063   do_checks = true;
00064   paranoid = false;
00065 
00066 }
00067 
00068 
00069 void MOG_generic::init(Array<vec> &means_in, bool full_in)
00070 {
00071   valid = false;
00072 
00073   K = means_in.size();
00074   D = means_in(0).size();
00075   full = full_in;
00076 
00077   it_assert(check_array_uniformity(means_in), "MOG_generic::init(): 'means' is empty or contains vectors of varying dimensionality");
00078   set_means(means_in);
00079   full ? set_full_covs_unity_internal() : set_diag_covs_unity_internal();
00080   set_weights_uniform_internal();
00081   setup_misc();
00082 
00083   valid = true;
00084   do_checks = true;
00085   paranoid = false;
00086 }
00087 
00088 
00089 void MOG_generic::init(Array<vec> &means_in, Array<vec> &diag_covs_in, vec &weights_in)
00090 {
00091   valid = false;
00092 
00093   K = means_in.size();
00094   D = means_in(0).size();
00095   full = false;
00096 
00097   it_assert(check_array_uniformity(means_in), "MOG_generic::init(): 'means' is empty or contains vectors of varying dimensionality");
00098 
00099   set_means_internal(means_in);
00100   set_diag_covs_internal(diag_covs_in);
00101   set_weights_internal(weights_in);
00102   setup_misc();
00103 
00104   valid = true;
00105   do_checks = true;
00106   paranoid = false;
00107 }
00108 
00109 
00110 void MOG_generic::init(Array<vec> &means_in, Array<mat> &full_covs_in, vec &weights_in)
00111 {
00112   valid = false;
00113 
00114   K = means_in.size();
00115   D = means_in(0).size();
00116   full = true;
00117 
00118   it_assert(check_array_uniformity(means_in), "MOG_generic::init(): 'means' is empty or contains vectors of varying dimensionality");
00119   set_means_internal(means_in);
00120   set_full_covs_internal(full_covs_in);
00121   set_weights_internal(weights_in);
00122   setup_misc();
00123 
00124   valid = true;
00125   do_checks = true;
00126   paranoid = false;
00127 }
00128 
00129 
00130 bool MOG_generic::check_size(const vec &x_in) const
00131 {
00132   if (x_in.size() == D) return true;
00133   return false;
00134 }
00135 
00136 
00137 bool MOG_generic::check_size(const Array<vec> &X_in) const
00138 {
00139   if (check_array_uniformity(X_in)) return check_size(X_in(0));
00140   return false;
00141 }
00142 
00143 
00144 bool MOG_generic::check_array_uniformity(const Array<vec> & A) const
00145 {
00146   int rows = A.size();
00147   int cols = A(0).size();
00148 
00149   if (!rows || !cols) return false;
00150 
00151   for (int row = 1; row < rows; row++)  if (A(row).size() != cols)  return false;
00152   return true;
00153 }
00154 
00155 
00156 void MOG_generic::set_means_zero_internal()
00157 {
00158   means.set_size(K);
00159   for (int k = 0; k < K; k++) { means(k).set_size(D);  means(k) = 0.0; }
00160   setup_means();
00161 }
00162 
00163 
00164 void MOG_generic::set_means_internal(Array<vec> &means_in)
00165 {
00166   it_assert((means_in.size() == K), "MOG_generic::set_means_internal(): number of vectors in 'means' is not equivalent to number of Gaussians");
00167 
00168   for (int k = 0; k < K; k++)
00169     it_assert((means_in(k).size() == D), "MOG_generic::set_means_internal(): dimensionality mismatch between model and one or more vectors in 'means'");
00170 
00171   for (int k = 0; k < K; k++)
00172     for (int d = 0; d < D; d++)
00173       it_assert(std::isfinite(means_in(k)(d)), "MOG_generic::set_means_internal(): 'means' has a non-finite value");
00174 
00175   means = means_in;
00176   setup_means();
00177 }
00178 
00179 
00180 void MOG_generic::set_diag_covs_internal(Array<vec> &diag_covs_in)
00181 {
00182   it_assert((diag_covs_in.size() == K), "MOG_generic::set_diag_covs_internal(): number of vectors in 'diag_covs' does not match number of Gaussians");
00183 
00184   for (int k = 0; k < K; k++)
00185     it_assert((diag_covs_in(k).size() == D), "MOG_generic::set_diag_covs_internal(): dimensionality mismatch between model and one or more vectors in 'diag_covs'");
00186 
00187   for (int k = 0; k < K; k++)
00188     for (int d = 0; d < D; d++) {
00189       it_assert((diag_covs_in(k)(d) > 0.0), "MOG_generic::set_diag_covs_internal(): 'diag_covs' has a zero or negative value");
00190       it_assert(std::isfinite(diag_covs_in(k)(d)), "MOG_generic::set_diag_covs_internal(): 'diag_covs' has a non-finite value");
00191     }
00192 
00193   full_covs.set_size(0);
00194   diag_covs = diag_covs_in;
00195   full = false;
00196   setup_covs();
00197 }
00198 
00199 
00200 void MOG_generic::set_full_covs_internal(Array<mat> &full_covs_in)
00201 {
00202   it_assert((full_covs_in.size() == K), "MOG_generic::set_full_covs_internal(): number of matrices in 'full_covs' does not match number of Gaussians");
00203 
00204   for (int k = 0; k < K; k++)
00205     it_assert(((full_covs_in(k).rows() == D) && (full_covs_in(k).cols() == D)), "MOG_generic::set_full_covs_internal(): dimensionality mismatch between model and one or more matrices in 'full_covs'");
00206 
00207   for (int k = 0; k < K; k++)
00208     for (int i = 0; i < D; i++)
00209       for (int j = 0; j < D; j++) {
00210         it_assert(std::isfinite(full_covs_in(k)(i, j)), "MOG_generic::set_full_covs_internal(): 'full_covs' has a non-finite value");
00211         if (i == j) {
00212           it_assert(full_covs_in(k)(i, j) > 0.0,
00213                     "MOG_generic::set_full_covs_internal(): 'full_covs' has "
00214                     "a zero or negative value on a diagonal");
00215         }
00216       }
00217 
00218   full_covs = full_covs_in;
00219   diag_covs.set_size(0);
00220   full = true;
00221   setup_covs();
00222 }
00223 
00224 
00225 void MOG_generic::set_weights_internal(vec &weights_in)
00226 {
00227 
00228   it_assert((weights_in.size() == K), "MOG_generic::set_weights_internal(): number of elements in 'weights' does not match number of Gaussians");
00229 
00230   for (int k = 0; k < K; k++) {
00231     it_assert((weights_in(k) >= 0), "MOG_generic::set_weights_internal(): 'weights' has a negative value");
00232     it_assert(std::isfinite(weights_in(k)), "MOG_generic::set_weights_internal(): 'weights' has a non-finite value");
00233   }
00234 
00235   weights = weights_in;
00236   setup_weights();
00237 
00238 }
00239 
00240 
00241 void MOG_generic::set_diag_covs_unity_internal()
00242 {
00243   full_covs.set_size(0);
00244   diag_covs.set_size(K);
00245 
00246   for (int k = 0; k < K; k++) { diag_covs(k).set_size(D);  diag_covs(k) = 1.0; }
00247 
00248   full = false;
00249   setup_covs();
00250 }
00251 
00252 
00253 void MOG_generic::set_full_covs_unity_internal()
00254 {
00255   full_covs.set_size(K);
00256   diag_covs.set_size(0);
00257 
00258   for (int k = 0; k < K; k++) {
00259     full_covs(k).set_size(D, D);
00260     full_covs(k) = 0.0;
00261     for (int d = 0;d < D;d++)  full_covs(k)(d, d) = 1.0;
00262   }
00263 
00264   full = true;
00265   setup_covs();
00266 }
00267 
00268 
00269 void MOG_generic::set_weights_uniform_internal()
00270 {
00271   weights.set_size(K);
00272   weights = 1.0 / K;
00273   setup_weights();
00274 }
00275 
00276 
00277 void MOG_generic::setup_means() { }
00278 
00279 void MOG_generic::setup_covs()
00280 {
00281 
00282   double Ddiv2_log_2pi = D / 2.0 * std::log(m_2pi);
00283   log_det_etc.set_size(K);
00284 
00285   if (full) {
00286     full_covs_inv.set_size(K);
00287     diag_covs_inv_etc.set_size(0);
00288     for (int k = 0;k < K;k++)  full_covs_inv(k) = inv(full_covs(k));
00289     for (int k = 0;k < K;k++)  log_det_etc(k) = -Ddiv2_log_2pi - 0.5 * std::log(det(full_covs(k)));
00290   }
00291   else {
00292     full_covs_inv.set_size(0);
00293     diag_covs_inv_etc.set_size(K);
00294     for (int k = 0;k < K;k++) diag_covs_inv_etc(k).set_size(D);
00295 
00296     for (int k = 0;k < K;k++) {
00297       double acc = 0.0;
00298       vec & diag_cov = diag_covs(k);
00299       vec & diag_cov_inv_etc = diag_covs_inv_etc(k);
00300 
00301       for (int d = 0;d < D;d++)  {
00302         double tmp = diag_cov(d);
00303         diag_cov_inv_etc(d) = 1.0 / (2.0 * tmp);
00304         acc += std::log(tmp);
00305       }
00306 
00307       log_det_etc(k) = -Ddiv2_log_2pi - 0.5 * acc;
00308 
00309     }
00310   }
00311 }
00312 
00313 
00314 void MOG_generic::setup_weights()
00315 {
00316   weights /= sum(weights);
00317   log_weights = log(weights);
00318 }
00319 
00320 
00321 void MOG_generic::setup_misc()
00322 {
00323   log_max_K = std::log(std::numeric_limits<double>::max() / K);
00324   tmpvecD.set_size(D);
00325   tmpvecK.set_size(K);
00326 }
00327 
00328 
00329 void MOG_generic::cleanup()
00330 {
00331 
00332   valid = false;
00333   do_checks = true;
00334   K = 0;
00335   D = 0;
00336 
00337   tmpvecD.set_size(0);
00338   tmpvecK.set_size(0);
00339   means.set_size(0);
00340   diag_covs.set_size(0);
00341   full_covs.set_size(0);
00342   weights.set_size(0);
00343   log_det_etc.set_size(0);
00344   log_weights.set_size(0);
00345   full_covs_inv.set_size(0);
00346   diag_covs_inv_etc.set_size(0);
00347 
00348 }
00349 
00350 
00351 void MOG_generic::set_means(Array<vec> &means_in)
00352 {
00353   if (!valid) return;
00354   set_means_internal(means_in);
00355 }
00356 
00357 
00358 void MOG_generic::set_means_zero()
00359 {
00360   if (!valid) return;
00361   set_means_zero_internal();
00362 }
00363 
00364 
00365 void MOG_generic::set_diag_covs(Array<vec> &diag_covs_in)
00366 {
00367   if (!valid) return;
00368   set_diag_covs_internal(diag_covs_in);
00369 }
00370 
00371 
00372 void MOG_generic::set_full_covs(Array<mat> &full_covs_in)
00373 {
00374   if (!valid) return;
00375   set_full_covs_internal(full_covs_in);
00376 }
00377 
00378 
00379 void MOG_generic::set_weights(vec &weights_in)
00380 {
00381   if (!valid) return;
00382   set_weights_internal(weights_in);
00383 }
00384 
00385 
00386 void MOG_generic::set_diag_covs_unity()
00387 {
00388   if (!valid) return;
00389   set_diag_covs_unity_internal();
00390 }
00391 
00392 
00393 void MOG_generic::set_full_covs_unity()
00394 {
00395   if (!valid) return;
00396   set_full_covs_unity_internal();
00397 }
00398 
00399 
00400 void MOG_generic::set_weights_uniform()
00401 {
00402   if (!valid) return;
00403   set_weights_uniform_internal();
00404 }
00405 
00406 
00407 void MOG_generic::load(const std::string &name_in)
00408 {
00409   valid = false;
00410 
00411   it_assert(exist(name_in), "MOG_generic::load(): couldn't access file '" + name_in + "'");
00412   it_file ff(name_in);
00413 
00414   bool contents = ff.exists("means") && (ff.exists("diag_covs") || ff.exists("full_covs")) && ff.exists("weights");
00415   it_assert(contents, "MOG_generic::load(): file '" + name_in + "' doesn't appear to be a model file");
00416 
00417   Array<vec> means_in;
00418   ff >> Name("means") >> means_in;
00419   vec weights_in;
00420   ff >> Name("weights") >> weights_in;
00421 
00422   if (ff.exists("full_covs")) {
00423     Array<mat> full_covs_in;
00424     ff >> Name("full_covs") >> full_covs_in;
00425     init(means_in, full_covs_in, weights_in);
00426   }
00427   else {
00428     Array<vec> diag_covs_in;
00429     ff >> Name("diag_covs") >> diag_covs_in;
00430     init(means_in, diag_covs_in, weights_in);
00431   }
00432 
00433   ff.close();
00434 
00435 }
00436 
00437 
00438 void MOG_generic::save(const std::string &name_in) const
00439 {
00440   if (!valid) return;
00441 
00442   it_file ff(name_in);
00443 
00444   ff << Name("means") << means;
00445   if (full) ff << Name("full_covs") << full_covs;
00446   else ff << Name("diag_covs") << diag_covs;
00447   ff << Name("weights") << weights;
00448 
00449   ff.close();
00450 
00451 }
00452 
00453 void MOG_generic::join(const MOG_generic &B_in)
00454 {
00455 
00456   if (!valid) return;
00457   if (!B_in.is_valid()) return;
00458 
00459   it_assert((full == B_in.is_full()), "MOG_generic::join(): given model must be of the same type");
00460   it_assert((B_in.get_D() == D), "MOG_generic::join(): given model has different dimensionality");
00461   it_assert((B_in.get_K() > 0),  "MOG_generic::join(): given model has no components");
00462 
00463   int new_K = K + B_in.get_K();
00464   vec new_weights(new_K);
00465   vec B_in_weights = B_in.get_weights();
00466 
00467   double alpha = double(K) / double(new_K);
00468   double beta = double(B_in.get_K()) / double(new_K);
00469 
00470   for (int k = 0;k < K;k++)  new_weights(k) = alpha * weights(k);
00471   for (int k = K;k < new_K;k++)  new_weights(k) = beta * B_in_weights(k);
00472 
00473   Array<vec> new_means = concat(means, B_in.get_means());
00474 
00475   if (full) {
00476     Array<mat> new_full_covs = concat(full_covs, B_in.get_full_covs());
00477     init(new_means, new_full_covs, new_weights);
00478   }
00479   else {
00480     Array<vec> new_diag_covs = concat(diag_covs, B_in.get_diag_covs());
00481     init(new_means, new_diag_covs, new_weights);
00482   }
00483 }
00484 
00485 
00486 void MOG_generic::convert_to_diag_internal()
00487 {
00488   if (!full) return;
00489 
00490   diag_covs.set_size(K);
00491   for (int k = 0;k < K;k++)  diag_covs(k) = diag(full_covs(k));
00492   full_covs.set_size(0);
00493 
00494   full = false;
00495   setup_covs();
00496 }
00497 
00498 
00499 void MOG_generic::convert_to_diag()
00500 {
00501   if (!valid) return;
00502   convert_to_diag_internal();
00503 }
00504 
00505 
00506 void MOG_generic::convert_to_full_internal()
00507 {
00508   if (full) return;
00509 
00510   full_covs.set_size(K);
00511   for (int k = 0;k < K;k++)  full_covs(k) = diag(diag_covs(k));
00512   diag_covs.set_size(0);
00513 
00514   full = true;
00515   setup_covs();
00516 }
00517 
00518 void MOG_generic::convert_to_full()
00519 {
00520   if (!valid) return;
00521   convert_to_full_internal();
00522 }
00523 
00524 
00525 
00526 double MOG_generic::log_lhood_single_gaus_internal(const vec &x_in, const int k)
00527 {
00528 
00529   const vec & mean = means(k);
00530 
00531   if (full) {
00532     for (int d = 0;d < D;d++)  tmpvecD[d] = x_in[d] - mean[d];
00533     double tmpval = tmpvecD * (full_covs_inv(k) * tmpvecD);
00534     return(log_det_etc[k] - 0.5*tmpval);
00535   }
00536   else {
00537     const vec & diag_cov_inv_etc = diag_covs_inv_etc(k);
00538 
00539     double acc = 0.0;
00540 
00541     for (int d = 0; d < D; d++) {
00542       double tmpval = x_in[d] - mean[d];
00543       acc += (tmpval * tmpval) * diag_cov_inv_etc[d];
00544     }
00545     return(log_det_etc[k] - acc);
00546   }
00547 
00548 }
00549 
00550 double MOG_generic::log_lhood_single_gaus(const vec &x_in, const int k)
00551 {
00552   if (do_checks) {
00553     it_assert(valid, "MOG_generic::log_lhood_single_gaus(): model not valid");
00554     it_assert(check_size(x_in), "MOG_generic::log_lhood_single_gaus(): x has wrong dimensionality");
00555     it_assert(((k >= 0) && (k < K)), "MOG_generic::log_lhood_single_gaus(): k specifies a non-existant Gaussian");
00556   }
00557   return log_lhood_single_gaus_internal(x_in, k);
00558 }
00559 
00560 
00561 double MOG_generic::log_lhood_internal(const vec &x_in)
00562 {
00563 
00564   bool danger = paranoid;
00565 
00566   for (int k = 0;k < K;k++)  {
00567     double tmp = log_weights[k] + log_lhood_single_gaus_internal(x_in, k);
00568     tmpvecK[k] = tmp;
00569 
00570     if (tmp >= log_max_K)  danger = true;
00571   }
00572 
00573   if (danger) {
00574     double log_sum = tmpvecK[0];
00575     for (int k = 1; k < K; k++)  log_sum = log_add(log_sum, tmpvecK[k]);
00576     return(log_sum);
00577   }
00578   else {
00579     double sum = 0.0;
00580     for (int k = 0;k < K;k++) sum += std::exp(tmpvecK[k]);
00581     return(std::log(sum));
00582   }
00583 }
00584 
00585 
00586 double MOG_generic::log_lhood(const vec &x_in)
00587 {
00588   if (do_checks) {
00589     it_assert(valid, "MOG_generic::log_lhood(): model not valid");
00590     it_assert(check_size(x_in), "MOG_generic::log_lhood(): x has wrong dimensionality");
00591   }
00592   return log_lhood_internal(x_in);
00593 }
00594 
00595 
00596 double MOG_generic::lhood_internal(const vec &x_in)
00597 {
00598 
00599   bool danger = paranoid;
00600 
00601   for (int k = 0;k < K;k++)  {
00602     double tmp = log_weights[k] + log_lhood_single_gaus_internal(x_in, k);
00603     tmpvecK[k] = tmp;
00604 
00605     if (tmp >= log_max_K)  danger = true;
00606   }
00607 
00608   if (danger) {
00609     double log_sum = tmpvecK[0];
00610     for (int k = 1; k < K; k++)  log_sum = log_add(log_sum, tmpvecK[k]);
00611     return(trunc_exp(log_sum));
00612   }
00613   else {
00614     double sum = 0.0;
00615     for (int k = 0;k < K;k++) sum += std::exp(tmpvecK[k]);
00616     return(sum);
00617   }
00618 }
00619 
00620 
00621 double MOG_generic::lhood(const vec &x_in)
00622 {
00623 
00624   if (do_checks) {
00625     it_assert(valid, "MOG_generic::lhood(): model not valid");
00626     it_assert(check_size(x_in), "MOG_generic::lhood(): x has wrong dimensionality");
00627   }
00628   return lhood_internal(x_in);
00629 }
00630 
00631 
00632 double MOG_generic::avg_log_lhood(const Array<vec> &X_in)
00633 {
00634   if (do_checks) {
00635     it_assert(valid, "MOG_generic::avg_log_lhood(): model not valid");
00636     it_assert(check_size(X_in), "MOG_generic::avg_log_lhood(): X is empty or at least one vector has the wrong dimensionality");
00637   }
00638 
00639   const int N = X_in.size();
00640   double acc = 0.0;
00641   for (int n = 0;n < N;n++)  acc += log_lhood_internal(X_in(n));
00642   return(acc / N);
00643 }
00644 
00645 
00646 
00647 } // namespace itpp
SourceForge Logo

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