op_median_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_median
00018 //! @{
00019 
00020 
00021 
00022 //! find the median value of a std::vector (contents is modified)
00023 template<typename eT>
00024 inline 
00025 eT
00026 op_median::direct_median(std::vector<eT>& X)
00027   {
00028   arma_extra_debug_sigprint();
00029   
00030   std::sort(X.begin(), X.end());
00031   
00032   const u32 n_elem = X.size();
00033   const u32 half   = n_elem/2;
00034   
00035   if((n_elem % 2) == 0)
00036     {
00037     return (X[half-1] + X[half]) / eT(2);
00038     }
00039   else
00040     {
00041     return X[half];
00042     }
00043   }
00044 
00045 
00046 
00047 template<typename eT>
00048 inline 
00049 eT
00050 op_median::direct_median(const eT* X, const u32 n_elem)
00051   {
00052   arma_extra_debug_sigprint();
00053   
00054   std::vector<eT> tmp(X, X+n_elem);
00055   return op_median::direct_median(tmp);
00056   }
00057 
00058 
00059 
00060 template<typename eT>
00061 inline 
00062 eT
00063 op_median::direct_median(const subview<eT>& X)
00064   {
00065   arma_extra_debug_sigprint();
00066   
00067   std::vector<eT> tmp(X.n_elem);
00068   
00069   for(u32 i=0; i<X.n_elem; ++i)
00070     {
00071     tmp[i] = X[i];
00072     }
00073   
00074   return op_median::direct_median(tmp);
00075   }
00076 
00077 
00078 
00079 template<typename eT>
00080 inline 
00081 eT
00082 op_median::direct_median(const diagview<eT>& X)
00083   {
00084   arma_extra_debug_sigprint();
00085   
00086   std::vector<eT> tmp(X.n_elem);
00087   
00088   for(u32 i=0; i<X.n_elem; ++i)
00089     {
00090     tmp[i] = X[i];
00091     }
00092   
00093   return op_median::direct_median(tmp);
00094   }
00095 
00096 
00097 
00098 //! \brief
00099 //! For each row or for each column, find the median value.
00100 //! The result is stored in a dense matrix that has either one column or one row.
00101 //! The dimension, for which the medians are found, is set via the median() function.
00102 template<typename T1>
00103 inline
00104 void
00105 op_median::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_median>& in)
00106   {
00107   arma_extra_debug_sigprint();
00108   
00109   typedef typename T1::elem_type eT;
00110   
00111   const unwrap_check<T1> tmp(in.m, out);
00112   const Mat<eT>& X = tmp.M;
00113   
00114   arma_debug_check( (X.n_elem == 0), "median(): given matrix has no elements" );
00115   
00116   const u32 dim = in.aux_u32_a;
00117   arma_debug_check( (dim > 1), "median(): incorrect usage. dim must be 0 or 1");
00118   
00119   
00120   if(dim == 0)  // column-wise
00121     {
00122     arma_extra_debug_print("op_median::apply(), dim = 0");
00123     
00124     out.set_size(1, X.n_cols);
00125     
00126     std::vector<eT> tmp_vec(X.n_rows);
00127     
00128     for(u32 col=0; col<X.n_cols; ++col)
00129       {
00130       const eT* colmem = X.colptr(col);
00131       
00132       for(u32 row=0; row<X.n_rows; ++row)
00133         {
00134         tmp_vec[row] = colmem[row];
00135         }
00136       
00137       out[col] = op_median::direct_median(tmp_vec);
00138       }
00139     }
00140   else
00141   if(dim == 1)  // row-wise
00142     {
00143     arma_extra_debug_print("op_median::apply(), dim = 1");
00144   
00145     out.set_size(X.n_rows, 1);
00146     
00147     std::vector<eT> tmp_vec(X.n_cols);
00148     
00149     for(u32 row=0; row<X.n_rows; ++row)
00150       {
00151       for(u32 col=0; col<X.n_cols; ++col)
00152         {
00153         tmp_vec[col] = X.at(row,col);
00154         }
00155   
00156       out[row] = op_median::direct_median(tmp_vec);
00157       }
00158     }
00159   
00160   }
00161 
00162 
00163 
00164 template<typename T>
00165 inline 
00166 void
00167 op_median::direct_cx_median_index(u32& out_index1, u32& out_index2, std::vector< arma_cx_median_packet<T> >& X)
00168   {
00169   arma_extra_debug_sigprint();
00170   
00171   std::sort(X.begin(), X.end());
00172   
00173   const u32 n_elem = X.size();
00174   const u32 half   = n_elem/2;
00175   
00176   if((n_elem % 2) == 0)
00177     {
00178     out_index1 = X[half-1].index;
00179     out_index2 = X[half].index;
00180     }
00181   else
00182     {
00183     out_index1 = X[half].index;
00184     out_index2 = X[half].index;
00185     }
00186   
00187   }
00188 
00189 
00190 
00191 template<typename T>
00192 inline 
00193 void
00194 op_median::direct_cx_median_index(u32& out_index1, u32& out_index2, const std::complex<T>* X, const u32 n_elem)
00195   {
00196   arma_extra_debug_sigprint();
00197   
00198   std::vector< arma_cx_median_packet<T> > tmp(n_elem);
00199   
00200   for(u32 i=0; i<n_elem; ++i)
00201     {
00202     tmp[i].val   = std::abs(X[i]);
00203     tmp[i].index = i;
00204     }
00205   
00206   op_median::direct_cx_median_index(out_index1, out_index2, tmp);
00207   }
00208 
00209 
00210 
00211 template<typename T>
00212 inline 
00213 void
00214 op_median::direct_cx_median_index(u32& out_index1, u32& out_index2, const subview< std::complex<T> >&X)
00215   {
00216   arma_extra_debug_sigprint();
00217   
00218   const u32 n_elem = X.n_elem;
00219   
00220   std::vector< arma_cx_median_packet<T> > tmp(n_elem);
00221   
00222   for(u32 i=0; i<n_elem; ++i)
00223     {
00224     tmp[i].val   = std::abs(X[i]);
00225     tmp[i].index = i;
00226     }
00227   
00228   op_median::direct_cx_median_index(out_index1, out_index2, tmp);
00229   }
00230 
00231 
00232 
00233 template<typename T>
00234 inline 
00235 void
00236 op_median::direct_cx_median_index(u32& out_index1, u32& out_index2, const diagview< std::complex<T> >&X)
00237   {
00238   arma_extra_debug_sigprint();
00239   
00240   const u32 n_elem = X.n_elem;
00241   
00242   std::vector< arma_cx_median_packet<T> > tmp(n_elem);
00243   
00244   for(u32 i=0; i<n_elem; ++i)
00245     {
00246     tmp[i].val   = std::abs(X[i]);
00247     tmp[i].index = i;
00248     }
00249   
00250   op_median::direct_cx_median_index(out_index1, out_index2, tmp);
00251   }
00252 
00253 
00254 
00255 //! Implementation for complex numbers
00256 template<typename T, typename T1>
00257 inline
00258 void
00259 op_median::apply(Mat< std::complex<T> >& out, const Op<T1,op_median>& in)
00260   {
00261   arma_extra_debug_sigprint();
00262   
00263   typedef typename std::complex<T> eT;
00264   isnt_same_type<eT, typename T1::elem_type>::check();
00265   
00266   const unwrap_check<T1> tmp(in.m, out);
00267   const Mat<eT>& X = tmp.M;
00268   
00269   arma_debug_check( (X.n_elem == 0), "median(): given matrix has no elements" );
00270   
00271   const u32 dim = in.aux_u32_a;
00272   arma_debug_check( (dim > 1), "median(): incorrect usage. dim must be 0 or 1");
00273   
00274   
00275   if(dim == 0)  // column-wise
00276     {
00277     arma_extra_debug_print("op_median::apply(), dim = 0");
00278     
00279     out.set_size(1, X.n_cols);
00280     
00281     std::vector< arma_cx_median_packet<T> > tmp_vec(X.n_rows);
00282     
00283     for(u32 col=0; col<X.n_cols; ++col)
00284       {
00285       const eT* colmem = X.colptr(col);
00286       
00287       for(u32 row=0; row<X.n_rows; ++row)
00288         {
00289         tmp_vec[row].val   = std::abs(colmem[row]);
00290         tmp_vec[row].index = row;
00291         }
00292       
00293       u32 index1;
00294       u32 index2;
00295       op_median::direct_cx_median_index(index1, index2, tmp_vec);
00296       
00297       out[col] = (colmem[index1] + colmem[index2]) / T(2);
00298       }
00299     }
00300   else
00301   if(dim == 1)  // row-wise
00302     {
00303     arma_extra_debug_print("op_median::apply(), dim = 1");
00304   
00305     out.set_size(X.n_rows, 1);
00306     
00307     std::vector< arma_cx_median_packet<T> > tmp_vec(X.n_cols);
00308     
00309     for(u32 row=0; row<X.n_rows; ++row)
00310       {
00311       for(u32 col=0; col<X.n_cols; ++col)
00312         {
00313         tmp_vec[col].val   = std::abs(X.at(row,col));
00314         tmp_vec[row].index = col;
00315         }
00316   
00317       u32 index1;
00318       u32 index2;
00319       op_median::direct_cx_median_index(index1, index2, tmp_vec);
00320       
00321       out[row] = ( X.at(row,index1) + X.at(row,index2) ) / T(2);
00322       }
00323     }
00324   
00325   }
00326 
00327 
00328 
00329 //! @}