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

inline
void
subview_each2<parent,mode,TB>::check_indices(const Mat<uword>& indices) const
  {
  if(mode == 0)
    {
    arma_debug_check( ((indices.is_vec() == false) && (indices.is_empty() == false)), "each_col(): list of indices must be a vector" );
    }
  else
    {
    arma_debug_check( ((indices.is_vec() == false) && (indices.is_empty() == false)), "each_row(): list of indices must be a vector" );
    }
  }

inline
void
glue_cov::direct_cov(Mat< std::complex<T> >& out, const Mat< std::complex<T> >& A, const Mat< std::complex<T> >& B, const u32 norm_type)
  {
  arma_extra_debug_sigprint();

  typedef typename std::complex<T> eT;

  if(A.is_vec() && B.is_vec())
    { 
    arma_debug_check( (A.n_elem != B.n_elem), "cov(): the number of elements in A and B must match" );

    const eT* A_ptr = A.memptr();
    const eT* B_ptr = B.memptr();        

    eT A_acc   = eT(0);
    eT B_acc   = eT(0);
    eT out_acc = eT(0);

    const u32 N = A.n_elem;

    for(u32 i=0; i<N; ++i)
      {
      const eT A_tmp = A_ptr[i];
      const eT B_tmp = B_ptr[i];

      A_acc += A_tmp;
      B_acc += B_tmp;

      out_acc += std::conj(A_tmp) * B_tmp;
      }

    out_acc -= (std::conj(A_acc) * B_acc)/eT(N);

    const eT norm_val = (norm_type == 0) ? ( (N > 1) ? eT(N-1) : eT(1) ) : eT(N);

    out.set_size(1,1);
    out[0] = out_acc/norm_val;
    }
  else
    {
    arma_debug_assert_same_size(A, B, "cov()");
  
    const u32 N = A.n_rows;
    const eT norm_val = (norm_type == 0) ? ( (N > 1) ? eT(N-1) : eT(1) ) : eT(N);
    
    out = trans(A) * B;                     // out = strans(conj(A)) * B;
    out -= (trans(sum(A)) * sum(B))/eT(N);  // out -= (strans(conj(sum(A))) * sum(B))/eT(N);
    out /= norm_val;
    }
  }

inline
void
op_cov::direct_cov(Mat< std::complex<T> >& out, const Mat< std::complex<T> >& A, const u32 norm_type)
  {
  arma_extra_debug_sigprint();
  
  typedef typename std::complex<T> eT;
  
  if(A.is_vec())
    {
    if(A.n_rows == 1)
      {
      const Mat<T> tmp_mat = var(trans(A), norm_type);
      out.set_size(1,1);
      out[0] = tmp_mat[0];
      }
    else
      {
      const Mat<T> tmp_mat = var(A, norm_type);
      out.set_size(1,1);
      out[0] = tmp_mat[0];
      }
    }
  else
    {
    const u32 N = A.n_rows;
    const eT norm_val = (norm_type == 0) ? ( (N > 1) ? eT(N-1) : eT(1) ) : eT(N);
    
    const Row<eT> acc = sum(A);
    
    out = trans(A) * A;               // out = strans(conj(A)) * A;
    out -= (trans(acc) * acc)/eT(N);  // out -= (strans(conj(acc)) * acc)/eT(N);
    out /= norm_val;
    }
  }

inline
void
op_cor::direct_cor(Mat< std::complex<T> >& out, const Mat< std::complex<T> >& A, const uword norm_type)
  {
  arma_extra_debug_sigprint();

  typedef typename std::complex<T> eT;

  if(A.is_empty())
    {
    out.reset();
    return;
    }
  
  if(A.is_vec())
    {
    out.set_size(1,1);
    out[0] = eT(1);
    }
  else
    {
    const uword N = A.n_rows;
    const eT norm_val = (norm_type == 0) ? ( (N > 1) ? eT(N-1) : eT(1) ) : eT(N);

    const Row<eT> acc = sum(A);
    const Row<T>  sd  = stddev(A);

    out = trans(A) * A;               // out = strans(conj(A)) * A;
    out -= (trans(acc) * acc)/eT(N);  // out -= (strans(conj(acc)) * acc)/eT(N);
    out /= norm_val;

    //out = out / (trans(sd) * sd);
    out /= conv_to< Mat<eT> >::from(trans(sd) * sd);
    }
  }

inline
void
glue_trapz::apply_noalias(Mat<eT>& out, const Mat<eT>& X, const Mat<eT>& Y, const uword dim)
  {
  arma_extra_debug_sigprint();
  
  arma_debug_check( (dim > 1), "trapz(): argument 'dim' must be 0 or 1" );
  
  arma_debug_check( ((X.is_vec() == false) && (X.is_empty() == false)), "trapz(): argument 'X' must be a vector" );
  
  const uword N = X.n_elem;
  
  if(dim == 0)
    {
    arma_debug_check( (N != Y.n_rows), "trapz(): length of X must equal the number of rows in Y when dim=0" );
    }
  else
  if(dim == 1)
    {
    arma_debug_check( (N != Y.n_cols), "trapz(): length of X must equal the number of columns in Y when dim=1" );
    }
  
  if(N <= 1)
    {
         if(dim == 0)  { out.zeros(1, Y.n_cols); }
    else if(dim == 1)  { out.zeros(Y.n_rows, 1); }
    
    return;
    }
  
  const Col<eT> vec_X( const_cast<eT*>(X.memptr()), X.n_elem, false, true );
  
  const Col<eT> diff_X = diff(vec_X);
  
  if(dim == 0)
    {
    const Row<eT> diff_X_t( const_cast<eT*>(diff_X.memptr()), diff_X.n_elem, false, true );
    
    out = diff_X_t * (0.5 * (Y.rows(0, N-2) + Y.rows(1, N-1)));
    }
  else
  if(dim == 1)
    {
    out = (0.5 * (Y.cols(0, N-2) + Y.cols(1, N-1))) * diff_X;
    }
  }

	Mat<elem_type> diff(const Mat<elem_type>& X, size_type n = 1, size_type dim = 0)
	{
		static_assert(ARMA_VERSION_MAJOR <= 5 && ARMA_VERSION_MINOR < 400, "This function is deprecated. Use arma::diff instead.");

		assert(n > 0);
		Mat<elem_type> y;

		if (X.empty() || dim > 1 || X.n_elem == 1) return y;

		if (n > 1) return diff(diff(X, n - 1, dim));

		if (X.is_vec()) {
			y.resize(std::max(X.n_rows - 1, (size_type)1), std::max(X.n_cols - 1, (size_type)1));
#ifdef _MSC_VER
			std::adjacent_difference(X.begin(), X.end(),
				stdext::unchecked_array_iterator<elem_type *>(y.begin()));
#else
            std::adjacent_difference(X.begin(), X.end(), y.begin());
#endif
			return y;
		}
		
		switch (dim) {
		case 0: // row difference
			y.resize(X.n_rows - 1, X.n_cols);
			for (size_type c = 0 ; c < X.n_cols ; c++)
#ifdef _MSC_VER
				std::adjacent_difference(X.begin_col(c), X.end_col(c),
				stdext::unchecked_array_iterator<elem_type *>(y.begin_col(c)));
#else
                std::adjacent_difference(X.begin_col(c), X.end_col(c), y.begin_col(c));
#endif
			break;
		case 1: // column difference
			y.resize(X.n_rows, X.n_cols - 1);
			for (size_type c = 0 ; c < X.n_cols - 1; c++)
				y.col(c) = X.col(c + 1) - X.col(c);
			break;
		default:
			throw std::invalid_argument("dim should be 0 or 1"); //
		}

		return y;
	}

inline
void
subview_cube_each2<eT,TB>::check_indices(const Mat<uword>& indices) const
  {
  arma_debug_check( ((indices.is_vec() == false) && (indices.is_empty() == false)), "each_slice(): list of indices must be a vector" );
  }

inline
void
interp1_helper(const Mat<eT>& X, const Mat<eT>& Y, const Mat<eT>& XI, Mat<eT>& YI, const uword sig, const eT extrap_val)
  {
  arma_extra_debug_sigprint();
  
  arma_debug_check( ((X.is_vec() == false) || (Y.is_vec() == false) || (XI.is_vec() == false)), "interp1(): currently only vectors are supported" );
  
  arma_debug_check( (X.n_elem != Y.n_elem), "interp1(): X and Y must have the same number of elements" );
  
  arma_debug_check( (X.n_elem < 2), "interp1(): X must have at least two unique elements" );
  
  // sig = 10: nearest neighbour
  // sig = 11: nearest neighbour, assume monotonic increase in X and XI
  // 
  // sig = 20: linear
  // sig = 21: linear, assume monotonic increase in X and XI
  
  if(sig == 11)  { interp1_helper_nearest(X, Y, XI, YI, extrap_val); return; }
  if(sig == 21)  { interp1_helper_linear (X, Y, XI, YI, extrap_val); return; }
  
  uvec X_indices;
  
  try { X_indices = find_unique(X,false); } catch(...) { }
  
  // NOTE: find_unique(X,false) provides indices of elements sorted in ascending order
  // NOTE: find_unique(X,false) will reset X_indices if X has NaN
  
  const uword N_subset = X_indices.n_elem;
  
  arma_debug_check( (N_subset < 2), "interp1(): X must have at least two unique elements" );
  
  Mat<eT> X_sanitised(N_subset,1);
  Mat<eT> Y_sanitised(N_subset,1);
  
  eT* X_sanitised_mem = X_sanitised.memptr();
  eT* Y_sanitised_mem = Y_sanitised.memptr();
  
  const eT* X_mem = X.memptr();
  const eT* Y_mem = Y.memptr();
  
  const uword* X_indices_mem = X_indices.memptr();
  
  for(uword i=0; i<N_subset; ++i)
    {
    const uword j = X_indices_mem[i];
    
    X_sanitised_mem[i] = X_mem[j];
    Y_sanitised_mem[i] = Y_mem[j];
    }
  
  
  Mat<eT> XI_tmp;
  uvec    XI_indices;
  
  const bool XI_is_sorted = XI.is_sorted();
  
  if(XI_is_sorted == false)
    {
    XI_indices = sort_index(XI);
    
    const uword N = XI.n_elem;
    
    XI_tmp.copy_size(XI);
    
    const uword* XI_indices_mem = XI_indices.memptr();
    
    const eT* XI_mem     = XI.memptr();
          eT* XI_tmp_mem = XI_tmp.memptr();
    
    for(uword i=0; i<N; ++i)
      {
      XI_tmp_mem[i] = XI_mem[ XI_indices_mem[i] ];
      }
    }
  
  const Mat<eT>& XI_sorted = (XI_is_sorted) ? XI : XI_tmp;
  
  
       if(sig == 10)  { interp1_helper_nearest(X_sanitised, Y_sanitised, XI_sorted, YI, extrap_val); }
  else if(sig == 20)  { interp1_helper_linear (X_sanitised, Y_sanitised, XI_sorted, YI, extrap_val); }
  
  
  if( (XI_is_sorted == false) && (YI.n_elem > 0) )
    {
    Mat<eT> YI_unsorted;
    
    YI_unsorted.copy_size(YI);
    
    const eT* YI_mem          = YI.memptr();
          eT* YI_unsorted_mem = YI_unsorted.memptr();
    
    const uword  N              = XI_sorted.n_elem;
    const uword* XI_indices_mem = XI_indices.memptr();
    
    for(uword i=0; i<N; ++i)
      {
      YI_unsorted_mem[ XI_indices_mem[i] ] = YI_mem[i];
      }
    
//.........这里部分代码省略.........

inline
void
running_stat_vec_aux::update_stats(running_stat_vec< std::complex<T> >& x, const Mat< std::complex<T> >& sample)
  {
  arma_extra_debug_sigprint();
  
  typedef typename std::complex<T> eT;
  
  const T N = x.counter.value();
  
  if(N > T(0))
    {
    arma_debug_assert_same_size(x.r_mean, sample, "running_stat_vec(): dimensionality mismatch");
    
    const u32 n_elem           = sample.n_elem;
    const eT* sample_mem       = sample.memptr();
          eT* r_mean_mem       = x.r_mean.memptr();
           T* r_var_mem        = x.r_var.memptr();
          eT* min_val_mem      = x.min_val.memptr();
          eT* max_val_mem      = x.max_val.memptr();
           T* min_val_norm_mem = x.min_val_norm.memptr();
           T* max_val_norm_mem = x.max_val_norm.memptr();
    
    const T  N_plus_1   = x.counter.value_plus_1();
    const T  N_minus_1  = x.counter.value_minus_1();
    
    if(x.calc_cov == true)
      {
      Mat<eT>& tmp1 = x.tmp1;
      Mat<eT>& tmp2 = x.tmp2;
      
      tmp1 = sample - x.r_mean;
      
      if(sample.n_cols == 1)
        {
        tmp2 = conj(tmp1)*trans(tmp1);
        }
      else
        {
        tmp2 = trans(conj(tmp1))*tmp1;
        }
      
      x.r_cov *= (N_minus_1/N);
      x.r_cov += tmp2 / N_plus_1;
      }
    
    
    for(u32 i=0; i<n_elem; ++i)
      {
      const eT& val      = sample_mem[i];
      const  T  val_norm = std::norm(val);
      
      if(val_norm < min_val_norm_mem[i])
        {
        min_val_norm_mem[i] = val_norm;
        min_val_mem[i]      = val;
        }
      
      if(val_norm > max_val_norm_mem[i])
        {
        max_val_norm_mem[i] = val_norm;
        max_val_mem[i]      = val;
        }
      
      const eT& r_mean_val = r_mean_mem[i];
      
      r_var_mem[i] = N_minus_1/N * r_var_mem[i] + std::norm(val - r_mean_val)/N_plus_1;
      
      r_mean_mem[i] = r_mean_val + (val - r_mean_val)/N_plus_1;
      }
    
    }
  else
    {
    arma_debug_check( (sample.is_vec() == false), "running_stat_vec(): given sample is not a vector");
    
    x.r_mean.set_size(sample.n_rows, sample.n_cols);
    
    x.r_var.zeros(sample.n_rows, sample.n_cols);
    
    if(x.calc_cov == true)
      {
      x.r_cov.zeros(sample.n_elem, sample.n_elem);
      }
    
    x.min_val.set_size(sample.n_rows, sample.n_cols);
    x.max_val.set_size(sample.n_rows, sample.n_cols);
    
    x.min_val_norm.set_size(sample.n_rows, sample.n_cols);
    x.max_val_norm.set_size(sample.n_rows, sample.n_cols);
    
    
    const u32 n_elem           = sample.n_elem;
    const eT* sample_mem       = sample.memptr();
          eT* r_mean_mem       = x.r_mean.memptr();
          eT* min_val_mem      = x.min_val.memptr();
          eT* max_val_mem      = x.max_val.memptr();
           T* min_val_norm_mem = x.min_val_norm.memptr();
           T* max_val_norm_mem = x.max_val_norm.memptr();
    
//.........这里部分代码省略.........

inline
void
glue_hist::apply_noalias(Mat<uword>& out, const Mat<eT>& X, const Mat<eT>& C, const uword dim)
  {
  arma_extra_debug_sigprint();
  
  arma_debug_check( ((C.is_vec() == false) && (C.is_empty() == false)), "hist(): parameter 'centers' must be a vector" );
  
  const uword X_n_rows = X.n_rows;
  const uword X_n_cols = X.n_cols;
  
  const uword C_n_elem = C.n_elem;
  
  if( C_n_elem == 0 )  { out.reset(); return; }
  
  arma_debug_check
    (
    ((Col<eT>(const_cast<eT*>(C.memptr()), C_n_elem, false, false)).is_sorted("strictascend") == false),
    "hist(): given 'centers' vector does not contain monotonically increasing values"
    );
  
  const eT* C_mem    = C.memptr();
  const eT  center_0 = C_mem[0];
  
  if(dim == 0)
    {
    out.zeros(C_n_elem, X_n_cols);
    
    for(uword col=0; col < X_n_cols; ++col)
      {
      const eT*    X_coldata   = X.colptr(col);
            uword* out_coldata = out.colptr(col);
      
      for(uword row=0; row < X_n_rows; ++row)
        {
        const eT val = X_coldata[row];
        
        if(arma_isfinite(val))
          {
          eT    opt_dist  = (center_0 >= val) ? (center_0 - val) : (val - center_0);
          uword opt_index = 0;
          
          for(uword j=1; j < C_n_elem; ++j)
            {
            const eT center = C_mem[j];
            const eT dist   = (center >= val) ? (center - val) : (val - center);
            
            if(dist < opt_dist)
              {
              opt_dist  = dist;
              opt_index = j;
              }
            else
              {
              break;
              }
            }
          
          out_coldata[opt_index]++;
          }
        else
          {
          // -inf
          if(val < eT(0)) { out_coldata[0]++; }
          
          // +inf
          if(val > eT(0)) { out_coldata[C_n_elem-1]++; }
          
          // ignore NaN
          }
        }
      }
    }
  else
  if(dim == 1)
    {
    out.zeros(X_n_rows, C_n_elem);
    
    if(X_n_rows == 1)
      {
      const uword  X_n_elem = X.n_elem;
      const eT*    X_mem    = X.memptr();
            uword* out_mem  = out.memptr();
      
      for(uword i=0; i < X_n_elem; ++i)
        {
        const eT val = X_mem[i];
        
        if(is_finite(val))
          {
          eT    opt_dist  = (val >= center_0) ? (val - center_0) : (center_0 - val);
          uword opt_index = 0;
          
          for(uword j=1; j < C_n_elem; ++j)
            {
            const eT center = C_mem[j];
            const eT dist   = (val >= center) ? (val - center) : (center - val);
            
            if(dist < opt_dist)
              {
//.........这里部分代码省略.........

inline
void
running_stat_vec_aux::update_stats
  (
  running_stat_vec<obj_type>& x,
  const                  Mat<typename running_stat_vec<obj_type>::eT>& sample,
  const typename arma_not_cx<typename running_stat_vec<obj_type>::eT>::result* junk
  )
  {
  arma_extra_debug_sigprint();
  arma_ignore(junk);
  
  typedef typename running_stat_vec<obj_type>::eT eT;
  typedef typename running_stat_vec<obj_type>::T   T;
  
  const T N = x.counter.value();
  
  if(N > T(0))
    {
    arma_debug_assert_same_size(x.r_mean, sample, "running_stat_vec(): dimensionality mismatch");
    
    const uword n_elem      = sample.n_elem;
    const eT*   sample_mem  = sample.memptr();
          eT*   r_mean_mem  = x.r_mean.memptr();
           T*   r_var_mem   = x.r_var.memptr();
          eT*   min_val_mem = x.min_val.memptr();
          eT*   max_val_mem = x.max_val.memptr();
    
    const T  N_plus_1   = x.counter.value_plus_1();
    const T  N_minus_1  = x.counter.value_minus_1();
    
    if(x.calc_cov == true)
      {
      Mat<eT>& tmp1 = x.tmp1;
      Mat<eT>& tmp2 = x.tmp2;
      
      tmp1 = sample - x.r_mean;
      
      if(sample.n_cols == 1)
        {
        tmp2 = tmp1*trans(tmp1);
        }
      else
        {
        tmp2 = trans(tmp1)*tmp1;
        }
      
      x.r_cov *= (N_minus_1/N);
      x.r_cov += tmp2 / N_plus_1;
      }
    
    
    for(uword i=0; i<n_elem; ++i)
      {
      const eT val = sample_mem[i];
      
      if(val < min_val_mem[i])
        {
        min_val_mem[i] = val;
        }
      
      if(val > max_val_mem[i])
        {
        max_val_mem[i] = val;
        }
        
      const eT r_mean_val = r_mean_mem[i];
      const eT tmp        = val - r_mean_val;
    
      r_var_mem[i] = N_minus_1/N * r_var_mem[i] + (tmp*tmp)/N_plus_1;
      
      r_mean_mem[i] = r_mean_val + (val - r_mean_val)/N_plus_1;
      }
    }
  else
    {
    arma_debug_check( (sample.is_vec() == false), "running_stat_vec(): given sample is not a vector");
    
    x.r_mean.set_size(sample.n_rows, sample.n_cols);
    
    x.r_var.zeros(sample.n_rows, sample.n_cols);
    
    if(x.calc_cov == true)
      {
      x.r_cov.zeros(sample.n_elem, sample.n_elem);
      }
    
    x.min_val.set_size(sample.n_rows, sample.n_cols);
    x.max_val.set_size(sample.n_rows, sample.n_cols);
    
    
    const uword n_elem      = sample.n_elem;
    const eT*   sample_mem  = sample.memptr();
          eT*   r_mean_mem  = x.r_mean.memptr();
          eT*   min_val_mem = x.min_val.memptr();
          eT*   max_val_mem = x.max_val.memptr();
          
    
    for(uword i=0; i<n_elem; ++i)
      {
//.........这里部分代码省略.........

inline
void
glue_histc::apply_noalias(Mat<uword>& C, const Mat<eT>& A, const Mat<eT>& B, const uword dim)
  {
  arma_extra_debug_sigprint();
  
  arma_debug_check( ((B.is_vec() == false) && (B.is_empty() == false)), "histc(): parameter 'edges' is not a vector" );
  
  const uword A_n_rows = A.n_rows;
  const uword A_n_cols = A.n_cols;
  
  const uword B_n_elem = B.n_elem;
  
  if( B_n_elem == uword(0) )  { C.reset(); return; }
  
  const eT*   B_mem       = B.memptr();
  const uword B_n_elem_m1 = B_n_elem - 1;
  
  if(dim == uword(0))
    {
    C.zeros(B_n_elem, A_n_cols);
    
    for(uword col=0; col < A_n_cols; ++col)
      {
      const eT*    A_coldata = A.colptr(col);
            uword* C_coldata = C.colptr(col);
      
      for(uword row=0; row < A_n_rows; ++row)
        {
        const eT x = A_coldata[row];
        
        for(uword i=0; i < B_n_elem_m1; ++i)
          {
               if( (B_mem[i]           <= x) && (x < B_mem[i+1]) )  { C_coldata[i]++;           break; }
          else if(  B_mem[B_n_elem_m1] == x                      )  { C_coldata[B_n_elem_m1]++; break; }    // for compatibility with Matlab
          }
        }
      }
    }
  else
  if(dim == uword(1))
    {
    C.zeros(A_n_rows, B_n_elem);
    
    if(A.n_rows == 1)
      {
      const uword  A_n_elem = A.n_elem;
      const eT*    A_mem    = A.memptr();
            uword* C_mem    = C.memptr();
      
      for(uword j=0; j < A_n_elem; ++j)
        {
        const eT x = A_mem[j];
        
        for(uword i=0; i < B_n_elem_m1; ++i)
          {
               if( (B_mem[i]           <= x) && (x < B_mem[i+1]) )  { C_mem[i]++;           break; }
          else if(  B_mem[B_n_elem_m1] == x                      )  { C_mem[B_n_elem_m1]++; break; }    // for compatibility with Matlab
          }
        }
      }
    else
      {
      for(uword row=0; row < A_n_rows; ++row)
      for(uword col=0; col < A_n_cols; ++col)
        {
        const eT x = A.at(row,col);
        
        for(uword i=0; i < B_n_elem_m1; ++i)
          {
               if( (B_mem[i]            <= x) && (x < B_mem[i+1]) )  { C.at(row,i)++;           break; }
          else if(  B_mem[B_n_elem_m1]  == x                      )  { C.at(row,B_n_elem_m1)++; break; }   // for compatibility with Matlab
          }
        }
      }
    }
  }