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

Generated on Sun Dec 9 17:38:50 2007 for IT++ by Doxygen 1.5.4