$search
00001 // Copyright (C) 2009-2011 NICTA (www.nicta.com.au) 00002 // Copyright (C) 2009-2011 Conrad Sanderson 00003 // 00004 // This file is part of the Armadillo C++ library. 00005 // It is provided without any warranty of fitness 00006 // for any purpose. You can redistribute this file 00007 // and/or modify it under the terms of the GNU 00008 // Lesser General Public License (LGPL) as published 00009 // by the Free Software Foundation, either version 3 00010 // of the License or (at your option) any later version. 00011 // (see http://www.opensource.org/licenses for more info) 00012 00013 00016 00017 00018 00019 template<typename eT> 00020 running_stat_vec<eT>::~running_stat_vec() 00021 { 00022 arma_extra_debug_sigprint_this(this); 00023 } 00024 00025 00026 00027 template<typename eT> 00028 running_stat_vec<eT>::running_stat_vec(const bool in_calc_cov) 00029 : calc_cov(in_calc_cov) 00030 { 00031 arma_extra_debug_sigprint_this(this); 00032 } 00033 00034 00035 00036 template<typename eT> 00037 running_stat_vec<eT>::running_stat_vec(const running_stat_vec<eT>& in_rsv) 00038 : calc_cov (in_rsv.calc_cov) 00039 , counter (in_rsv.counter) 00040 , r_mean (in_rsv.r_mean) 00041 , r_var (in_rsv.r_var) 00042 , r_cov (in_rsv.r_cov) 00043 , min_val (in_rsv.min_val) 00044 , max_val (in_rsv.max_val) 00045 , min_val_norm(in_rsv.min_val_norm) 00046 , max_val_norm(in_rsv.max_val_norm) 00047 { 00048 arma_extra_debug_sigprint_this(this); 00049 } 00050 00051 00052 00053 template<typename eT> 00054 const running_stat_vec<eT>& 00055 running_stat_vec<eT>::operator=(const running_stat_vec<eT>& in_rsv) 00056 { 00057 arma_extra_debug_sigprint(); 00058 00059 access::rw(calc_cov) = in_rsv.calc_cov; 00060 00061 counter = in_rsv.counter; 00062 r_mean = in_rsv.r_mean; 00063 r_var = in_rsv.r_var; 00064 r_cov = in_rsv.r_cov; 00065 min_val = in_rsv.min_val; 00066 max_val = in_rsv.max_val; 00067 min_val_norm = in_rsv.min_val_norm; 00068 max_val_norm = in_rsv.max_val_norm; 00069 00070 return *this; 00071 } 00072 00073 00074 00076 template<typename eT> 00077 template<typename T1> 00078 arma_hot 00079 inline 00080 void 00081 running_stat_vec<eT>::operator() (const Base<typename get_pod_type<eT>::result, T1>& X) 00082 { 00083 arma_extra_debug_sigprint(); 00084 00085 //typedef typename get_pod_type<eT>::result T; 00086 00087 const unwrap<T1> tmp(X.get_ref()); 00088 const Mat<eT>& sample = tmp.M; 00089 00090 if( sample.is_empty() ) 00091 { 00092 return; 00093 } 00094 00095 if( sample.is_finite() == false ) 00096 { 00097 arma_warn(true, "running_stat_vec: sample ignored as it has non-finite elements"); 00098 return; 00099 } 00100 00101 running_stat_vec_aux::update_stats(*this, sample); 00102 } 00103 00104 00105 00107 template<typename eT> 00108 template<typename T1> 00109 arma_hot 00110 inline 00111 void 00112 running_stat_vec<eT>::operator() (const Base<std::complex<typename get_pod_type<eT>::result>, T1>& X) 00113 { 00114 arma_extra_debug_sigprint(); 00115 00116 //typedef typename std::complex<typename get_pod_type<eT>::result> eT; 00117 00118 const unwrap<T1> tmp(X.get_ref()); 00119 const Mat<eT>& sample = tmp.M; 00120 00121 if( sample.is_empty() ) 00122 { 00123 return; 00124 } 00125 00126 if( sample.is_finite() == false ) 00127 { 00128 arma_warn(true, "running_stat_vec: sample ignored as it has non-finite elements"); 00129 return; 00130 } 00131 00132 running_stat_vec_aux::update_stats(*this, sample); 00133 } 00134 00135 00136 00138 template<typename eT> 00139 inline 00140 void 00141 running_stat_vec<eT>::reset() 00142 { 00143 arma_extra_debug_sigprint(); 00144 00145 counter.reset(); 00146 00147 r_mean.reset(); 00148 r_var.reset(); 00149 r_cov.reset(); 00150 00151 min_val.reset(); 00152 max_val.reset(); 00153 00154 min_val_norm.reset(); 00155 max_val_norm.reset(); 00156 00157 r_var_dummy.reset(); 00158 r_cov_dummy.reset(); 00159 00160 tmp1.reset(); 00161 tmp2.reset(); 00162 } 00163 00164 00165 00167 template<typename eT> 00168 inline 00169 const Mat<eT>& 00170 running_stat_vec<eT>::mean() const 00171 { 00172 arma_extra_debug_sigprint(); 00173 00174 return r_mean; 00175 } 00176 00177 00178 00180 template<typename eT> 00181 inline 00182 const Mat<typename get_pod_type<eT>::result>& 00183 running_stat_vec<eT>::var(const uword norm_type) 00184 { 00185 arma_extra_debug_sigprint(); 00186 00187 const T N = counter.value(); 00188 00189 if(N > T(1)) 00190 { 00191 if(norm_type == 0) 00192 { 00193 return r_var; 00194 } 00195 else 00196 { 00197 const T N_minus_1 = counter.value_minus_1(); 00198 00199 r_var_dummy = (N_minus_1/N) * r_var; 00200 00201 return r_var_dummy; 00202 } 00203 } 00204 else 00205 { 00206 r_var_dummy.zeros(r_mean.n_rows, r_mean.n_cols); 00207 00208 return r_var_dummy; 00209 } 00210 00211 } 00212 00213 00214 00216 template<typename eT> 00217 inline 00218 Mat<typename get_pod_type<eT>::result> 00219 running_stat_vec<eT>::stddev(const uword norm_type) const 00220 { 00221 arma_extra_debug_sigprint(); 00222 00223 const T N = counter.value(); 00224 00225 if(N > T(1)) 00226 { 00227 if(norm_type == 0) 00228 { 00229 return sqrt(r_var); 00230 } 00231 else 00232 { 00233 const T N_minus_1 = counter.value_minus_1(); 00234 00235 return sqrt( (N_minus_1/N) * r_var ); 00236 } 00237 } 00238 else 00239 { 00240 return Mat<T>(); 00241 } 00242 } 00243 00244 00245 00247 template<typename eT> 00248 inline 00249 const Mat<eT>& 00250 running_stat_vec<eT>::cov(const uword norm_type) 00251 { 00252 arma_extra_debug_sigprint(); 00253 00254 if(calc_cov == true) 00255 { 00256 const T N = counter.value(); 00257 00258 if(N > T(1)) 00259 { 00260 if(norm_type == 0) 00261 { 00262 return r_cov; 00263 } 00264 else 00265 { 00266 const T N_minus_1 = counter.value_minus_1(); 00267 00268 r_cov_dummy = (N_minus_1/N) * r_cov; 00269 00270 return r_cov_dummy; 00271 } 00272 } 00273 else 00274 { 00275 r_cov_dummy.zeros(r_mean.n_rows, r_mean.n_cols); 00276 00277 return r_cov_dummy; 00278 } 00279 } 00280 else 00281 { 00282 r_cov_dummy.reset(); 00283 00284 return r_cov_dummy; 00285 } 00286 00287 } 00288 00289 00290 00292 template<typename eT> 00293 inline 00294 const Mat<eT>& 00295 running_stat_vec<eT>::min() const 00296 { 00297 arma_extra_debug_sigprint(); 00298 00299 return min_val; 00300 } 00301 00302 00303 00305 template<typename eT> 00306 inline 00307 const Mat<eT>& 00308 running_stat_vec<eT>::max() const 00309 { 00310 arma_extra_debug_sigprint(); 00311 00312 return max_val; 00313 } 00314 00315 00316 00318 template<typename eT> 00319 inline 00320 typename get_pod_type<eT>::result 00321 running_stat_vec<eT>::count() const 00322 { 00323 arma_extra_debug_sigprint(); 00324 00325 return counter.value(); 00326 } 00327 00328 00329 00330 // 00331 00332 00333 00335 template<typename eT> 00336 inline 00337 void 00338 running_stat_vec_aux::update_stats(running_stat_vec<eT>& x, const Mat<eT>& sample) 00339 { 00340 arma_extra_debug_sigprint(); 00341 00342 typedef typename running_stat_vec<eT>::T T; 00343 00344 const T N = x.counter.value(); 00345 00346 if(N > T(0)) 00347 { 00348 arma_debug_assert_same_size(x.r_mean, sample, "running_stat_vec(): dimensionality mismatch"); 00349 00350 const uword n_elem = sample.n_elem; 00351 const eT* sample_mem = sample.memptr(); 00352 eT* r_mean_mem = x.r_mean.memptr(); 00353 T* r_var_mem = x.r_var.memptr(); 00354 eT* min_val_mem = x.min_val.memptr(); 00355 eT* max_val_mem = x.max_val.memptr(); 00356 00357 const T N_plus_1 = x.counter.value_plus_1(); 00358 const T N_minus_1 = x.counter.value_minus_1(); 00359 00360 if(x.calc_cov == true) 00361 { 00362 Mat<eT>& tmp1 = x.tmp1; 00363 Mat<eT>& tmp2 = x.tmp2; 00364 00365 tmp1 = sample - x.r_mean; 00366 00367 if(sample.n_cols == 1) 00368 { 00369 tmp2 = tmp1*trans(tmp1); 00370 } 00371 else 00372 { 00373 tmp2 = trans(tmp1)*tmp1; 00374 } 00375 00376 x.r_cov *= (N_minus_1/N); 00377 x.r_cov += tmp2 / N_plus_1; 00378 } 00379 00380 00381 for(uword i=0; i<n_elem; ++i) 00382 { 00383 const eT val = sample_mem[i]; 00384 00385 if(val < min_val_mem[i]) 00386 { 00387 min_val_mem[i] = val; 00388 } 00389 00390 if(val > max_val_mem[i]) 00391 { 00392 max_val_mem[i] = val; 00393 } 00394 00395 const eT r_mean_val = r_mean_mem[i]; 00396 const eT tmp = val - r_mean_val; 00397 00398 r_var_mem[i] = N_minus_1/N * r_var_mem[i] + (tmp*tmp)/N_plus_1; 00399 00400 r_mean_mem[i] = r_mean_val + (val - r_mean_val)/N_plus_1; 00401 } 00402 } 00403 else 00404 { 00405 arma_debug_check( (sample.is_vec() == false), "running_stat_vec(): given sample is not a vector"); 00406 00407 x.r_mean.set_size(sample.n_rows, sample.n_cols); 00408 00409 x.r_var.zeros(sample.n_rows, sample.n_cols); 00410 00411 if(x.calc_cov == true) 00412 { 00413 x.r_cov.zeros(sample.n_elem, sample.n_elem); 00414 } 00415 00416 x.min_val.set_size(sample.n_rows, sample.n_cols); 00417 x.max_val.set_size(sample.n_rows, sample.n_cols); 00418 00419 00420 const uword n_elem = sample.n_elem; 00421 const eT* sample_mem = sample.memptr(); 00422 eT* r_mean_mem = x.r_mean.memptr(); 00423 eT* min_val_mem = x.min_val.memptr(); 00424 eT* max_val_mem = x.max_val.memptr(); 00425 00426 00427 for(uword i=0; i<n_elem; ++i) 00428 { 00429 const eT val = sample_mem[i]; 00430 00431 r_mean_mem[i] = val; 00432 min_val_mem[i] = val; 00433 max_val_mem[i] = val; 00434 } 00435 } 00436 00437 x.counter++; 00438 } 00439 00440 00441 00443 template<typename T> 00444 inline 00445 void 00446 running_stat_vec_aux::update_stats(running_stat_vec< std::complex<T> >& x, const Mat<T>& sample) 00447 { 00448 arma_extra_debug_sigprint(); 00449 00450 const Mat< std::complex<T> > tmp = conv_to< Mat< std::complex<T> > >::from(sample); 00451 00452 running_stat_vec_aux::update_stats(x, tmp); 00453 } 00454 00455 00456 00458 template<typename T> 00459 inline 00460 void 00461 running_stat_vec_aux::update_stats(running_stat_vec< std::complex<T> >& x, const Mat< std::complex<T> >& sample) 00462 { 00463 arma_extra_debug_sigprint(); 00464 00465 typedef typename std::complex<T> eT; 00466 00467 const T N = x.counter.value(); 00468 00469 if(N > T(0)) 00470 { 00471 arma_debug_assert_same_size(x.r_mean, sample, "running_stat_vec(): dimensionality mismatch"); 00472 00473 const uword n_elem = sample.n_elem; 00474 const eT* sample_mem = sample.memptr(); 00475 eT* r_mean_mem = x.r_mean.memptr(); 00476 T* r_var_mem = x.r_var.memptr(); 00477 eT* min_val_mem = x.min_val.memptr(); 00478 eT* max_val_mem = x.max_val.memptr(); 00479 T* min_val_norm_mem = x.min_val_norm.memptr(); 00480 T* max_val_norm_mem = x.max_val_norm.memptr(); 00481 00482 const T N_plus_1 = x.counter.value_plus_1(); 00483 const T N_minus_1 = x.counter.value_minus_1(); 00484 00485 if(x.calc_cov == true) 00486 { 00487 Mat<eT>& tmp1 = x.tmp1; 00488 Mat<eT>& tmp2 = x.tmp2; 00489 00490 tmp1 = sample - x.r_mean; 00491 00492 if(sample.n_cols == 1) 00493 { 00494 tmp2 = arma::conj(tmp1)*strans(tmp1); 00495 } 00496 else 00497 { 00498 tmp2 = trans(tmp1)*tmp1; //tmp2 = strans(conj(tmp1))*tmp1; 00499 } 00500 00501 x.r_cov *= (N_minus_1/N); 00502 x.r_cov += tmp2 / N_plus_1; 00503 } 00504 00505 00506 for(uword i=0; i<n_elem; ++i) 00507 { 00508 const eT& val = sample_mem[i]; 00509 const T val_norm = std::norm(val); 00510 00511 if(val_norm < min_val_norm_mem[i]) 00512 { 00513 min_val_norm_mem[i] = val_norm; 00514 min_val_mem[i] = val; 00515 } 00516 00517 if(val_norm > max_val_norm_mem[i]) 00518 { 00519 max_val_norm_mem[i] = val_norm; 00520 max_val_mem[i] = val; 00521 } 00522 00523 const eT& r_mean_val = r_mean_mem[i]; 00524 00525 r_var_mem[i] = N_minus_1/N * r_var_mem[i] + std::norm(val - r_mean_val)/N_plus_1; 00526 00527 r_mean_mem[i] = r_mean_val + (val - r_mean_val)/N_plus_1; 00528 } 00529 00530 } 00531 else 00532 { 00533 arma_debug_check( (sample.is_vec() == false), "running_stat_vec(): given sample is not a vector"); 00534 00535 x.r_mean.set_size(sample.n_rows, sample.n_cols); 00536 00537 x.r_var.zeros(sample.n_rows, sample.n_cols); 00538 00539 if(x.calc_cov == true) 00540 { 00541 x.r_cov.zeros(sample.n_elem, sample.n_elem); 00542 } 00543 00544 x.min_val.set_size(sample.n_rows, sample.n_cols); 00545 x.max_val.set_size(sample.n_rows, sample.n_cols); 00546 00547 x.min_val_norm.set_size(sample.n_rows, sample.n_cols); 00548 x.max_val_norm.set_size(sample.n_rows, sample.n_cols); 00549 00550 00551 const uword n_elem = sample.n_elem; 00552 const eT* sample_mem = sample.memptr(); 00553 eT* r_mean_mem = x.r_mean.memptr(); 00554 eT* min_val_mem = x.min_val.memptr(); 00555 eT* max_val_mem = x.max_val.memptr(); 00556 T* min_val_norm_mem = x.min_val_norm.memptr(); 00557 T* max_val_norm_mem = x.max_val_norm.memptr(); 00558 00559 for(uword i=0; i<n_elem; ++i) 00560 { 00561 const eT& val = sample_mem[i]; 00562 const T val_norm = std::norm(val); 00563 00564 r_mean_mem[i] = val; 00565 min_val_mem[i] = val; 00566 max_val_mem[i] = val; 00567 00568 min_val_norm_mem[i] = val_norm; 00569 max_val_norm_mem[i] = val_norm; 00570 } 00571 } 00572 00573 x.counter++; 00574 } 00575 00576 00577