op_var_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 op_var
00018 //! @{
00019 
00020 
00021 //! find the variance of an array
00022 template<typename eT>
00023 inline
00024 eT
00025 op_var::direct_var(const eT* const X, const u32 n_elem, const u32 norm_type)
00026   {
00027   arma_extra_debug_sigprint();
00028   
00029   eT acc1 = eT(0);
00030   
00031   for(u32 i=0; i<n_elem; ++i)
00032     {
00033     acc1 += X[i];
00034     }
00035   
00036   const eT div_val = (n_elem > 0) ? eT(n_elem) : eT(1);
00037   acc1 /= div_val;
00038   
00039   eT acc2 = eT(0);
00040   eT acc3 = eT(0);
00041 
00042   for(u32 i=0; i<n_elem; ++i)
00043     {
00044     const eT tmp = acc1 - X[i];
00045   
00046     acc2 += tmp*tmp;
00047     acc3 += tmp;
00048     }
00049   
00050   
00051   const eT norm_val = (norm_type == 0) ? ( (n_elem > 1) ? eT(n_elem-1) : eT(1) ) : eT(n_elem);
00052   const eT var_val  = (acc2 - acc3*acc3/div_val) / norm_val;
00053   
00054   return var_val;
00055   }
00056 
00057 
00058 
00059 //! find the variance of an array (version for complex numbers)
00060 template<typename T>
00061 inline
00062 T
00063 op_var::direct_var(const std::complex<T>* const X, const u32 n_elem, const u32 norm_type)
00064   {
00065   arma_extra_debug_sigprint();
00066   
00067   typedef typename std::complex<T> eT;
00068   
00069   eT acc1 = eT(0);
00070   
00071   for(u32 i=0; i<n_elem; ++i)
00072     {
00073     acc1 += X[i];
00074     }
00075   
00076   const T div_val = (n_elem > 0) ? T(n_elem) : T(1);
00077   acc1 /= div_val;
00078   
00079   T  acc2 =  T(0);
00080   eT acc3 = eT(0);
00081   
00082   for(u32 i=0; i<n_elem; ++i)
00083     {
00084     const eT tmp = acc1 - X[i];
00085     
00086     acc2 += std::norm(tmp);
00087     acc3 += tmp;
00088     }
00089   
00090   const T norm_val = (norm_type == 0) ? ( (n_elem > 1) ? T(n_elem-1) : T(1) ) : T(n_elem);
00091   const T var_val  = (acc2 - std::norm(acc3)/div_val) / norm_val;
00092   
00093   return var_val;
00094   }
00095 
00096 
00097 
00098 //! find the variance of a subview_row
00099 template<typename eT>
00100 inline 
00101 typename get_pod_type<eT>::result
00102 op_var::direct_var(const subview_row<eT>& X, const u32 norm_type)
00103   {
00104   arma_extra_debug_sigprint();
00105   
00106   const u32 n_elem = X.n_elem;
00107   
00108   podarray<eT> tmp(n_elem);
00109   
00110   eT* tmp_mem = tmp.memptr();
00111   
00112   for(u32 i=0; i<n_elem; ++i)
00113     {
00114     tmp_mem[i] = X[i];
00115     }
00116   
00117   return op_var::direct_var(tmp_mem, n_elem, norm_type);
00118   }
00119 
00120 
00121 
00122 //! find the variance of a subview_col
00123 template<typename eT>
00124 inline 
00125 typename get_pod_type<eT>::result
00126 op_var::direct_var(const subview_col<eT>& X, const u32 norm_type)
00127   {
00128   arma_extra_debug_sigprint();
00129   
00130   return op_var::direct_var(X.colptr(0), X.n_elem, norm_type);
00131   }
00132 
00133 
00134 
00135 //! find the variance of a diagview
00136 template<typename eT>
00137 inline 
00138 typename get_pod_type<eT>::result
00139 op_var::direct_var(const diagview<eT>& X, const u32 norm_type)
00140   {
00141   arma_extra_debug_sigprint();
00142   
00143   const u32 n_elem = X.n_elem;
00144   
00145   podarray<eT> tmp(n_elem);
00146   
00147   eT* tmp_mem = tmp.memptr();
00148   
00149   for(u32 i=0; i<n_elem; ++i)
00150     {
00151     tmp_mem[i] = X[i];
00152     }
00153   
00154   return op_var::direct_var(tmp_mem, n_elem, norm_type);
00155   }
00156 
00157 
00158 
00159 //! \brief
00160 //! For each row or for each column, find the variance.
00161 //! The result is stored in a dense matrix that has either one column or one row.
00162 //! The dimension, for which the variances are found, is set via the var() function.
00163 template<typename eT>
00164 inline
00165 void
00166 op_var::apply(Mat< typename get_pod_type<eT>::result >& out, const Mat<eT>& X, const u32 norm_type, const u32 dim)
00167   {
00168   arma_extra_debug_sigprint();
00169   
00170   arma_debug_check( (X.n_elem == 0), "var(): given matrix has no elements" );
00171   
00172   arma_debug_check( (norm_type > 1), "var(): incorrect usage. norm_type must be 0 or 1");
00173   arma_debug_check( (dim > 1),       "var(): incorrect usage. dim must be 0 or 1"      );
00174   
00175   
00176   if(dim == 0)
00177     {
00178     arma_extra_debug_print("op_var::apply(), dim = 0");
00179     
00180     out.set_size(1, X.n_cols);
00181     
00182     for(u32 col=0; col<X.n_cols; ++col)
00183       {
00184       out[col] = op_var::direct_var( X.colptr(col), X.n_rows, norm_type );
00185       }
00186     }
00187   else
00188   if(dim == 1)
00189     {
00190     arma_extra_debug_print("op_var::apply(), dim = 1");
00191     
00192     const u32 n_rows = X.n_rows;
00193     const u32 n_cols = X.n_cols;
00194     
00195     out.set_size(n_rows, 1);
00196     
00197     podarray<eT> tmp(n_cols);
00198     
00199     eT* tmp_mem = tmp.memptr();
00200     
00201     for(u32 row=0; row<n_rows; ++row)
00202       {
00203       for(u32 col=0; col<n_cols; ++col)
00204         {
00205         tmp_mem[col] = X.at(row,col);
00206         }
00207       
00208       out[row] = op_var::direct_var(tmp_mem, n_cols, norm_type);
00209       }
00210     
00211     }
00212   
00213   }
00214 
00215 
00216 
00217 //! @}