这篇博文的目的是展示R语言中下标操作矩阵的潜在问题。R语言提供了多种方法提取一个矩阵的单个或者部分元素,不同方法对应的效率在Hadley Wickham的Advance R中已有讨论。这些方法中,使用最广泛的是通过下标(行或者列)取值,即操作符[。然而,这种方法存在潜在问题,即内存中会拷贝原始对象。

举例:首先建立一个矩阵,之后取这个矩阵除了第一行之外的部分,接下来操作这个部分矩阵。

manipulate
1
2
3
4
5
6
7
8
9
10
11
12
## step1: build matrix
n <- 8000
tmp1 <- matrix(rnorm(n * n), nrow = n, ncol = n)
gc()

## step2: manipulate a subset of matrix
sink('/dev/null')
apply(tmp1[2:n, ], 1, function(x) x[1])
sink()

## step3: garbage collection
gc()

1. C语言指针基础

C语言的指针设计是一致和优雅的。C语言中“指针(pointer)”就是地址(所以不能用普通整数储存地址),“指针变量(pointer variable)”是存储地址的变量。一个指针变量,只能指向一个特定类型的变量,比如整数、浮点数、字符或者指针。

Initiate a pointer
1
2
3
4
5
6
7
8
9
10
11
12
13
14
int tmp1 = 1, tmp2;

/* "=" does not mean "assignment", it just means "initiating" */
/* p is the address of tmp1, *p is equal to the value of tmp1*/
int *p = &tmp1;

int *q;
q = &tmp2;

/* p points to tmp1, q points to tmp2, now the value of tmp2 is 1*/
*p = *q;

/* p and q now both points to tmp1*/
q = p;

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

1. Rcpp

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

2. RcppArmadillo

基本类型是matveccolvec)和rowvec

  • 属性

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

    • 全是0ones<mat>(3, 4);/vec(10, fill::ones);/;全是1zeros<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)

参考网址

更新记录

2017年1月15日

XPath提供了一种对XML节点、节点属性和内容快速查询的规则。在各种编程语言中都有实现,比如C语言的libxml2和对应的R包 xml2

1. 查询规则

XPath查询集中在三个对象:节点、节点属性和节点内容。

1.1 选择节点

  • /nodeA/nodeB:nodeA为根节点,nodeA下的所有nodeB节点;等价于nodeB

  • //nodeB:所有nodeB节点,在R包xml2中(比如函数xml_find_all()),//nodeB搜索范围是整个文档,忽略当前节点;而.//nodeB搜索范围是当前节点之下。