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