Chris@49: // Copyright (C) 2008-2012 NICTA (www.nicta.com.au)
Chris@49: // Copyright (C) 2008-2012 Conrad Sanderson
Chris@49: // Copyright (C) 2012 Ryan Curtin
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 fn_accu
Chris@49: //! @{
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename T1>
Chris@49: arma_hot
Chris@49: inline
Chris@49: typename T1::elem_type
Chris@49: accu_proxy_linear(const Proxy<T1>& P)
Chris@49:   {
Chris@49:   typedef typename T1::elem_type      eT;
Chris@49:   typedef typename Proxy<T1>::ea_type ea_type;
Chris@49:   
Chris@49:         ea_type A      = P.get_ea();
Chris@49:   const uword   n_elem = P.get_n_elem();
Chris@49:   
Chris@49:   eT val1 = eT(0);
Chris@49:   eT val2 = eT(0);
Chris@49:   
Chris@49:   uword i,j;
Chris@49:   for(i=0, j=1; j < n_elem; i+=2, j+=2)
Chris@49:     {
Chris@49:     val1 += A[i];
Chris@49:     val2 += A[j];
Chris@49:     }
Chris@49:   
Chris@49:   if(i < n_elem)
Chris@49:     {
Chris@49:     val1 += A[i];   // equivalent to: val1 += A[n_elem-1];
Chris@49:     }
Chris@49:   
Chris@49:   return (val1 + val2);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename T1>
Chris@49: arma_hot
Chris@49: inline
Chris@49: typename T1::elem_type
Chris@49: accu_proxy_at(const Proxy<T1>& P)
Chris@49:   {
Chris@49:   typedef typename T1::elem_type eT;
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:   eT val = eT(0);
Chris@49:     
Chris@49:   if(n_rows != 1)
Chris@49:     {
Chris@49:     for(uword col=0; col < n_cols; ++col)
Chris@49:     for(uword row=0; row < n_rows; ++row)
Chris@49:       {
Chris@49:       val += P.at(row,col);
Chris@49:       }
Chris@49:     }
Chris@49:   else
Chris@49:     {
Chris@49:     for(uword col=0; col < n_cols; ++col)
Chris@49:       {
Chris@49:       val += P.at(0,col);
Chris@49:       }
Chris@49:     }
Chris@49:   
Chris@49:   return val;
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: //! accumulate the elements of a matrix
Chris@49: template<typename T1>
Chris@49: arma_hot
Chris@49: inline
Chris@49: typename enable_if2< is_arma_type<T1>::value, typename T1::elem_type >::result
Chris@49: accu(const T1& X)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   const Proxy<T1> P(X);
Chris@49:   
Chris@49:   return (Proxy<T1>::prefer_at_accessor == false) ? accu_proxy_linear(P) : accu_proxy_at(P);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: //! explicit handling of Hamming norm (also known as zero norm)
Chris@49: template<typename T1>
Chris@49: inline
Chris@49: arma_warn_unused
Chris@49: uword
Chris@49: accu(const mtOp<uword,T1,op_rel_noteq>& X)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   typedef typename T1::elem_type eT;
Chris@49:   
Chris@49:   const eT val = X.aux;
Chris@49:   
Chris@49:   const Proxy<T1> P(X.m);
Chris@49:   
Chris@49:   uword n_nonzero = 0;
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:     const uword   n_elem = P.get_n_elem();
Chris@49:     
Chris@49:     for(uword i=0; i<n_elem; ++i)
Chris@49:       {
Chris@49:       if(A[i] != val) { ++n_nonzero; }
Chris@49:       }
Chris@49:     }
Chris@49:   else
Chris@49:     {
Chris@49:     const uword P_n_cols = P.get_n_cols();
Chris@49:     const uword P_n_rows = P.get_n_rows();
Chris@49:     
Chris@49:     if(P_n_rows == 1)
Chris@49:       {
Chris@49:       for(uword col=0; col < P_n_cols; ++col)
Chris@49:         {
Chris@49:         if(P.at(0,col) != val) { ++n_nonzero; }
Chris@49:         }
Chris@49:       }
Chris@49:     else
Chris@49:       {
Chris@49:       for(uword col=0; col < P_n_cols; ++col)
Chris@49:       for(uword row=0; row < P_n_rows; ++row)
Chris@49:         {
Chris@49:         if(P.at(row,col) != val) { ++n_nonzero; }
Chris@49:         }
Chris@49:       }
Chris@49:     }
Chris@49:   
Chris@49:   return n_nonzero;
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: //! accumulate the elements of a subview (submatrix)
Chris@49: template<typename eT>
Chris@49: arma_hot
Chris@49: arma_pure
Chris@49: arma_warn_unused
Chris@49: inline
Chris@49: eT
Chris@49: accu(const subview<eT>& X)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();  
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:   eT val = eT(0);
Chris@49:   
Chris@49:   if(X_n_rows == 1)
Chris@49:     {
Chris@49:     const Mat<eT>& A = X.m;
Chris@49:     
Chris@49:     const uword start_row = X.aux_row1;
Chris@49:     const uword start_col = X.aux_col1;
Chris@49:     
Chris@49:     const uword end_col_p1 = start_col + X_n_cols;
Chris@49:     
Chris@49:     uword i,j;
Chris@49:     for(i=start_col, j=start_col+1; j < end_col_p1; i+=2, j+=2)
Chris@49:       {
Chris@49:       val += A.at(start_row, i);
Chris@49:       val += A.at(start_row, j);
Chris@49:       }
Chris@49:     
Chris@49:     if(i < end_col_p1)
Chris@49:       {
Chris@49:       val += A.at(start_row, i);
Chris@49:       }
Chris@49:     }
Chris@49:   else
Chris@49:   if(X_n_cols == 1)
Chris@49:     {
Chris@49:     val = arrayops::accumulate( X.colptr(0), X_n_rows );
Chris@49:     }
Chris@49:   else
Chris@49:     {
Chris@49:     for(uword col=0; col < X_n_cols; ++col)
Chris@49:       {
Chris@49:       val += arrayops::accumulate( X.colptr(col), X_n_rows );
Chris@49:       }
Chris@49:     }
Chris@49:   
Chris@49:   return val;
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: arma_hot
Chris@49: arma_pure
Chris@49: arma_warn_unused
Chris@49: inline
Chris@49: eT
Chris@49: accu(const subview_col<eT>& X)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();  
Chris@49:   
Chris@49:   return arrayops::accumulate( X.colptr(0), X.n_rows );
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: //! accumulate the elements of a cube
Chris@49: template<typename T1>
Chris@49: arma_hot
Chris@49: arma_warn_unused
Chris@49: inline
Chris@49: typename T1::elem_type
Chris@49: accu(const BaseCube<typename T1::elem_type,T1>& X)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   typedef typename T1::elem_type          eT;
Chris@49:   typedef typename ProxyCube<T1>::ea_type ea_type;
Chris@49:   
Chris@49:   const ProxyCube<T1> A(X.get_ref());
Chris@49:   
Chris@49:   if(ProxyCube<T1>::prefer_at_accessor == false)
Chris@49:     {
Chris@49:           ea_type P      = A.get_ea();
Chris@49:     const uword   n_elem = A.get_n_elem();
Chris@49:     
Chris@49:     eT val1 = eT(0);
Chris@49:     eT val2 = eT(0);
Chris@49:     
Chris@49:     uword i,j;
Chris@49:     
Chris@49:     for(i=0, j=1; j<n_elem; i+=2, j+=2)
Chris@49:       {
Chris@49:       val1 += P[i];
Chris@49:       val2 += P[j];
Chris@49:       }
Chris@49:     
Chris@49:     if(i < n_elem)
Chris@49:       {
Chris@49:       val1 += P[i];
Chris@49:       }
Chris@49:     
Chris@49:     return val1 + val2;
Chris@49:     }
Chris@49:   else
Chris@49:     {
Chris@49:     const uword n_rows   = A.get_n_rows();
Chris@49:     const uword n_cols   = A.get_n_cols();
Chris@49:     const uword n_slices = A.get_n_slices();
Chris@49:     
Chris@49:     eT val = eT(0);
Chris@49:     
Chris@49:     for(uword slice=0; slice<n_slices; ++slice)
Chris@49:     for(uword col=0; col<n_cols; ++col)
Chris@49:     for(uword row=0; row<n_rows; ++row)
Chris@49:       {
Chris@49:       val += A.at(row,col,slice);
Chris@49:       }
Chris@49:     
Chris@49:     return val;
Chris@49:     }
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename T>
Chris@49: arma_inline
Chris@49: arma_warn_unused
Chris@49: const typename arma_scalar_only<T>::result &
Chris@49: accu(const T& x)
Chris@49:   {
Chris@49:   return x;
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: //! accumulate values in a sparse object
Chris@49: template<typename T1>
Chris@49: arma_hot
Chris@49: inline
Chris@49: arma_warn_unused
Chris@49: typename enable_if2<is_arma_sparse_type<T1>::value, typename T1::elem_type>::result
Chris@49: accu(const T1& x)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   typedef typename T1::elem_type eT;
Chris@49:   
Chris@49:   const SpProxy<T1> p(x);
Chris@49:   
Chris@49:   if(SpProxy<T1>::must_use_iterator == false)
Chris@49:     {
Chris@49:     // direct counting
Chris@49:     return arrayops::accumulate(p.get_values(), p.get_n_nonzero());
Chris@49:     }
Chris@49:   else
Chris@49:     {
Chris@49:     typename SpProxy<T1>::const_iterator_type it     = p.begin();
Chris@49:     typename SpProxy<T1>::const_iterator_type it_end = p.end();
Chris@49:     
Chris@49:     eT result = eT(0);
Chris@49:     
Chris@49:     while(it != it_end)
Chris@49:       {
Chris@49:       result += (*it);
Chris@49:       ++it;
Chris@49:       }
Chris@49:     
Chris@49:     return result;
Chris@49:     }
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: //! @}