running_stat_vec_meat.hpp

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