Chris@49: // Copyright (C) 2012-2013 NICTA (www.nicta.com.au)
Chris@49: // Copyright (C) 2012-2013 Conrad Sanderson
Chris@49: // Copyright (C) 2012 Boris Sabanin
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: 
Chris@49: template<typename T1, typename T2>
Chris@49: inline
Chris@49: void
Chris@49: glue_histc::apply(Mat<uword>& out, const mtGlue<uword,T1,T2,glue_histc>& 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;
Chris@49:   
Chris@49:   const unwrap_check_mixed<T1> tmp1(in.A, out);
Chris@49:   const unwrap_check_mixed<T2> tmp2(in.B, out);
Chris@49:   
Chris@49:   const Mat<eT>& X = tmp1.M;
Chris@49:   const Mat<eT>& E = tmp2.M;
Chris@49:   
Chris@49:   arma_debug_check
Chris@49:     (
Chris@49:     ((E.is_vec() == false) && (E.is_empty() == false)),
Chris@49:     "histc(): parameter 'edges' must be a vector"
Chris@49:     );
Chris@49:   
Chris@49:   arma_debug_check
Chris@49:     (
Chris@49:     (dim > 1),
Chris@49:     "histc(): parameter 'dim' must be 0 or 1"
Chris@49:     );
Chris@49:   
Chris@49:   const uword X_n_elem = X.n_elem;
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 E_n_elem = E.n_elem;
Chris@49:   
Chris@49:   if( E_n_elem == 0 )
Chris@49:     {
Chris@49:     out.reset();
Chris@49:     return;
Chris@49:     }
Chris@49:   
Chris@49:   
Chris@49:   // for vectors we are currently ignoring the "dim" parameter
Chris@49:   
Chris@49:   uword out_n_rows = 0;
Chris@49:   uword out_n_cols = 0;
Chris@49:   
Chris@49:   if(X.is_vec())
Chris@49:     {
Chris@49:     if(X.is_rowvec())
Chris@49:       {
Chris@49:       out_n_rows = 1;
Chris@49:       out_n_cols = E_n_elem;
Chris@49:       }
Chris@49:     else
Chris@49:     if(X.is_colvec())
Chris@49:       {
Chris@49:       out_n_rows = E_n_elem;
Chris@49:       out_n_cols = 1;
Chris@49:       }
Chris@49:     }
Chris@49:   else
Chris@49:     {
Chris@49:     if(dim == 0)
Chris@49:       {
Chris@49:       out_n_rows = E_n_elem;
Chris@49:       out_n_cols = X_n_cols;
Chris@49:       }
Chris@49:     else
Chris@49:     if(dim == 1)
Chris@49:       {
Chris@49:       out_n_rows = X_n_rows;
Chris@49:       out_n_cols = E_n_elem;
Chris@49:       }
Chris@49:     }
Chris@49:   
Chris@49:   out.zeros(out_n_rows, out_n_cols);
Chris@49:   
Chris@49:   const eT* E_mem = E.memptr();
Chris@49: 
Chris@49:   if(X.is_vec() == true)
Chris@49:     {
Chris@49:           uword* out_mem = out.memptr();
Chris@49:     const eT*    X_mem   = X.memptr();
Chris@49:     
Chris@49:     for(uword j=0; j<X_n_elem; ++j)
Chris@49:       {
Chris@49:       const eT val = X_mem[j];
Chris@49:       
Chris@49:       for(uword i=0; i<E_n_elem-1; ++i)
Chris@49:         {
Chris@49:         if( (E_mem[i] <= val) && (val < E_mem[i+1]) )
Chris@49:           {
Chris@49:           out_mem[i]++;
Chris@49:           break;
Chris@49:           }
Chris@49:         else
Chris@49:         if(val == E_mem[E_n_elem-1])
Chris@49:           {
Chris@49:           // in general, the above == operation doesn't make sense for floating point values (due to precision issues),
Chris@49:           // but is included for compatibility with Matlab and Octave.
Chris@49:           // Matlab folks must have been smoking something strong.
Chris@49:           out_mem[E_n_elem-1]++;
Chris@49:           break;
Chris@49:           }
Chris@49:         }
Chris@49:       }
Chris@49:     }
Chris@49:   else
Chris@49:   if(dim == 0)
Chris@49:     {
Chris@49:     for(uword col=0; col<X_n_cols; ++col)
Chris@49:       {
Chris@49:             uword* out_coldata = out.colptr(col);
Chris@49:       const eT*    X_coldata   = X.colptr(col);
Chris@49:       
Chris@49:       for(uword row=0; row<X_n_rows; ++row)
Chris@49:         {
Chris@49:         const eT val = X_coldata[row];
Chris@49:         
Chris@49:         for(uword i=0; i<E_n_elem-1; ++i)
Chris@49:           {
Chris@49:           if( (E_mem[i] <= val) && (val < E_mem[i+1]) )
Chris@49:             {
Chris@49:             out_coldata[i]++;
Chris@49:             break;
Chris@49:             }
Chris@49:           else
Chris@49:           if(val == E_mem[E_n_elem-1])
Chris@49:             {
Chris@49:             out_coldata[E_n_elem-1]++;
Chris@49:             break;
Chris@49:             }
Chris@49:           }
Chris@49:         }
Chris@49:       }
Chris@49:     }
Chris@49:   else
Chris@49:   if(dim == 1)
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:         const eT val = X.at(row,col);
Chris@49:         
Chris@49:         for(uword i=0; i<E_n_elem-1; ++i)
Chris@49:           {
Chris@49:           if( (E_mem[i] <= val) && (val < E_mem[i+1]) )
Chris@49:             {
Chris@49:             out.at(row,i)++;
Chris@49:             break;
Chris@49:             }
Chris@49:           else
Chris@49:           if(val == E_mem[E_n_elem-1])
Chris@49:             {
Chris@49:             out.at(row,E_n_elem-1)++;
Chris@49:             break;
Chris@49:             }
Chris@49:           }
Chris@49:         }
Chris@49:       }
Chris@49:     }
Chris@49:   }