C++ Mat::at方法代码示例

inline 
eT
op_min::direct_min(const Mat<eT>& X, const uword row)
  {
  arma_extra_debug_sigprint();
  
  const uword X_n_cols = X.n_cols;
  
  eT min_val = priv::most_pos<eT>();
  
  uword i,j;
  for(i=0, j=1; j < X_n_cols; i+=2, j+=2)
    {
    const eT tmp_i = X.at(row,i);
    const eT tmp_j = X.at(row,j);
    
    if(tmp_i < min_val) { min_val = tmp_i; }
    if(tmp_j < min_val) { min_val = tmp_j; }
    }
  
  if(i < X_n_cols)
    {
    const eT tmp_i = X.at(row,i);
    
    if(tmp_i < min_val) { min_val = tmp_i; }
    }
  
  return min_val;
  }

arma_hot
inline
eT
op_dot::dot_and_copy_row(eT* out, const Mat<eT>& A, const uword row, const eT* B_mem, const uword N)
  {
  eT acc1 = eT(0);
  eT acc2 = eT(0);
  
  uword i,j;
  for(i=0, j=1; j < N; i+=2, j+=2)
    {
    const eT val_i = A.at(row, i);
    const eT val_j = A.at(row, j);
    
    out[i] = val_i;
    out[j] = val_j;
    
    acc1 += val_i * B_mem[i];
    acc2 += val_j * B_mem[j];
    }
  
  if(i < N)
    {
    const eT val_i = A.at(row, i);
    
    out[i] = val_i;
    
    acc1 += val_i * B_mem[i];
    }
  
  return acc1 + acc2;
  }

inline
bool
op_logmat::apply_direct(Mat< std::complex<typename T1::elem_type> >& out, const Op<T1,op_diagmat>& expr, const uword)
  {
  arma_extra_debug_sigprint();
  
  typedef typename T1::elem_type T;
  
  const diagmat_proxy<T1> P(expr.m);
  
  arma_debug_check( (P.n_rows != P.n_cols), "logmat(): given matrix must be square sized" );
  
  const uword N = P.n_rows;
  
  out.zeros(N,N);
  
  for(uword i=0; i<N; ++i)
    {
    const T val = P[i];
    
    if(val >= T(0))
      {
      out.at(i,i) = std::log(val);
      }
    else
      {
      out.at(i,i) = std::log( std::complex<T>(val) );
      }
    }
  
  return true;
  }

inline
eT
op_mean::direct_mean(const Mat<eT>& X, const uword row)
  {
  arma_extra_debug_sigprint();
  
  typedef typename get_pod_type<eT>::result T;
  
  const uword X_n_cols = X.n_cols;
  
  eT val = eT(0);
  
  uword i,j;
  for(i=0, j=1; j < X_n_cols; i+=2, j+=2)
    {
    val += X.at(row,i);
    val += X.at(row,j);
    }
  
  if(i < X_n_cols)
    {
    val += X.at(row,i);
    }
  
  const eT result = val / T(X_n_cols);
  
  return arma_isfinite(result) ? result : op_mean::direct_mean_robust(X, row);
  }

inline
void
op_diagmat2::apply(Mat<typename T1::elem_type>& out, const Proxy<T1>& P, const uword row_offset, const uword col_offset)
  {
  arma_extra_debug_sigprint();
  
  const uword n_rows = P.get_n_rows();
  const uword n_cols = P.get_n_cols();
  const uword n_elem = P.get_n_elem();
  
  if(n_elem == 0)  { out.reset(); return; }
  
  const bool P_is_vec = (T1::is_row) || (T1::is_col) || (n_rows == 1) || (n_cols == 1);
  
  if(P_is_vec)
    {
    const uword n_pad = (std::max)(row_offset, col_offset);
    
    out.zeros(n_elem + n_pad, n_elem + n_pad);
    
    if(Proxy<T1>::prefer_at_accessor == false)
      {
      typename Proxy<T1>::ea_type Pea = P.get_ea();
      
      for(uword i=0; i < n_elem; ++i)
        {
        out.at(row_offset + i, col_offset + i) = Pea[i];
        }
      }
    else
      {
      const unwrap<typename Proxy<T1>::stored_type> U(P.Q);
      
      const Proxy<typename unwrap<typename Proxy<T1>::stored_type>::stored_type> PP(U.M);
      
      op_diagmat2::apply(out, PP, row_offset, col_offset);
      }
    }
  else  // P represents a matrix 
    {
    arma_debug_check
      (
      ((row_offset > 0) && (row_offset >= n_rows)) || ((col_offset > 0) && (col_offset >= n_cols)),
      "diagmat(): requested diagonal out of bounds"
      );
    
    out.zeros(n_rows, n_cols);
    
    const uword N = (std::min)(n_rows - row_offset, n_cols - col_offset);
    
    for(uword i=0; i<N; ++i)
      {
      const uword row = i + row_offset;
      const uword col = i + col_offset;
      
      out.at(row,col) = P.at(row,col);
      }
    }
  }

inline
void
arma_ostream::print(std::ostream& o, const Mat<eT>& m, const bool modify)
  {
  arma_extra_debug_sigprint();
  
  const arma_ostream_state stream_state(o);
  
  const std::streamsize cell_width = modify ? arma_ostream::modify_stream(o, m.memptr(), m.n_elem) : o.width();
  
  const uword m_n_rows = m.n_rows;
  const uword m_n_cols = m.n_cols;
  
  if(m.is_empty() == false)
    {
    if(m_n_cols > 0)
      {
      if(cell_width > 0)
        {
        for(uword row=0; row < m_n_rows; ++row)
          {
          for(uword col=0; col < m_n_cols; ++col)
            {
            // the cell width appears to be reset after each element is printed,
            // hence we need to restore it
            o.width(cell_width);
            arma_ostream::print_elem(o, m.at(row,col), modify);
            }
        
          o << '\n';
          }
        }
      else
        {
        for(uword row=0; row < m_n_rows; ++row)
          {
          for(uword col=0; col < m_n_cols-1; ++col)
            {
            arma_ostream::print_elem(o, m.at(row,col), modify);
            o << ' ';
            }
        
          arma_ostream::print_elem(o, m.at(row, m_n_cols-1), modify);
          o << '\n';
          }
        }
      }
    }
  else
    {
    o << "[matrix size: " << m_n_rows << 'x' << m_n_cols << "]\n";
    }
  
  o.flush();
  stream_state.restore(o);
  }

arma_hot
inline
void
op_strans2::apply(Mat<eT>& out, const TA& A, const eT val)
  {
  arma_extra_debug_sigprint();
  
  if(&out != &A)
    {
    op_strans2::apply_noalias(out, A, val);
    }
  else
    {
    const uword n_rows = out.n_rows;
    const uword n_cols = out.n_cols;
    
    if(n_rows == n_cols)
      {
      arma_extra_debug_print("op_strans2::apply(): doing in-place transpose of a square matrix");
      
      const uword N = n_rows;
      
      // TODO: do multiplication while swapping
      
      for(uword k=0; k < N; ++k)
        {
        eT* colptr = out.colptr(k);
        
        uword i,j;
        
        for(i=(k+1), j=(k+2); j < N; i+=2, j+=2)
          {
          std::swap(out.at(k,i), colptr[i]);
          std::swap(out.at(k,j), colptr[j]);
          }
        
        if(i < N)
          {
          std::swap(out.at(k,i), colptr[i]);
          }
        }
      
      arrayops::inplace_mul( out.memptr(), val, out.n_elem );
      }
    else
      {
      Mat<eT> tmp;
      op_strans2::apply_noalias(tmp, A, val);
      
      out.steal_mem(tmp);
      }
    }
  }

inline
void
op_trimat::apply_htrans
  (
        Mat<eT>& out,
  const Mat<eT>& A,
  const bool     upper,
  const typename arma_cx_only<eT>::result* junk
  )
  {
  arma_extra_debug_sigprint();
  arma_ignore(junk);
  
  arma_debug_check( (A.is_square() == false), "trimatu()/trimatl(): given matrix must be square sized" );
  
  const uword N = A.n_rows;
  
  if(&out != &A)
    {
    out.copy_size(A);
    }
  
  if(upper)
    {
    // Upper triangular: but since we're transposing, we're taking the lower
    // triangular and putting it in the upper half.
    for(uword row = 0; row < N; ++row)
      {
      eT* out_colptr = out.colptr(row);
      
      for(uword col = 0; col <= row; ++col)
        {
        //out.at(col, row) = std::conj( A.at(row, col) );
        out_colptr[col] = std::conj( A.at(row, col) );
        }
      }
    }
  else
    {
    // Lower triangular: but since we're transposing, we're taking the upper
    // triangular and putting it in the lower half.
    for(uword row = 0; row < N; ++row)
      {
      for(uword col = row; col < N; ++col)
        {
        out.at(col, row) = std::conj( A.at(row, col) );
        }
      }
    }
  
  op_trimat::fill_zeros(out, upper);
  }

arma_hot
inline
void
podarray<eT>::copy_row(const Mat<eT>& A, const uword row)
  {
  const uword cols = A.n_cols;
  
  // note: this function assumes that the podarray has been set to the correct size beforehand
  eT* out = memptr();
  
  switch(cols)
    {
    default:
      {
      uword i,j;
      for(i=0, j=1; j < cols; i+=2, j+=2)
        {
        const eT tmp_i = A.at(row, i);
        const eT tmp_j = A.at(row, j);
        
        out[i] = tmp_i;
        out[j] = tmp_j;
        }
      
      if(i < cols)
        {
        out[i] = A.at(row, i);
        }
      }
      break;
    
    case 8:  out[7] = A.at(row, 7);
    // fallthrough
    case 7:  out[6] = A.at(row, 6);
    // fallthrough
    case 6:  out[5] = A.at(row, 5);
    // fallthrough
    case 5:  out[4] = A.at(row, 4);
    // fallthrough
    case 4:  out[3] = A.at(row, 3);
    // fallthrough
    case 3:  out[2] = A.at(row, 2);
    // fallthrough
    case 2:  out[1] = A.at(row, 1);
    // fallthrough
    case 1:  out[0] = A.at(row, 0);
    // fallthrough
    case 0:  ;
    // fallthrough
    }
  }

arma_hot
inline
void
eop_core<eop_type>::apply_unwrap(Mat<typename T1::elem_type>& out, const eOp<T1, eop_type>& x)
  {
  arma_extra_debug_sigprint();
  
  typedef typename T1::elem_type eT;
  
  const Proxy<T1>& P = x.P;
  
  out.set_size(P.n_rows, P.n_cols);
  
        eT* out_mem = out.memptr();
  const u32 n_elem  = P.n_elem;
    
  const unwrap<typename Proxy<T1>::stored_type> tmp(P.Q);
  const Mat<eT>& A = tmp.M;
  
  if(is_generator<eop_type>::value == true)
    {
    for(u32 i=0; i<n_elem; ++i)
      {
      out_mem[i] = eop_aux::generate<eT,eop_type>();
      }
    }
  else
    {
    if(is_same_type<eop_type, eop_ones_diag>::value == true)
      {
      for(u32 col=0; col<P.n_rows; ++col)
        {
        for(u32 row=0; row<col; ++row)          { out.at(row,col) = eT(0); }
          
        out.at(col,col) = eT(1);
        
        for(u32 row=col+1; row<P.n_rows; ++row) { out.at(row,col) = eT(0); }
        }
      }
    else
      {
      const eT* A_mem = A.memptr();
      
      for(u32 i=0; i<n_elem; ++i)
        {
        out_mem[i] = eop_core<eop_type>::process(x, A_mem[i]);
        }
      }
    }
  }

arma_hot
inline
void
eop_core<eop_type>::apply_proxy(Mat<typename T1::elem_type>& out, const eOp<T1, eop_type>& x)
  {
  arma_extra_debug_sigprint();
  
  // eop_type::apply_proxy() function is not allowed to unwrap things
  // (in order to get the input into a common format).
  // the proxy class is already providing objects with element access
  
  typedef typename T1::elem_type eT;
  
  const Proxy<T1>& P = x.P;
  
  out.set_size(P.n_rows, P.n_cols);
  
        eT* out_mem = out.memptr();
  const u32 n_elem  = P.n_elem;
    
  if(is_generator<eop_type>::value == true)
    {
    for(u32 i=0; i<n_elem; ++i)
      {
      out_mem[i] = eop_aux::generate<eT,eop_type>();
      }
    }
  else
    {
    if(is_same_type<eop_type, eop_ones_diag>::value == true)
      {
      for(u32 col=0; col<P.n_rows; ++col)
        {
        for(u32 row=0; row<col; ++row)          { out.at(row,col) = eT(0); }
          
        out.at(col,col) = eT(1);
        
        for(u32 row=col+1; row<P.n_rows; ++row) { out.at(row,col) = eT(0); }
        }
      }
    else
      {
      for(u32 i=0; i<n_elem; ++i)
        {
        out_mem[i] = eop_core<eop_type>::process(x, P[i]);
        }
      }
    }
  }

Mat<N>
Mat<N>::operator*(const Mat<N>& r) const
{
    Mat<N> ret;
    for (uint8_t x = 0; x < N; x++) {
        for (uint8_t y = 0; y < N; y++) {
            float sum = 0.0f;
            for (uint8_t i = 0; i < N; i++) {
                sum += at(i,y) * r.at(x,i);
            }
            ret.at(x,y) = sum;
        }
    }
    return ret;
}

inline
bool
op_sqrtmat_cx::apply_direct_noalias(Mat<typename T1::elem_type>& out, const diagmat_proxy<T1>& P)
  {
  arma_extra_debug_sigprint();
  
  typedef typename T1::elem_type eT;
  
  arma_debug_check( (P.n_rows != P.n_cols), "sqrtmat(): given matrix must be square sized" );
  
  const uword N = P.n_rows;
  
  out.zeros(N,N);
  
  const eT zero = eT(0);
  
  bool singular = false;
  
  for(uword i=0; i<N; ++i)
    {
    const eT val = P[i];
    
    singular = (singular || (val == zero));
    
    out.at(i,i) = std::sqrt(val);
    }
  
  return (singular) ? false : true;
  }

inline
bool
diskio::save_pgm_binary(const Mat<eT>& x, std::ostream& f)
  {
  arma_extra_debug_sigprint();
  
  f << "P5" << '\n';
  f << x.n_cols << ' ' << x.n_rows << '\n';
  f << 255 << '\n';
  
  const u32 n_elem = x.n_rows * x.n_cols;
  podarray<u8> tmp(n_elem);
  
  u32 i = 0;
  
  for(u32 row=0; row < x.n_rows; ++row)
    {
    for(u32 col=0; col < x.n_cols; ++col)
      {
      tmp[i] = u8( x.at(row,col) );  // TODO: add round() ?
      ++i;
      }
    }
  
  f.write(reinterpret_cast<const char*>(tmp.mem), n_elem);
  
  return f.good();
  }

inline
void
glue_times::apply_inplace(Mat<typename T1::elem_type>& out, const T1& X)
  {
  arma_extra_debug_sigprint();
  
  typedef typename T1::elem_type eT;
  
  const unwrap_check<T1> tmp(X, out);
  const Mat<eT>& B     = tmp.M;
  
  arma_debug_assert_mul_size(out, B, "matrix multiply");
  
  if(out.n_cols == B.n_cols)
    {
    podarray<eT> tmp(out.n_cols);
    eT* tmp_rowdata = tmp.memptr();
    
    for(u32 out_row=0; out_row < out.n_rows; ++out_row)
      {
      for(u32 out_col=0; out_col < out.n_cols; ++out_col)
        {
        tmp_rowdata[out_col] = out.at(out_row,out_col);
        }
      
      for(u32 B_col=0; B_col < B.n_cols; ++B_col)
        {
        const eT* B_coldata = B.colptr(B_col);
        
        eT val = eT(0);
        for(u32 i=0; i < B.n_rows; ++i)
          {
          val += tmp_rowdata[i] * B_coldata[i];
          }
        
        out.at(out_row,B_col) = val;
        }
      }
    
    }
  else
    {
    const Mat<eT> tmp(out);
    glue_times::apply(out, tmp, B, eT(1), false, false, false);
    }
  
  }