收集和记录Rcpp或者RcppArmadillo操作矩阵和向量。

1. Rcpp

  • 可以使用逻辑下标(LogicalVector)对向量和列表取值

2. RcppArmadillo

基本类型是matveccolvec)和rowvec

  • 属性

    • 对于矩阵,行数:m.n_rows;;列数:m.n_cols;;维度:m.size();size(m);。对于向量,元素数:v.n_elem;
  • 特殊向量或矩阵

    • 全是1ones<mat>(3, 4);/vec(10, fill::ones);/;全是0zeros<vec>(10);/mat(3, 4, fill::zeros);;全是某个数mat a(4, 5); a.fill(123.4);

    • 连续向量,规定长度linspace<vec>(0, 5, 6);;连续向量,规定间距regspace<vec>(0, 2, 9);

  • 取值

    • 对于向量,连续取值:v.subvec(1stIdx, lastIdx);;非连续,可以考虑使用find()函数,比如:v.elem(find(v > 0));

    • 对于矩阵,连续取值:m.col(Idx);/m.row(Idx);/m.cols(Idx);/m.rows(Idx);/m.submat(1stRowIdx, lastRowIdx, 1stColIdx, lastColIdx);;非连续,m.submat(vecRowIdx, vecColIdx);

  • Rcpp矩阵转换为RcppArmadillo矩阵,可以避免拷贝矩阵,以提升效率,mat(ptr_aux_mem, n_rows, n_cols, copy_aux_mem = true, strict = false)。同样道理,可以转化向量。例如:

transfer matrix and vector
1
2
3
4
5
6
7
8
9
10
11
12
13
arma::mat TransferMatArma(Rcpp::NumericMatrix x, Rcpp::NumericVector y) {
    mat tx(x.begin(), x.nrow(), x.ncol(), false);
    vec ty(y.begin(), y.size(), false);
    return tx;
}

Rcpp::NumericVector TransferMatRcpp(arma::mat x, arma::vec y) {
    NumericMatrix tx(x.n_rows, x.n_cols, x.begin());
    NumericVector ty(y.begin(), y.end());
    return ty;

// 不要使用as<IntegerVector>(wrap(y)),会有内存泄露。
}
  • 使用.each_col()/.each_row()/.for_each()替代apply()
replace apply()
1
2
3
4
5
6
7
8
9
10
11
12
arma::mat TestMat(arma::mat M, double a) {

  M.for_each([a](mat::elem_type& val) {
      val = val > 0 ? val : a;
    });

  M.each_row([](rowvec& r) {
      r /= r.max();
    });

  return M;
}
  • 使用sum(M, 0);sum(M, 1);分别替代colSums(M)rowSums(M)

3. bigmemory

bigmemory包提供了四种数据类型的矩阵,即double(默认)、integershortchar。对于big.matrix对象pMat,四种类型通过通过MatrixAccessor<double> macc(*pMat)MatrixAccessor<int> macc(*pMat)MatrixAccessor<short> macc(*pMat)MatrixAccessor<char> macc(*pMat)提取元素。pMat有三种属性nrow()ncol()matrix_type()可以使用。一下代码示例展示了将big.matrix转换为matrix

manipulate big.matrix with cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
#include <Rcpp.h>
// [[Rcpp::depends(BH, bigmemory)]]
#include <bigmemory/MatrixAccessor.hpp>

#include <numeric>

using namespace Rcpp;


// [[Rcpp::export]]
Rcpp::NumericMatrix TestBigMat(XPtr<BigMatrix> pMat) {

  MatrixAccessor<int> macc(*pMat);

  int n = pMat->nrow();
  int m = pMat->ncol();

  NumericMatrix resMat(n, m);

  for (int i = 0; i < n; ++i) {
    for (int j = 0; j < m; ++j) {
      resMat(i, j) = macc[j][i];
    }
  }

  return resMat;
}

注意事项:

  • 获取元素为列-行形式,因为big.matrix按照列存储矩阵。例如macc[j][i]表示i-1行的j-1元素。

  • 调用函数使用big.matrix的地址,例如TestBigMat(bigmat@address)

同样,RcppArmadillo也能与bigmemory结合,例如:

manipulate big.matrix with armadillo
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
#include <RcppArmadillo.h>
#include <bigmemory/MatrixAccessor.hpp>

#include <numeric>

// [[Rcpp::depends(RcppArmadillo, bigmemory)]]

using namespace Rcpp;
using namespace arma;

// [[Rcpp::export]]
Rcpp::NumericMatrix TestBigMatArma(SEXP pMat) {

  XPtr<BigMatrix> xpMat(pMat);

  Mat<int> macc = Mat<int>((int *)xpMat->matrix(), xpMat->nrow(), xpMat->ncol(), false);

  int n = xpMat->nrow();
  int m = xpMat->ncol();

  NumericMatrix resMat(n, m);

  for (int i = 0; i < n; ++i) {
    for (int j = 0; j < m; ++j) {
      resMat(i, j) = macc(i, j);
    }
  }

  return resMat;
}

参考网址

更新记录

2018年9月17日

Comments