这里有一种使用
RcppArmadillo
的方法。很多部分与
RcppGallery 示例 非常相似。这将返回一个带有相关成对(按行)欧几里得距离的
big.matrix
。我喜欢用一个包装器函数来包装我的
big.matrix
函数,以创建更清晰的语法(即避免
@address
和其他初始化)。
注意-由于我们使用的是bigmemory(因此关注RAM使用情况),我让这个示例返回了仅包含下三角元素的 N-1 x N-1 矩阵。你可以修改这个,但这就是我拼凑出来的东西。
euc_dist.cpp
#define ARMA_NO_DEBUG
#include <RcppArmadillo.h>
using namespace Rcpp;
using namespace arma;
#include <bigmemory/BigMatrix.h>
template <typename T>
void BigArmaEuclidean(const Mat<T>& inBigMat, Mat<T> outBigMat) {
int W = inBigMat.n_rows;
for(int i = 0; i < W - 1; i++){
for(int j=i+1; j < W; j++){
outBigMat(j-1,i) = sqrt(sum(pow((inBigMat.row(i) - inBigMat.row(j)),2)));
}
}
}
void BigArmaEuc(SEXP pInBigMat, SEXP pOutBigMat) {
XPtr<BigMatrix> xpMat(pInBigMat);
XPtr<BigMatrix> xpOutMat(pOutBigMat);
int type = xpMat->matrix_type();
switch(type) {
case 1:
BigArmaEuclidean(
arma::Mat<char>((char *)xpMat->matrix(), xpMat->nrow(), xpMat->ncol(), false),
arma::Mat<char>((char *)xpOutMat->matrix(), xpOutMat->nrow(), xpOutMat->ncol(), false)
);
return;
case 2:
BigArmaEuclidean(
arma::Mat<short>((short *)xpMat->matrix(), xpMat->nrow(), xpMat->ncol(), false),
arma::Mat<short>((short *)xpOutMat->matrix(), xpOutMat->nrow(), xpOutMat->ncol(), false)
);
return;
case 4:
BigArmaEuclidean(
arma::Mat<int>((int *)xpMat->matrix(), xpMat->nrow(), xpMat->ncol(), false),
arma::Mat<int>((int *)xpOutMat->matrix(), xpOutMat->nrow(), xpOutMat->ncol(), false)
);
return;
case 8:
BigArmaEuclidean(
arma::Mat<double>((double *)xpMat->matrix(), xpMat->nrow(), xpMat->ncol(), false),
arma::Mat<double>((double *)xpOutMat->matrix(), xpOutMat->nrow(), xpOutMat->ncol(), false)
);
return;
default:
throw Rcpp::exception("Undefined type for provided big.matrix");
}
}
我的小包装
bigMatrixEuc <- function(bigMat){
zeros <- big.matrix(nrow = nrow(bigMat)-1,
ncol = nrow(bigMat)-1,
init = 0,
type = typeof(bigMat))
BigArmaEuc(bigMat@address, zeros@address)
return(zeros)
}
测试
library(Rcpp)
sourceCpp("euc_dist.cpp")
library(bigmemory)
set.seed(123)
mat <- matrix(rnorm(16), 4)
bm <- as.big.matrix(mat)
bm_out <- bigMatrixEuc(bm)[]
distMat <- as.matrix(dist(mat))
distMat[upper.tri(distMat, diag=TRUE)] <- 0
distMat <- distMat[2:4, 1:3]
all.equal(bm_out, distMat, check.attributes = FALSE)
[1] TRUE
bm_out
的矩阵。在阅读包装函数时,我原本以为bm_out
应该是一个big.matrix
。难道我错了吗?这个例子确实应该返回一个matrix
吗?有没有办法让bm_out
直接成为一个big.matrix
(而不是先得到一个matrix
再通过as.big.matrix
转换)? - Ricky