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 //! update statistics to reflect new sample
00040 template<typename eT>
00041 template<typename T1>
00042 inline
00043 void
00044 running_stat_vec<eT>::operator() (const Base<typename get_pod_type<eT>::result, T1>& X)
00045   {
00046   arma_extra_debug_sigprint();
00047   
00048   //typedef typename get_pod_type<eT>::result T;
00049   
00050   const unwrap<T1>        tmp(X.get_ref());
00051   const Mat<eT>& sample = tmp.M;
00052   
00053   arma_check( (sample.is_finite() == false), "running_stat_vec: given sample has non-finite elements" );
00054   
00055   running_stat_vec_aux::update_stats(*this, sample);
00056   }
00057 
00058 
00059 
00060 //! update statistics to reflect new sample (version for complex numbers)
00061 template<typename eT>
00062 template<typename T1>
00063 inline
00064 void
00065 running_stat_vec<eT>::operator() (const Base<std::complex<typename get_pod_type<eT>::result>, T1>& X)
00066   {
00067   arma_extra_debug_sigprint();
00068   
00069   //typedef typename std::complex<typename get_pod_type<eT>::result> eT;
00070   
00071   const unwrap<T1>        tmp(X.get_ref());
00072   const Mat<eT>& sample = tmp.M;
00073   
00074   arma_check( (sample.is_finite() == false), "running_stat_vec: given sample has non-finite elements" );
00075   
00076   running_stat_vec_aux::update_stats(*this, sample);
00077   }
00078 
00079 
00080 
00081 //! set all statistics to zero
00082 template<typename eT>
00083 inline
00084 void
00085 running_stat_vec<eT>::reset()
00086   {
00087   arma_extra_debug_sigprint();
00088   
00089   counter.reset();
00090   
00091   r_mean.reset();
00092   r_var.reset();
00093   r_cov.reset();
00094   
00095   min_val.reset();
00096   max_val.reset();
00097   
00098   min_val_norm.reset();
00099   max_val_norm.reset();
00100   }
00101 
00102 
00103 
00104 //! mean or average value
00105 template<typename eT>
00106 inline
00107 Mat<eT>
00108 running_stat_vec<eT>::mean()
00109   const
00110   {
00111   arma_extra_debug_sigprint();
00112   
00113   return r_mean;
00114   }
00115 
00116 
00117 
00118 //! variance
00119 template<typename eT>
00120 inline
00121 Mat<typename get_pod_type<eT>::result>
00122 running_stat_vec<eT>::var(const u32 norm_type)
00123   const
00124   {
00125   arma_extra_debug_sigprint();
00126   
00127   const T N = counter.value();
00128   
00129   if(N > T(1))
00130     {
00131     if(norm_type == 0)
00132       {
00133       return r_var;
00134       }
00135     else
00136       {
00137       const T N_minus_1 = counter.value_minus_1();
00138       return (N_minus_1/N) * r_var;
00139       }
00140     }
00141   else
00142     {
00143     return zeros< Mat<typename get_pod_type<eT>::result> >(r_mean.n_rows, r_mean.n_cols);
00144     }
00145   
00146   }
00147 
00148 
00149 
00150 //! standard deviation
00151 template<typename eT>
00152 inline
00153 Mat<typename get_pod_type<eT>::result>
00154 running_stat_vec<eT>::stddev(const u32 norm_type)
00155   const
00156   {
00157   arma_extra_debug_sigprint();
00158   
00159   return sqrt( (*this).var(norm_type) );
00160   }
00161 
00162 
00163 
00164 //! covariance
00165 template<typename eT>
00166 inline
00167 Mat<eT>
00168 running_stat_vec<eT>::cov(const u32 norm_type)
00169   const
00170   {
00171   arma_extra_debug_sigprint();
00172   
00173   if(calc_cov == true)
00174     {
00175     const T N = counter.value();
00176     
00177     if(N > T(1))
00178       {
00179       if(norm_type == 0)
00180         {
00181         return r_cov;
00182         }
00183       else
00184         {
00185         const T N_minus_1 = counter.value_minus_1();
00186         return (N_minus_1/N) * r_cov;
00187         }
00188       }
00189     else
00190       {
00191       return zeros< Mat<eT> >(r_mean.n_rows, r_mean.n_cols);
00192       }
00193     }
00194   else
00195     {
00196     return Mat<eT>();
00197     }
00198   
00199   }
00200 
00201 
00202 
00203 //! vector with minimum values
00204 template<typename eT>
00205 inline
00206 Mat<eT>
00207 running_stat_vec<eT>::min()
00208 const
00209   {
00210   arma_extra_debug_sigprint();
00211 
00212   return min_val;
00213   }
00214 
00215 
00216 
00217 //! vector with maximum values
00218 template<typename eT>
00219 inline
00220 Mat<eT>
00221 running_stat_vec<eT>::max()
00222 const
00223   {
00224   arma_extra_debug_sigprint();
00225 
00226   return max_val;
00227   }
00228 
00229 
00230 
00231 //
00232 
00233 
00234 
00235 //! update statistics to reflect new sample
00236 template<typename eT>
00237 inline
00238 void
00239 running_stat_vec_aux::update_stats(running_stat_vec<eT>& x, const Mat<eT>& sample)
00240   {
00241   arma_extra_debug_sigprint();
00242   
00243   typedef typename running_stat_vec<eT>::T T;
00244   
00245   const T N = x.counter.value();
00246   
00247   if(N > T(0))
00248     {
00249     arma_debug_assert_same_size(x.r_mean, sample, "running_stat_vec(): dimensionality mismatch");
00250     
00251     const u32 n_elem      = sample.n_elem;
00252     const eT* sample_mem  = sample.memptr();
00253           eT* r_mean_mem  = x.r_mean.memptr();
00254            T* r_var_mem   = x.r_var.memptr();
00255           eT* min_val_mem = x.min_val.memptr();
00256           eT* max_val_mem = x.max_val.memptr();
00257     
00258     const T  N_plus_1   = x.counter.value_plus_1();
00259     const T  N_minus_1  = x.counter.value_minus_1();
00260     
00261     if(x.calc_cov == true)
00262       {
00263       const Mat<eT> tmp1 = sample - x.r_mean;
00264       
00265       Mat<eT> tmp2;
00266       
00267       if(sample.n_cols == 1)
00268         {
00269         tmp2 = tmp1*trans(tmp1);
00270         }
00271       else
00272         {
00273         tmp2 = trans(tmp1)*tmp1;
00274         }
00275       
00276       x.r_cov *= (N_minus_1/N);
00277       x.r_cov += tmp2 / N_plus_1;
00278       }
00279     
00280     
00281     for(u32 i=0; i<n_elem; ++i)
00282       {
00283       const eT val = sample_mem[i];
00284       
00285       if(val < min_val_mem[i])
00286         {
00287         min_val_mem[i] = val;
00288         }
00289       
00290       if(val > max_val_mem[i])
00291         {
00292         max_val_mem[i] = val;
00293         }
00294         
00295       const eT r_mean_val = r_mean_mem[i];
00296       const eT tmp        = val - r_mean_val;
00297     
00298       r_var_mem[i] = N_minus_1/N * r_var_mem[i] + (tmp*tmp)/N_plus_1;
00299       
00300       r_mean_mem[i] = r_mean_val + (val - r_mean_val)/N_plus_1;
00301       }
00302     }
00303   else
00304     {
00305     arma_debug_check( (sample.is_vec() == false), "running_stat_vec(): given sample is not a vector");
00306     
00307     x.r_mean.set_size(sample.n_rows, sample.n_cols);
00308     
00309     x.r_var.zeros(sample.n_rows, sample.n_cols);
00310     
00311     if(x.calc_cov == true)
00312       {
00313       x.r_cov.zeros(sample.n_elem, sample.n_elem);
00314       }
00315     
00316     x.min_val.set_size(sample.n_rows, sample.n_cols);
00317     x.max_val.set_size(sample.n_rows, sample.n_cols);
00318     
00319     
00320     const u32 n_elem      = sample.n_elem;
00321     const eT* sample_mem  = sample.memptr();
00322           eT* r_mean_mem  = x.r_mean.memptr();
00323           eT* min_val_mem = x.min_val.memptr();
00324           eT* max_val_mem = x.max_val.memptr();
00325           
00326     
00327     for(u32 i=0; i<n_elem; ++i)
00328       {
00329       const eT val = sample_mem[i];
00330       
00331       r_mean_mem[i]  = val;
00332       min_val_mem[i] = val;
00333       max_val_mem[i] = val;
00334       }
00335     }
00336   
00337   x.counter++;
00338   }
00339 
00340 
00341 
00342 //! update statistics to reflect new sample (version for complex numbers)
00343 template<typename T>
00344 inline
00345 void
00346 running_stat_vec_aux::update_stats(running_stat_vec< std::complex<T> >& x, const Mat<T>& sample)
00347   {
00348   arma_extra_debug_sigprint();
00349   
00350   const Mat< std::complex<T> > tmp = conv_to< Mat< std::complex<T> > >::from(sample);
00351   
00352   running_stat_vec_aux::update_stats(x, tmp);
00353   }
00354 
00355 
00356 
00357 //! alter statistics to reflect new sample (version for complex numbers)
00358 template<typename T>
00359 inline
00360 void
00361 running_stat_vec_aux::update_stats(running_stat_vec< std::complex<T> >& x, const Mat< std::complex<T> >& sample)
00362   {
00363   arma_extra_debug_sigprint();
00364   
00365   typedef typename std::complex<T> eT;
00366   
00367   const T N = x.counter.value();
00368   
00369   if(N > T(0))
00370     {
00371     arma_debug_assert_same_size(x.r_mean, sample, "running_stat_vec(): dimensionality mismatch");
00372     
00373     const u32 n_elem           = sample.n_elem;
00374     const eT* sample_mem       = sample.memptr();
00375           eT* r_mean_mem       = x.r_mean.memptr();
00376            T* r_var_mem        = x.r_var.memptr();
00377           eT* min_val_mem      = x.min_val.memptr();
00378           eT* max_val_mem      = x.max_val.memptr();
00379            T* min_val_norm_mem = x.min_val_norm.memptr();
00380            T* max_val_norm_mem = x.max_val_norm.memptr();
00381     
00382     const T  N_plus_1   = x.counter.value_plus_1();
00383     const T  N_minus_1  = x.counter.value_minus_1();
00384     
00385     if(x.calc_cov == true)
00386       {
00387       const Mat<eT> tmp1 = sample - x.r_mean;
00388       
00389       Mat<eT> tmp2;
00390       
00391       if(sample.n_cols == 1)
00392         {
00393         tmp2 = conj(tmp1)*trans(tmp1);
00394         }
00395       else
00396         {
00397         tmp2 = trans(conj(tmp1))*tmp1;
00398         }
00399       
00400       x.r_cov *= (N_minus_1/N);
00401       x.r_cov += tmp2 / N_plus_1;
00402       }
00403     
00404     
00405     for(u32 i=0; i<n_elem; ++i)
00406       {
00407       const eT& val      = sample_mem[i];
00408       const  T  val_norm = std::norm(val);
00409       
00410       if(val_norm < min_val_norm_mem[i])
00411         {
00412         min_val_norm_mem[i] = val_norm;
00413         min_val_mem[i]      = val;
00414         }
00415       
00416       if(val_norm > max_val_norm_mem[i])
00417         {
00418         max_val_norm_mem[i] = val_norm;
00419         max_val_mem[i]      = val;
00420         }
00421       
00422       const eT& r_mean_val = r_mean_mem[i];
00423       
00424       r_var_mem[i] = N_minus_1/N * r_var_mem[i] + std::norm(val - r_mean_val)/N_plus_1;
00425       
00426       r_mean_mem[i] = r_mean_val + (val - r_mean_val)/N_plus_1;
00427       }
00428     
00429     }
00430   else
00431     {
00432     arma_debug_check( (sample.is_vec() == false), "running_stat_vec(): given sample is not a vector");
00433     
00434     x.r_mean.set_size(sample.n_rows, sample.n_cols);
00435     
00436     x.r_var.zeros(sample.n_rows, sample.n_cols);
00437     
00438     if(x.calc_cov == true)
00439       {
00440       x.r_cov.zeros(sample.n_elem, sample.n_elem);
00441       }
00442     
00443     x.min_val.set_size(sample.n_rows, sample.n_cols);
00444     x.max_val.set_size(sample.n_rows, sample.n_cols);
00445     
00446     x.min_val_norm.set_size(sample.n_rows, sample.n_cols);
00447     x.max_val_norm.set_size(sample.n_rows, sample.n_cols);
00448     
00449     
00450     const u32 n_elem           = sample.n_elem;
00451     const eT* sample_mem       = sample.memptr();
00452           eT* r_mean_mem       = x.r_mean.memptr();
00453           eT* min_val_mem      = x.min_val.memptr();
00454           eT* max_val_mem      = x.max_val.memptr();
00455            T* min_val_norm_mem = x.min_val_norm.memptr();
00456            T* max_val_norm_mem = x.max_val_norm.memptr();
00457     
00458     for(u32 i=0; i<n_elem; ++i)
00459       {
00460       const eT& val      = sample_mem[i];
00461       const  T  val_norm = std::norm(val);
00462       
00463       r_mean_mem[i]  = val;
00464       min_val_mem[i] = val;
00465       max_val_mem[i] = val;
00466       
00467       min_val_norm_mem[i] = val_norm;
00468       max_val_norm_mem[i] = val_norm;
00469       }
00470     }
00471   
00472   x.counter++;
00473   }
00474 
00475 
00476 
00477 //! @}