Chris@49: // Copyright (C) 2009-2013 NICTA (www.nicta.com.au)
Chris@49: // Copyright (C) 2009-2013 Conrad Sanderson
Chris@49: // Copyright (C) 2013 Ruslan Shestopalyuk
Chris@49: // 
Chris@49: // This Source Code Form is subject to the terms of the Mozilla Public
Chris@49: // License, v. 2.0. If a copy of the MPL was not distributed with this
Chris@49: // file, You can obtain one at http://mozilla.org/MPL/2.0/.
Chris@49: 
Chris@49: 
Chris@49: //! \addtogroup op_median
Chris@49: //! @{
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: //! \brief
Chris@49: //! For each row or for each column, find the median value.
Chris@49: //! The result is stored in a dense matrix that has either one column or one row.
Chris@49: //! The dimension, for which the medians are found, is set via the median() function.
Chris@49: template<typename T1>
Chris@49: inline
Chris@49: void
Chris@49: op_median::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_median>& in)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   typedef typename T1::elem_type eT;
Chris@49:   
Chris@49:   const uword dim = in.aux_uword_a;
Chris@49:   arma_debug_check( (dim > 1), "median(): incorrect usage. dim must be 0 or 1");
Chris@49:   
Chris@49:   const Proxy<T1> P(in.m);
Chris@49:   
Chris@49:   typedef typename Proxy<T1>::stored_type P_stored_type;
Chris@49:   
Chris@49:   const bool is_alias = P.is_alias(out);
Chris@49:   
Chris@49:   if( (is_Mat<P_stored_type>::value == true) || is_alias )
Chris@49:     {
Chris@49:     const unwrap_check<P_stored_type> tmp(P.Q, is_alias);
Chris@49:     
Chris@49:     const typename unwrap_check<P_stored_type>::stored_type& X = tmp.M;
Chris@49:     
Chris@49:     const uword X_n_rows = X.n_rows;
Chris@49:     const uword X_n_cols = X.n_cols;
Chris@49:     
Chris@49:     if(dim == 0)  // in each column
Chris@49:       {
Chris@49:       arma_extra_debug_print("op_median::apply(), dim = 0");
Chris@49:       
Chris@49:       arma_debug_check( (X_n_rows == 0), "median(): given object has zero rows" );
Chris@49:       
Chris@49:       out.set_size(1, X_n_cols);
Chris@49:       
Chris@49:       std::vector<eT> tmp_vec(X_n_rows);
Chris@49:       
Chris@49:       for(uword col=0; col < X_n_cols; ++col)
Chris@49:         {
Chris@49:         arrayops::copy( &(tmp_vec[0]), X.colptr(col), X_n_rows );
Chris@49:         
Chris@49:         out[col] = op_median::direct_median(tmp_vec);
Chris@49:         }
Chris@49:       }
Chris@49:     else  // in each row
Chris@49:       {
Chris@49:       arma_extra_debug_print("op_median::apply(), dim = 1");
Chris@49:       
Chris@49:       arma_debug_check( (X_n_cols == 0), "median(): given object has zero columns" );
Chris@49:       
Chris@49:       out.set_size(X_n_rows, 1);
Chris@49:       
Chris@49:       std::vector<eT> tmp_vec(X_n_cols);
Chris@49:         
Chris@49:       for(uword row=0; row < X_n_rows; ++row)
Chris@49:         {
Chris@49:         for(uword col=0; col < X_n_cols; ++col)  { tmp_vec[col] = X.at(row,col); }
Chris@49:         
Chris@49:         out[row] = op_median::direct_median(tmp_vec);
Chris@49:         }
Chris@49:       }
Chris@49:     }
Chris@49:   else
Chris@49:     {
Chris@49:     const uword P_n_rows = P.get_n_rows();
Chris@49:     const uword P_n_cols = P.get_n_cols();
Chris@49:     
Chris@49:     if(dim == 0)  // in each column
Chris@49:       {
Chris@49:       arma_extra_debug_print("op_median::apply(), dim = 0");
Chris@49:       
Chris@49:       arma_debug_check( (P_n_rows == 0), "median(): given object has zero rows" );
Chris@49:       
Chris@49:       out.set_size(1, P_n_cols);
Chris@49:       
Chris@49:       std::vector<eT> tmp_vec(P_n_rows);
Chris@49:       
Chris@49:       for(uword col=0; col < P_n_cols; ++col)
Chris@49:         {
Chris@49:         for(uword row=0; row < P_n_rows; ++row)  { tmp_vec[row] = P.at(row,col); }
Chris@49:         
Chris@49:         out[col] = op_median::direct_median(tmp_vec);
Chris@49:         }
Chris@49:       }
Chris@49:     else  // in each row
Chris@49:       {
Chris@49:       arma_extra_debug_print("op_median::apply(), dim = 1");
Chris@49:       
Chris@49:       arma_debug_check( (P_n_cols == 0), "median(): given object has zero columns" );
Chris@49:       
Chris@49:       out.set_size(P_n_rows, 1);
Chris@49:       
Chris@49:       std::vector<eT> tmp_vec(P_n_cols);
Chris@49:         
Chris@49:       for(uword row=0; row < P_n_rows; ++row)
Chris@49:         {
Chris@49:         for(uword col=0; col < P_n_cols; ++col)  { tmp_vec[col] = P.at(row,col); }
Chris@49:         
Chris@49:         out[row] = op_median::direct_median(tmp_vec);
Chris@49:         }
Chris@49:       }
Chris@49:     }
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: //! Implementation for complex numbers
Chris@49: template<typename T, typename T1>
Chris@49: inline
Chris@49: void
Chris@49: op_median::apply(Mat< std::complex<T> >& out, const Op<T1,op_median>& in)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   typedef typename std::complex<T> eT;
Chris@49:   
Chris@49:   arma_type_check(( is_same_type<eT, typename T1::elem_type>::value == false ));
Chris@49:   
Chris@49:   const unwrap_check<T1> tmp(in.m, out);
Chris@49:   const Mat<eT>&     X = tmp.M;
Chris@49:   
Chris@49:   const uword X_n_rows = X.n_rows;
Chris@49:   const uword X_n_cols = X.n_cols;
Chris@49:   
Chris@49:   const uword dim = in.aux_uword_a;
Chris@49:   arma_debug_check( (dim > 1), "median(): incorrect usage. dim must be 0 or 1");
Chris@49:   
Chris@49:   if(dim == 0)  // in each column
Chris@49:     {
Chris@49:     arma_extra_debug_print("op_median::apply(), dim = 0");
Chris@49:     
Chris@49:     arma_debug_check( (X_n_rows == 0), "median(): given object has zero rows" );
Chris@49: 
Chris@49:     out.set_size(1, X_n_cols);
Chris@49:     
Chris@49:     std::vector< arma_cx_median_packet<T> > tmp_vec(X_n_rows);
Chris@49:     
Chris@49:     for(uword col=0; col<X_n_cols; ++col)
Chris@49:       {
Chris@49:       const eT* colmem = X.colptr(col);
Chris@49:       
Chris@49:       for(uword row=0; row<X_n_rows; ++row)
Chris@49:         {
Chris@49:         tmp_vec[row].val   = std::abs(colmem[row]);
Chris@49:         tmp_vec[row].index = row;
Chris@49:         }
Chris@49:       
Chris@49:       uword index1;
Chris@49:       uword index2;
Chris@49:       op_median::direct_cx_median_index(index1, index2, tmp_vec);
Chris@49:         
Chris@49:       out[col] = op_mean::robust_mean(colmem[index1], colmem[index2]);
Chris@49:       }
Chris@49:     }
Chris@49:   else
Chris@49:   if(dim == 1)  // in each row
Chris@49:     {
Chris@49:     arma_extra_debug_print("op_median::apply(), dim = 1");
Chris@49:     
Chris@49:     arma_debug_check( (X_n_cols == 0), "median(): given object has zero columns" );
Chris@49:     
Chris@49:     out.set_size(X_n_rows, 1);
Chris@49:     
Chris@49:     std::vector< arma_cx_median_packet<T> > tmp_vec(X_n_cols);
Chris@49:     
Chris@49:     for(uword row=0; row<X_n_rows; ++row)
Chris@49:       {
Chris@49:       for(uword col=0; col<X_n_cols; ++col)
Chris@49:         {
Chris@49:         tmp_vec[col].val   = std::abs(X.at(row,col));
Chris@49:         tmp_vec[col].index = col;
Chris@49:         }
Chris@49:       
Chris@49:       uword index1;
Chris@49:       uword index2;
Chris@49:       op_median::direct_cx_median_index(index1, index2, tmp_vec);
Chris@49:       
Chris@49:       out[row] = op_mean::robust_mean( X.at(row,index1), X.at(row,index2) );
Chris@49:       }
Chris@49:     }
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename T1>
Chris@49: inline
Chris@49: typename T1::elem_type
Chris@49: op_median::median_vec
Chris@49:   (
Chris@49:   const T1& X,
Chris@49:   const typename arma_not_cx<typename T1::elem_type>::result* junk
Chris@49:   )
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   arma_ignore(junk);
Chris@49:   
Chris@49:   typedef typename T1::elem_type eT;
Chris@49:   
Chris@49:   typedef typename Proxy<T1>::stored_type P_stored_type;
Chris@49:     
Chris@49:   const Proxy<T1> P(X);
Chris@49:   
Chris@49:   const uword n_elem = P.get_n_elem();
Chris@49:   
Chris@49:   arma_debug_check( (n_elem == 0), "median(): given object has no elements" );
Chris@49:   
Chris@49:   std::vector<eT> tmp_vec(n_elem);
Chris@49:   
Chris@49:   if(is_Mat<P_stored_type>::value == true)
Chris@49:     {
Chris@49:     const unwrap<P_stored_type> tmp(P.Q);
Chris@49:     
Chris@49:     const typename unwrap_check<P_stored_type>::stored_type& Y = tmp.M;
Chris@49:     
Chris@49:     arrayops::copy( &(tmp_vec[0]), Y.memptr(), n_elem );
Chris@49:     }
Chris@49:   else
Chris@49:     {
Chris@49:     if(Proxy<T1>::prefer_at_accessor == false)
Chris@49:       {
Chris@49:       typedef typename Proxy<T1>::ea_type ea_type;
Chris@49:       
Chris@49:       ea_type A = P.get_ea();
Chris@49:       
Chris@49:       for(uword i=0; i<n_elem; ++i)  { tmp_vec[i] = A[i]; }
Chris@49:       }
Chris@49:     else
Chris@49:       {
Chris@49:       const uword n_rows = P.get_n_rows();
Chris@49:       const uword n_cols = P.get_n_cols();
Chris@49:       
Chris@49:       if(n_cols == 1)
Chris@49:         {
Chris@49:         for(uword row=0; row < n_rows; ++row)  { tmp_vec[row] = P.at(row,0); }
Chris@49:         }
Chris@49:       else
Chris@49:       if(n_rows == 1)
Chris@49:         {
Chris@49:         for(uword col=0; col < n_cols; ++col)  { tmp_vec[col] = P.at(0,col); }
Chris@49:         }
Chris@49:       else
Chris@49:         {
Chris@49:         arma_stop("op_median::median_vec(): expected a vector" );
Chris@49:         }
Chris@49:       }
Chris@49:     }
Chris@49:   
Chris@49:   return op_median::direct_median(tmp_vec);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename T1>
Chris@49: inline
Chris@49: typename T1::elem_type
Chris@49: op_median::median_vec
Chris@49:   (
Chris@49:   const T1& X,
Chris@49:   const typename arma_cx_only<typename T1::elem_type>::result* junk
Chris@49:   )
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   arma_ignore(junk);
Chris@49:   
Chris@49:   typedef typename T1::elem_type eT;
Chris@49:   typedef typename T1::pod_type   T;
Chris@49:   
Chris@49:   const Proxy<T1> P(X);
Chris@49:   
Chris@49:   const uword n_elem = P.get_n_elem();
Chris@49:   
Chris@49:   arma_debug_check( (n_elem == 0), "median(): given object has no elements" );
Chris@49:   
Chris@49:   std::vector< arma_cx_median_packet<T> > tmp_vec(n_elem);
Chris@49:   
Chris@49:   if(Proxy<T1>::prefer_at_accessor == false)
Chris@49:     {
Chris@49:     typedef typename Proxy<T1>::ea_type ea_type;
Chris@49:     
Chris@49:     ea_type A = P.get_ea();
Chris@49:     
Chris@49:     for(uword i=0; i<n_elem; ++i)
Chris@49:       {
Chris@49:       tmp_vec[i].val   = std::abs( A[i] );
Chris@49:       tmp_vec[i].index = i;
Chris@49:       }
Chris@49:     
Chris@49:     uword index1;
Chris@49:     uword index2;
Chris@49:     op_median::direct_cx_median_index(index1, index2, tmp_vec);
Chris@49:     
Chris@49:     return op_mean::robust_mean( A[index1], A[index2] );
Chris@49:     }
Chris@49:   else
Chris@49:     {
Chris@49:     const uword n_rows = P.get_n_rows();
Chris@49:     const uword n_cols = P.get_n_cols();
Chris@49:     
Chris@49:     if(n_cols == 1)
Chris@49:       {
Chris@49:       for(uword row=0; row < n_rows; ++row)
Chris@49:         {
Chris@49:         tmp_vec[row].val   = std::abs( P.at(row,0) );
Chris@49:         tmp_vec[row].index = row;
Chris@49:         }
Chris@49:       
Chris@49:       uword index1;
Chris@49:       uword index2;
Chris@49:       op_median::direct_cx_median_index(index1, index2, tmp_vec);
Chris@49:       
Chris@49:       return op_mean::robust_mean( P.at(index1,0), P.at(index2,0) );
Chris@49:       }
Chris@49:     else
Chris@49:     if(n_rows == 1)
Chris@49:       {
Chris@49:       for(uword col=0; col < n_cols; ++col)
Chris@49:         {
Chris@49:         tmp_vec[col].val   = std::abs( P.at(0,col) );
Chris@49:         tmp_vec[col].index = col;
Chris@49:         }
Chris@49:       
Chris@49:       uword index1;
Chris@49:       uword index2;
Chris@49:       op_median::direct_cx_median_index(index1, index2, tmp_vec);
Chris@49:       
Chris@49:       return op_mean::robust_mean( P.at(0,index1), P.at(0,index2) );
Chris@49:       }
Chris@49:     else
Chris@49:       {
Chris@49:       arma_stop("op_median::median_vec(): expected a vector" );
Chris@49:       
Chris@49:       return eT(0);
Chris@49:       }
Chris@49:     }
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: //! find the median value of a std::vector (contents is modified)
Chris@49: template<typename eT>
Chris@49: inline 
Chris@49: eT
Chris@49: op_median::direct_median(std::vector<eT>& X)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   const uword n_elem = uword(X.size());
Chris@49:   const uword half   = n_elem/2;
Chris@49:   
Chris@49:   std::nth_element(X.begin(), X.begin() + half, X.end());
Chris@49:   
Chris@49:   if((n_elem % 2) == 0)
Chris@49:     {
Chris@49:     return op_mean::robust_mean(*(std::max_element(X.begin(), X.begin() + half)), X[half]);
Chris@49:     }
Chris@49:   else
Chris@49:     {
Chris@49:     return X[half];
Chris@49:     }
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename T>
Chris@49: inline 
Chris@49: void
Chris@49: op_median::direct_cx_median_index
Chris@49:   (
Chris@49:   uword& out_index1, 
Chris@49:   uword& out_index2, 
Chris@49:   std::vector< arma_cx_median_packet<T> >& X
Chris@49:   )
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   const uword n_elem = uword(X.size());
Chris@49:   const uword half   = n_elem/2;
Chris@49:   
Chris@49:   std::nth_element(X.begin(), X.begin() + half, X.end());
Chris@49:   
Chris@49:   if((n_elem % 2) == 0)
Chris@49:     {
Chris@49:     out_index1 = std::max_element(X.begin(), X.begin() + half)->index;
Chris@49:     out_index2 = X[half].index;
Chris@49:     }
Chris@49:   else
Chris@49:     {
Chris@49:     out_index1 = X[half].index;
Chris@49:     out_index2 = out_index1;
Chris@49:     }
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: //! @}
Chris@49: