op_var_meat.hpp

Go to the documentation of this file.
00001 // Copyright (C) 2009 NICTA
00002 // 
00003 // Authors:
00004 // - Conrad Sanderson (conradsand at ieee dot org)
00005 // 
00006 // This file is part of the Armadillo C++ library.
00007 // It is provided without any warranty of fitness
00008 // for any purpose. You can redistribute this file
00009 // and/or modify it under the terms of the GNU
00010 // Lesser General Public License (LGPL) as published
00011 // by the Free Software Foundation, either version 3
00012 // of the License or (at your option) any later version.
00013 // (see http://www.opensource.org/licenses for more info)
00014 
00015 
00016 //! \addtogroup op_var
00017 //! @{
00018 
00019 
00020 //! find the variance of an array
00021 template<typename eT>
00022 inline 
00023 eT
00024 op_var::direct_var(const eT* const X, const u32 n_elem, const u32 norm_type)
00025   {
00026   arma_extra_debug_sigprint();
00027   
00028   eT acc1 = eT(0);
00029   eT acc2 = eT(0);
00030   
00031   for(u32 i=0; i<n_elem; ++i)
00032     {
00033     const eT tmp_val = X[i];
00034     acc1 += tmp_val;
00035     acc2 += tmp_val*tmp_val;
00036     }
00037   
00038   const eT norm_val = (norm_type == 0) ? eT(n_elem-1) : eT(n_elem);
00039   const eT var_val  = (acc2 - acc1*acc1/eT(n_elem)) / norm_val;
00040   
00041   return var_val;
00042   }
00043 
00044 
00045 
00046 //! find the variance of an array (version for complex numbers)
00047 template<typename T>
00048 inline 
00049 T
00050 op_var::direct_var(const std::complex<T>* const X, const u32 n_elem, const u32 norm_type)
00051   {
00052   arma_extra_debug_sigprint();
00053   
00054   typedef typename std::complex<T> eT;
00055   
00056   eT acc1 = eT(0);
00057   T  acc2 = T(0);
00058   
00059   for(u32 i=0; i<n_elem; ++i)
00060     {
00061     acc1 += X[i];
00062     acc2 += std::norm(X[i]);
00063     }
00064   
00065   const T norm_val = (norm_type == 0) ? T(n_elem-1) : T(n_elem);
00066   const T var_val  = (acc2 - std::norm(acc1)/T(n_elem)) / norm_val;
00067   
00068   return var_val;
00069   }
00070 
00071 
00072 
00073 //! \brief
00074 //! For each row or for each column, find the variance.
00075 //! The result is stored in a dense matrix that has either one column or one row.
00076 //! The dimension, for which the variances are found, is set via the var() function.
00077 template<typename eT>
00078 inline
00079 void
00080 op_var::apply(Mat<eT>& out, const Mat<eT>& X, const u32 norm_type, const u32 dim)
00081   {
00082   arma_extra_debug_sigprint();
00083   
00084   arma_debug_check( (X.n_elem == 0), "op_var::apply(): given matrix has no elements" );
00085   
00086   arma_debug_check( (norm_type > 1), "op_var::apply(): incorrect usage. norm_type must be 0 or 1");
00087   arma_debug_check( (dim > 1),       "op_var::apply(): incorrect usage. dim must be 0 or 1"      );
00088   
00089   
00090   if(dim == 0)
00091     {
00092     arma_extra_debug_print("op_var::apply(), dim = 0");
00093     
00094     out.set_size(1, X.n_cols);
00095     
00096     for(u32 col=0; col<X.n_cols; ++col)
00097       {
00098       out[col] = op_var::direct_var( X.colptr(col), X.n_rows, norm_type );
00099       }
00100     }
00101   else
00102   if(dim == 1)
00103     {
00104     arma_extra_debug_print("op_var::apply(), dim = 1");
00105     
00106     out.set_size(X.n_rows, 1);
00107     
00108     const eT norm_val = (norm_type == 0) ? eT(X.n_cols-1) : eT(X.n_cols);
00109     
00110     for(u32 row=0; row<X.n_rows; ++row)
00111       {
00112       eT acc1 = eT(0);
00113       eT acc2 = eT(0);
00114   
00115       for(u32 col=0; col<X.n_cols; ++col)
00116         {
00117         const eT tmp_val = X.at(row,col);
00118         acc1 += tmp_val;
00119         acc2 += tmp_val*tmp_val;
00120         }
00121       
00122       const eT var_val = (acc2 - acc1*acc1/eT(X.n_cols)) / norm_val;
00123       
00124       out[row] = var_val;
00125       }
00126     
00127     }
00128   
00129   }
00130 
00131 
00132 
00133 //! implementation for complex numbers
00134 template<typename T>
00135 inline
00136 void
00137 op_var::apply(Mat<T>& out, const Mat< std::complex<T> >& X, const u32 norm_type, const u32 dim)
00138   {
00139   arma_extra_debug_sigprint();
00140   
00141   typedef typename std::complex<T> eT;
00142   
00143   arma_debug_check( (X.n_elem == 0), "op_var::apply(): given matrix has no elements" );
00144   
00145   arma_debug_check( (norm_type > 1), "op_var::apply(): incorrect usage. norm_type must be 0 or 1");
00146   arma_debug_check( (dim > 1),       "op_var::apply(): incorrect usage. dim must be 0 or 1"      );
00147   
00148   
00149   if(dim == 0)
00150     {
00151     arma_extra_debug_print("op_var::apply(), dim = 0");
00152     
00153     out.set_size(1, X.n_cols);
00154     
00155     for(u32 col=0; col<X.n_cols; ++col)
00156       {
00157       out[col] = op_var::direct_var( X.colptr(col), X.n_rows, norm_type );
00158       }
00159     }
00160   else
00161   if(dim == 1)
00162     {
00163     arma_extra_debug_print("op_var::apply(), dim = 1");
00164     
00165     out.set_size(X.n_rows, 1);
00166     
00167     const T norm_val = (norm_type == 0) ? T(X.n_cols-1) : T(X.n_cols);
00168     
00169     for(u32 row=0; row<X.n_rows; ++row)
00170       {
00171       eT acc1 = eT(0);
00172       T  acc2 = T(0);
00173   
00174       for(u32 col=0; col<X.n_cols; ++col)
00175         {
00176         acc1 += X.at(row,col);
00177         acc2 += std::norm(X.at(row,col));
00178         }
00179       
00180       const T var_val = (acc2 - std::norm(acc1)/T(X.n_cols)) / norm_val;
00181       
00182       out[row] = var_val;
00183       }
00184     
00185     }
00186   
00187   }
00188 
00189 
00190 
00191 //! find the variance of a subview
00192 template<typename eT>
00193 inline 
00194 eT
00195 op_var::direct_var(const subview<eT>& X, const u32 norm_type)
00196   {
00197   arma_extra_debug_sigprint();
00198   
00199   const u32 n_elem = X.n_elem;
00200   
00201   eT acc1 = eT(0);
00202   eT acc2 = eT(0);
00203   
00204   for(u32 i=0; i<n_elem; ++i)
00205     {
00206     const eT tmp_val = X[i];
00207     acc1 += tmp_val;
00208     acc2 += tmp_val*tmp_val;
00209     }
00210   
00211   const eT norm_val = (norm_type == 0) ? eT(n_elem-1) : eT(n_elem);
00212   const eT var_val  = (acc2 - acc1*acc1/eT(n_elem)) / norm_val;
00213   
00214   return var_val;
00215   }
00216 
00217 
00218 
00219 //! find the variance of a subview (version for complex numbers)
00220 template<typename T>
00221 inline 
00222 T
00223 op_var::direct_var(const subview< std::complex<T> >& X, const u32 norm_type)
00224   {
00225   arma_extra_debug_sigprint();
00226   
00227   typedef typename std::complex<T> eT;
00228   
00229   const u32 n_elem = X.n_elem;
00230   
00231   eT acc1 = eT(0);
00232   T  acc2 = T(0);
00233   
00234   for(u32 i=0; i<n_elem; ++i)
00235     {
00236     acc1 += X[i];
00237     acc2 += std::norm(X[i]);
00238     }
00239   
00240   const T norm_val = (norm_type == 0) ? T(n_elem-1) : T(n_elem);
00241   const T var_val  = (acc2 - std::norm(acc1)/T(n_elem)) / norm_val;
00242   
00243   return var_val;
00244   }
00245 
00246 
00247 
00248 //! find the variance of a diagview
00249 template<typename eT>
00250 inline 
00251 eT
00252 op_var::direct_var(const diagview<eT>& X, const u32 norm_type)
00253   {
00254   arma_extra_debug_sigprint();
00255   
00256   const u32 n_elem = X.n_elem;
00257   
00258   eT acc1 = eT(0);
00259   eT acc2 = eT(0);
00260   
00261   for(u32 i=0; i<n_elem; ++i)
00262     {
00263     const eT tmp_val = X[i];
00264     acc1 += tmp_val;
00265     acc2 += tmp_val*tmp_val;
00266     }
00267   
00268   const eT norm_val = (norm_type == 0) ? eT(n_elem-1) : eT(n_elem);
00269   const eT var_val  = (acc2 - acc1*acc1/eT(n_elem)) / norm_val;
00270   
00271   return var_val;
00272   }
00273 
00274 
00275 
00276 //! find the variance of a diagview (version for complex numbers)
00277 template<typename T>
00278 inline 
00279 T
00280 op_var::direct_var(const diagview< std::complex<T> >& X, const u32 norm_type)
00281   {
00282   arma_extra_debug_sigprint();
00283   
00284   typedef typename std::complex<T> eT;
00285   
00286   const u32 n_elem = X.n_elem;
00287   
00288   eT acc1 = eT(0);
00289   T  acc2 = T(0);
00290   
00291   for(u32 i=0; i<n_elem; ++i)
00292     {
00293     acc1 += X[i];
00294     acc2 += std::norm(X[i]);
00295     }
00296   
00297   const T norm_val = (norm_type == 0) ? T(n_elem-1) : T(n_elem);
00298   const T var_val  = (acc2 - std::norm(acc1)/T(n_elem)) / norm_val;
00299   
00300   return var_val;
00301   }
00302 
00303 
00304 
00305 //! @}