00001 00030 #include <itpp/base/math/log_exp.h> 00031 #include <itpp/stat/mog_diag.h> 00032 #include <cstdlib> 00033 00034 00035 namespace itpp { 00036 00037 double MOG_diag::log_lhood_single_gaus_internal(const double * c_x_in, const int k) const { 00038 00039 const double * c_mean = c_means[k]; 00040 const double * c_diag_cov_inv_etc = c_diag_covs_inv_etc[k]; 00041 00042 double acc = 0.0; 00043 00044 for(int d=0; d<D; d++) { 00045 double tmp_val = c_x_in[d] - c_mean[d]; 00046 acc += (tmp_val*tmp_val) * c_diag_cov_inv_etc[d]; 00047 } 00048 return(c_log_det_etc[k] - acc); 00049 } 00050 00051 00052 double MOG_diag::log_lhood_single_gaus_internal(const vec &x_in, const int k) const { 00053 return log_lhood_single_gaus_internal(x_in._data(), k); 00054 } 00055 00056 00057 double MOG_diag::log_lhood_single_gaus(const double * c_x_in, const int k) const { 00058 if(do_checks) { 00059 it_assert(valid, "MOG_diag::log_lhood_single_gaus(): model not valid"); 00060 it_assert( ( (k>=0) && (k<K) ), "MOG::log_lhood_single_gaus(): k specifies a non-existant Gaussian"); 00061 } 00062 return log_lhood_single_gaus_internal(c_x_in,k); 00063 } 00064 00065 00066 double MOG_diag::log_lhood_single_gaus(const vec &x_in, const int k) const { 00067 if(do_checks) { 00068 it_assert(valid, "MOG_diag::log_lhood_single_gaus(): model not valid"); 00069 it_assert(check_size(x_in), "MOG_diag::log_lhood_single_gaus(): x has wrong dimensionality"); 00070 it_assert( ( (k>=0) && (k<K) ), "MOG::log_lhood_single_gaus(): k specifies a non-existant Gaussian"); 00071 } 00072 return log_lhood_single_gaus_internal(x_in._data(),k); 00073 } 00074 00075 00076 double MOG_diag::log_lhood_internal(const double * c_x_in) { 00077 00078 bool danger = paranoid; 00079 00080 for(int k=0;k<K;k++) { 00081 double tmp = c_log_weights[k] + log_lhood_single_gaus_internal(c_x_in,k); 00082 c_tmpvecK[k] = tmp; 00083 00084 if(tmp >= log_max_K) danger = true; 00085 } 00086 00087 00088 if(danger) { 00089 double log_sum = c_tmpvecK[0]; for(int k=1; k<K; k++) log_sum = log_add( log_sum, c_tmpvecK[k] ); 00090 return(log_sum); 00091 } 00092 else { 00093 double sum = 0.0; for(int k=0;k<K;k++) sum += std::exp(c_tmpvecK[k]); 00094 return(std::log(sum)); 00095 } 00096 } 00097 00098 00099 double MOG_diag::log_lhood_internal(const vec &x_in) { 00100 return log_lhood_internal(x_in._data()); 00101 } 00102 00103 00104 double MOG_diag::log_lhood(const vec &x_in) { 00105 if(do_checks) { 00106 it_assert(valid, "MOG_diag::log_lhood(): model not valid"); 00107 it_assert(check_size(x_in), "MOG_diag::log_lhood(): x has wrong dimensionality"); 00108 } 00109 return log_lhood_internal(x_in._data()); 00110 } 00111 00112 00113 double MOG_diag::log_lhood(const double * c_x_in) { 00114 if(do_checks) { 00115 it_assert(valid, "MOG_diag::log_lhood(): model not valid"); 00116 it_assert( (c_x_in != 0), "MOG_diag::log_lhood(): c_x_in is a null pointer"); 00117 } 00118 00119 return log_lhood_internal(c_x_in); 00120 } 00121 00122 00123 double MOG_diag::lhood_internal(const double * c_x_in) { 00124 00125 bool danger = paranoid; 00126 00127 for(int k=0;k<K;k++) { 00128 double tmp = c_log_weights[k] + log_lhood_single_gaus_internal(c_x_in,k); 00129 c_tmpvecK[k] = tmp; 00130 00131 if(tmp >= log_max_K) danger = true; 00132 } 00133 00134 00135 if(danger) { 00136 double log_sum = c_tmpvecK[0]; for(int k=1; k<K; k++) log_sum = log_add( log_sum, c_tmpvecK[k] ); 00137 return(trunc_exp(log_sum)); 00138 } 00139 else { 00140 double sum = 0.0; for(int k=0;k<K;k++) sum += std::exp(c_tmpvecK[k]); 00141 return(sum); 00142 } 00143 } 00144 00145 double MOG_diag::lhood_internal(const vec &x_in) { return lhood_internal(x_in._data()); } 00146 00147 double MOG_diag::lhood(const vec &x_in) { 00148 if(do_checks) { 00149 it_assert(valid, "MOG_diag::lhood(): model not valid"); 00150 it_assert(check_size(x_in), "MOG_diag::lhood(): x has wrong dimensionality"); 00151 } 00152 return lhood_internal(x_in._data()); 00153 } 00154 00155 00156 double MOG_diag::lhood(const double * c_x_in) { 00157 if(do_checks) { 00158 it_assert(valid, "MOG_diag::lhood(): model not valid"); 00159 it_assert( (c_x_in != 0), "MOG_diag::lhood(): c_x_in is a null pointer"); 00160 } 00161 00162 return lhood_internal(c_x_in); 00163 } 00164 00165 00166 double MOG_diag::avg_log_lhood(const double ** c_x_in, const int N) { 00167 if(do_checks) { 00168 it_assert(valid, "MOG_diag::avg_log_lhood(): model not valid"); 00169 it_assert( (c_x_in != 0), "MOG_diag::avg_log_lhood(): c_x_in is a null pointer"); 00170 it_assert( (N >= 0), "MOG_diag::avg_log_lhood(): N is zero or negative"); 00171 } 00172 00173 double acc = 0.0; for(int n=0;n<N;n++) acc += log_lhood_internal(c_x_in[n]); 00174 return(acc/N); 00175 } 00176 00177 00178 double MOG_diag::avg_log_lhood(const Array<vec> &X_in) { 00179 if(do_checks) { 00180 it_assert(valid, "MOG_diag::avg_log_lhood(): model not valid"); 00181 it_assert(check_size(X_in), "MOG_diag::avg_log_lhood(): X is empty or at least one vector has the wrong dimensionality"); 00182 } 00183 const int N = X_in.size(); 00184 double acc = 0.0; 00185 for(int n=0;n<N;n++) acc += log_lhood_internal(X_in(n)._data()); 00186 return(acc/N); 00187 } 00188 00189 void MOG_diag::zero_all_ptrs() { 00190 c_means = 0; 00191 c_diag_covs = 0; 00192 c_diag_covs_inv_etc = 0; 00193 c_weights = 0; 00194 c_log_weights = 0; 00195 c_log_det_etc = 0; 00196 c_tmpvecK = 0; 00197 } 00198 00199 00200 void MOG_diag::free_all_ptrs() { 00201 c_means = disable_c_access(c_means); 00202 c_diag_covs = disable_c_access(c_diag_covs); 00203 c_diag_covs_inv_etc = disable_c_access(c_diag_covs_inv_etc); 00204 c_weights = disable_c_access(c_weights); 00205 c_log_weights = disable_c_access(c_log_weights); 00206 c_log_det_etc = disable_c_access(c_log_det_etc); 00207 c_tmpvecK = disable_c_access(c_tmpvecK); 00208 } 00209 00210 00211 void MOG_diag::setup_means() { 00212 MOG_generic::setup_means(); 00213 disable_c_access(c_means); 00214 c_means = enable_c_access(means); 00215 } 00216 00217 00218 void MOG_diag::setup_covs() { 00219 MOG_generic::setup_covs(); 00220 if(full) return; 00221 00222 disable_c_access(c_diag_covs); 00223 disable_c_access(c_diag_covs_inv_etc); 00224 disable_c_access(c_log_det_etc); 00225 00226 c_diag_covs = enable_c_access(diag_covs); 00227 c_diag_covs_inv_etc = enable_c_access(diag_covs_inv_etc); 00228 c_log_det_etc = enable_c_access(log_det_etc); 00229 } 00230 00231 00232 void MOG_diag::setup_weights() { 00233 MOG_generic::setup_weights(); 00234 00235 disable_c_access(c_weights); 00236 disable_c_access(c_log_weights); 00237 00238 c_weights = enable_c_access(weights); 00239 c_log_weights = enable_c_access(log_weights); 00240 } 00241 00242 00243 void MOG_diag::setup_misc() { 00244 disable_c_access(c_tmpvecK); 00245 tmpvecK.set_size(K); 00246 c_tmpvecK = enable_c_access(tmpvecK); 00247 00248 MOG_generic::setup_misc(); 00249 if(full) convert_to_diag_internal(); 00250 } 00251 00252 00253 void MOG_diag::load(const std::string &name_in) { 00254 MOG_generic::load(name_in); 00255 if(full) convert_to_diag(); 00256 } 00257 00258 00259 double ** MOG_diag::enable_c_access(Array<vec> & A_in) { 00260 int rows = A_in.size(); 00261 double ** A = (double **)std::malloc(rows*sizeof(double *)); 00262 if(A) for(int row=0;row<rows;row++) A[row] = A_in(row)._data(); 00263 return(A); 00264 } 00265 00266 int ** MOG_diag::enable_c_access(Array<ivec> & A_in) { 00267 int rows = A_in.size(); 00268 int ** A = (int **)std::malloc(rows*sizeof(int *)); 00269 if(A) for(int row=0;row<rows;row++) A[row] = A_in(row)._data(); 00270 return(A); 00271 } 00272 00273 double ** MOG_diag::disable_c_access(double ** A_in) { if(A_in) std::free(A_in); return(0); } 00274 int ** MOG_diag::disable_c_access(int ** A_in) { if(A_in) std::free(A_in); return(0); } 00275 00276 double * MOG_diag::enable_c_access(vec & v_in) { return v_in._data(); } 00277 int * MOG_diag::enable_c_access(ivec & v_in) { return v_in._data(); } 00278 00279 double * MOG_diag::disable_c_access(double * v_in) { return(0); } 00280 int * MOG_diag::disable_c_access(int * v_in) { return(0); } 00281 00282 }
Generated on Sun Apr 20 12:40:07 2008 for IT++ by Doxygen 1.5.5