op_stddev_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_stddev
00017 //! @{
00018 
00019 
00020 //! \brief
00021 //! For each row or for each column, find the standard deviation.
00022 //! The result is stored in a dense matrix that has either one column or one row.
00023 //! The dimension for which the standard deviations are found is set via the stddev() function.
00024 template<typename eT>
00025 inline
00026 void
00027 op_stddev::apply(Mat<eT>& out, const Mat<eT>& X, const u32 norm_type, const u32 dim)
00028   {
00029   arma_extra_debug_sigprint();
00030   
00031   arma_debug_check( (X.n_elem == 0), "op_stddev::apply(): given matrix has no elements" );
00032   
00033   arma_debug_check( (norm_type > 1), "op_stddev::apply(): incorrect usage. norm_type must be 0 or 1");
00034   arma_debug_check( (dim > 1),       "op_stddev::apply(): incorrect usage. dim must be 0 or 1"      );
00035   
00036   if(dim == 0)
00037     {
00038     arma_extra_debug_print("op_stddev::apply(), dim = 0");
00039     
00040     out.set_size(1, X.n_cols);
00041     
00042     for(u32 col=0; col<X.n_cols; ++col)
00043       {
00044       out[col] = std::sqrt( op_var::direct_var( X.colptr(col), X.n_rows, norm_type ) );
00045       }
00046     }
00047   else
00048   if(dim == 1)
00049     {
00050     arma_extra_debug_print("op_stddev::apply(), dim = 1");
00051     
00052     out.set_size(X.n_rows, 1);
00053     
00054     const eT norm_val = (norm_type == 0) ? eT(X.n_cols-1) : eT(X.n_cols);
00055     
00056     for(u32 row=0; row<X.n_rows; ++row)
00057       {
00058       eT acc1 = eT(0);
00059       eT acc2 = eT(0);
00060   
00061       for(u32 col=0; col<X.n_cols; ++col)
00062         {
00063         const eT tmp_val = X.at(row,col);
00064         acc1 += tmp_val;
00065         acc2 += tmp_val*tmp_val;
00066         }
00067       
00068       const eT sd_val = std::sqrt( (acc2 - acc1*acc1/eT(X.n_cols)) / norm_val );
00069       
00070       out[row] = sd_val;
00071       }
00072     
00073     }
00074   
00075   }
00076 
00077 
00078 
00079 //! implementation for complex numbers
00080 template<typename T>
00081 inline
00082 void
00083 op_stddev::apply(Mat<T>& out, const Mat< std::complex<T> >& X, const u32 norm_type, const u32 dim)
00084   {
00085   arma_extra_debug_sigprint();
00086   
00087   typedef typename std::complex<T> eT;
00088   
00089   arma_debug_check( (X.n_elem == 0), "op_stddev::apply(): given matrix has no elements" );
00090   
00091   arma_debug_check( (norm_type > 1), "op_stddev::apply(): incorrect usage. norm_type must be 0 or 1");
00092   arma_debug_check( (dim > 1),       "op_stddev::apply(): incorrect usage. dim must be 0 or 1"      );
00093   
00094   
00095   if(dim == 0)
00096     {
00097     arma_extra_debug_print("op_stddev::apply(), dim = 0");
00098     
00099     out.set_size(1, X.n_cols);
00100     
00101     for(u32 col=0; col<X.n_cols; ++col)
00102       {
00103       out[col] = std::sqrt( op_var::direct_var( X.colptr(col), X.n_rows, norm_type ) );
00104       }
00105     }
00106   else
00107   if(dim == 1)
00108     {
00109     arma_extra_debug_print("op_stddev::apply(), dim = 1");
00110     
00111     out.set_size(X.n_rows, 1);
00112     
00113     const T norm_val = (norm_type == 0) ? T(X.n_cols-1) : T(X.n_cols);
00114     
00115     for(u32 row=0; row<X.n_rows; ++row)
00116       {
00117       eT acc1 = eT(0);
00118       T  acc2 = T(0);
00119   
00120       for(u32 col=0; col<X.n_cols; ++col)
00121         {
00122         acc1 += X.at(row,col);
00123         acc2 += std::norm(X.at(row,col));
00124         }
00125       
00126       const T var_val = (acc2 - std::norm(acc1)/T(X.n_cols)) / norm_val;
00127       
00128       out[row] = std::sqrt(var_val);
00129       }
00130     
00131     }
00132   
00133   }
00134 
00135 
00136 
00137 //! @}