Chris@49: // Copyright (C) 2012-2013 NICTA (www.nicta.com.au)
Chris@49: // Copyright (C) 2012-2013 Conrad Sanderson
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_hist::apply(Mat<uword>& out, const mtGlue<uword,T1,T2,glue_hist>& 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>& C = tmp2.M;
Chris@49:   
Chris@49:   
Chris@49:   arma_debug_check
Chris@49:     (
Chris@49:     ((C.is_vec() == false) && (C.is_empty() == false)),
Chris@49:     "hist(): parameter 'centers' must be a vector"
Chris@49:     );
Chris@49:   
Chris@49:   arma_debug_check
Chris@49:     (
Chris@49:     (dim > 1),
Chris@49:     "hist(): parameter 'dim' must be 0 or 1"
Chris@49:     );
Chris@49:   
Chris@49:   const uword X_n_rows = X.n_rows;
Chris@49:   const uword X_n_cols = X.n_cols;
Chris@49:   const uword X_n_elem = X.n_elem;
Chris@49:   
Chris@49:   const uword C_n_elem = C.n_elem;
Chris@49:   
Chris@49:   if( C_n_elem == 0 )
Chris@49:     {
Chris@49:     out.reset();
Chris@49:     return;
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 = C_n_elem;
Chris@49:       }
Chris@49:     else
Chris@49:     if(X.is_colvec())
Chris@49:       {
Chris@49:       out_n_rows = C_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 = C_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 = C_n_elem;
Chris@49:       }
Chris@49:     }
Chris@49:   
Chris@49:   out.zeros(out_n_rows, out_n_cols);
Chris@49:   
Chris@49:   
Chris@49:   const eT* C_mem    = C.memptr();
Chris@49:   const eT  center_0 = C_mem[0];
Chris@49:   
Chris@49:   if(X.is_vec())
Chris@49:     {
Chris@49:     const eT*    X_mem   = X.memptr();
Chris@49:           uword* out_mem = out.memptr();
Chris@49:     
Chris@49:     for(uword i=0; i < X_n_elem; ++i)
Chris@49:       {
Chris@49:       const eT val = X_mem[i];
Chris@49:       
Chris@49:       if(is_finite(val))
Chris@49:         {
Chris@49:         eT    opt_dist  = (val >= center_0) ? (val - center_0) : (center_0 - val);
Chris@49:         uword opt_index = 0;
Chris@49:         
Chris@49:         for(uword j=1; j < C_n_elem; ++j)
Chris@49:           {
Chris@49:           const eT center = C_mem[j];
Chris@49:           const eT dist   = (val >= center) ? (val - center) : (center - val);
Chris@49:           
Chris@49:           if(dist < opt_dist)
Chris@49:             {
Chris@49:             opt_dist  = dist;
Chris@49:             opt_index = j;
Chris@49:             }
Chris@49:           else
Chris@49:             {
Chris@49:             break;
Chris@49:             }
Chris@49:           }
Chris@49:         
Chris@49:         out_mem[opt_index]++;
Chris@49:         }
Chris@49:       else
Chris@49:         {
Chris@49:         // -inf
Chris@49:         if(val < eT(0)) { out_mem[0]++; }
Chris@49:         
Chris@49:         // +inf
Chris@49:         if(val > eT(0)) { out_mem[C_n_elem-1]++; }
Chris@49:         
Chris@49:         // ignore NaN
Chris@49:         }
Chris@49:       }
Chris@49:     }
Chris@49:   else
Chris@49:     {
Chris@49:     if(dim == 0)
Chris@49:       {
Chris@49:       for(uword col=0; col < X_n_cols; ++col)
Chris@49:         {
Chris@49:         const eT*    X_coldata   = X.colptr(col);
Chris@49:               uword* out_coldata = out.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:           if(arma_isfinite(val))
Chris@49:             {
Chris@49:             eT    opt_dist  = (center_0 >= val) ? (center_0 - val) : (val - center_0);
Chris@49:             uword opt_index = 0;
Chris@49:             
Chris@49:             for(uword j=1; j < C_n_elem; ++j)
Chris@49:               {
Chris@49:               const eT center = C_mem[j];
Chris@49:               const eT dist   = (center >= val) ? (center - val) : (val - center);
Chris@49:               
Chris@49:               if(dist < opt_dist)
Chris@49:                 {
Chris@49:                 opt_dist  = dist;
Chris@49:                 opt_index = j;
Chris@49:                 }
Chris@49:               else
Chris@49:                 {
Chris@49:                 break;
Chris@49:                 }
Chris@49:               }
Chris@49:             
Chris@49:             out_coldata[opt_index]++;
Chris@49:             }
Chris@49:           else
Chris@49:             {
Chris@49:             // -inf
Chris@49:             if(val < eT(0)) { out_coldata[0]++; }
Chris@49:             
Chris@49:             // +inf
Chris@49:             if(val > eT(0)) { out_coldata[C_n_elem-1]++; }
Chris@49:             
Chris@49:             // ignore NaN
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:           if(arma_isfinite(val))
Chris@49:             {
Chris@49:             eT    opt_dist  = (center_0 >= val) ? (center_0 - val) : (val - center_0);
Chris@49:             uword opt_index = 0;
Chris@49:             
Chris@49:             for(uword j=1; j < C_n_elem; ++j)
Chris@49:               {
Chris@49:               const eT center = C_mem[j];
Chris@49:               const eT dist   = (center >= val) ? (center - val) : (val - center);
Chris@49:               
Chris@49:               if(dist < opt_dist)
Chris@49:                 {
Chris@49:                 opt_dist  = dist;
Chris@49:                 opt_index = j;
Chris@49:                 }
Chris@49:               else
Chris@49:                 {
Chris@49:                 break;
Chris@49:                 }
Chris@49:               }
Chris@49:             
Chris@49:             out.at(row,opt_index)++;
Chris@49:             }
Chris@49:           else
Chris@49:             {
Chris@49:             // -inf
Chris@49:             if(val < eT(0)) { out.at(row,0)++; }
Chris@49:             
Chris@49:             // +inf
Chris@49:             if(val > eT(0)) { out.at(row,C_n_elem-1)++; }
Chris@49:             
Chris@49:             // ignore NaN
Chris@49:             }
Chris@49:           }
Chris@49:         }
Chris@49:       }
Chris@49:     }
Chris@49:   }