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
Generated on Sun Dec 9 17:26:19 2007 for IT++ by Doxygen 1.5.4