Rcpp中按行和列名的NumericMatrix子集

Question

我正试图在Rcpp中创建一个函数，该函数将接受成对数字矩阵以及向量列表作为输入，每个元素都是行/列名的子集。我希望这个函数识别与这些名称匹配的矩阵的子集，并返回值的平均值。

下面我生成了一些类似于我所拥有的数据类型的伪数据，并尝试使用Rcpp函数。

library(Rcpp)
dat <- c(spA = 4, spB = 10, spC = 8, spD = 1, spE = 5, spF = 9)
pdist <- as.matrix(dist(dat))
pdist[upper.tri(pdist, diag = TRUE)] <- NA

这里我有一个由pdist中行/列名的不同子集的字符向量组成的列表

subsetList <- replicate(10, sample(names(dat), 4), simplify=FALSE)

对于这些名称集合中的每一个，我想确定成对矩阵的子集，并取值的平均值

这是我迄今为止所拥有的，但并不奏效，但我认为它说明了我正在努力实现的目标。

cppFunction('
List meanDistByCell(List input, NumericMatrix pairmat) {
int n = input.size();
List out(n);
List dimnames = pairmat.attr( "dimnames" );
CharacterVector colnames = dimnames[1];
for (int i = 0; i < n; i++) {
CharacterVector sp = as< CharacterVector >(input[i]);
if (sp.size() > 0) {
out[i] = double(mean(pairmat(sp, sp)));
} else {
out[i] = NA_REAL;
}
}
return out;
}
')

如有任何帮助，我们将不胜感激！谢谢

Answer 1

尽管可以使用(连续的)基于范围的子集设置(例如x(Range(first_row, last_row), Range(first_col, last_col)))，但正如coatless所指出的，目前不支持通过CharacterVector进行子集设置，因此您暂时必须自己进行。一种通用的方法可能看起来像这样：

template <int RTYPE> inline Matrix<RTYPE>
Subset2D(const Matrix<RTYPE>& x, CharacterVector crows, CharacterVector ccols) {
R_xlen_t i = 0, j = 0, rr = crows.length(), rc = ccols.length(), pos;
Matrix<RTYPE> res(rr, rc);
CharacterVector xrows = rownames(x), xcols = colnames(x);
IntegerVector rows = match(crows, xrows), cols = match(ccols, xcols);
for (; j < rc; j++) {
// NB: match returns 1-based indices
pos = cols[j] - 1;
for (i = 0; i < rr; i++) {
res(i, j) = x(rows[i] - 1, pos);
}
}
rownames(res) = crows;
colnames(res) = ccols;
return res;
}
// [[Rcpp::export]]
NumericMatrix subset2d(NumericMatrix x, CharacterVector rows, CharacterVector cols) {
return Subset2D(x, rows, cols);
}

这假设输入矩阵同时具有行和列名，并且行和列查找向量是这些dimname的有效子集；可以添加额外的防御代码以使其更加健壮。为了证明，

subset2d(pdist, subsetList[[1]], subsetList[[1]])
#     spB spD spE spC
# spB  NA  NA  NA  NA
# spD   9  NA  NA   7
# spE   5   4  NA   3
# spC   2  NA  NA  NA
pdist[subsetList[[1]], subsetList[[1]]]
#     spB spD spE spC
# spB  NA  NA  NA  NA
# spD   9  NA  NA   7
# spE   5   4  NA   3
# spC   2  NA  NA  NA

---

Subset2D负责实现meanDistByCell所涉及的大部分样板；剩下的就是在输入列表上循环，将其应用于每个列表元素，并将结果的平均值存储在输出列表中：

// [[Rcpp::export]]
List meanDistByCell(List keys, NumericMatrix x, bool na_rm = false) {
R_xlen_t i = 0, sz = keys.size();
List res(sz);
if (!na_rm) {
for (; i < sz; i++) {
res[i] = NumericVector::create(
mean(Subset2D(x, keys[i], keys[i]))
);
}
} else {
for (; i < sz; i++) {
res[i] = NumericVector::create(
mean(na_omit(Subset2D(x, keys[i], keys[i])))
);
}
}
return res;
}

all.equal(
lapply(subsetList, function(x) mean(pdist[x, x], na.rm = TRUE)),
meanDistByCell2(subsetList, pdist, TRUE)
)
# [1] TRUE

---

尽管使用Subset2D可以实现更干净的meanDistByCell，但在这种情况下，它的效率低下，至少有几个原因：

• 它设置返回对象(rownames(res) = crows;，colnames(res) = ccols;)的dimname，您在这里不需要它们
• 它调用match以获得rownames和colnames的每个的索引，这是不必要的，因为您事先知道rownames(x) == colnames(x)

对于长度为k的输入列表，您将承担这两个点的成本k次。

一种更有效但不那么简洁的方法是基本上只实现Subset2D所需的方面，在meanDistByCell内部内联：

// [[Rcpp::export]]
List meanDistByCell2(List keys, NumericMatrix x, bool na_rm = false) {
R_xlen_t k = 0, sz = keys.size(), i = 0, j = 0, nidx, pos;
List res(sz);
CharacterVector cx = colnames(x);
if (!na_rm) {
for (; k < sz; k++) {
// NB: match returns 1-based indices
IntegerVector idx = match(as<CharacterVector>(keys[k]), cx) - 1;
nidx = idx.size();
NumericVector tmp(nidx * nidx);
for (j = 0; j < nidx; j++) {
pos = idx[j];
for (i = 0; i < nidx; i++) {
tmp[nidx * j + i] = x(idx[i], pos);
}
}
res[k] = NumericVector::create(mean(tmp));
}
} else {
for (; k < sz; k++) {
IntegerVector idx = match(as<CharacterVector>(keys[k]), cx) - 1;
nidx = idx.size();
NumericVector tmp(nidx * nidx);
for (j = 0; j < nidx; j++) {
pos = idx[j];
for (i = 0; i < nidx; i++) {
tmp[nidx * j + i] = x(idx[i], pos);
}
}
res[k] = NumericVector::create(mean(na_omit(tmp)));
}
}
return res;
}

all.equal(
meanDistByCell(subsetList, pdist, TRUE),
meanDistByCell2(subsetList, pdist, TRUE)
)
# [1] TRUE

---