R-在Rcpp中存储距离矩阵

debugcn 发表于 Dev

贾克

我很难弄清楚如何在中存储距离矩阵Rcpp。

让我们想象一下，我想将以下函数存储在n*n个人的距离矩阵中（sum由于不确定如何在中进行运算，所以不对进行平方）rcpp。

distxy = function(x,y) sum (x - y)

在此示例中，我想成对比较3个人

      [,1] [,2] [,3] [,4]
[1,]   24   24   22   20
[2,]   21   24   30   20
[3,]   44   34   41   13

在R我会通过这样的矩阵循环功能

mat = matrix(0, nrow(d), nrow(d))

  len = nrow(d)
  mat = matrix(0, len, len)

  for(j in 1:len){
    for(i in 1:len){
      mat[j,i] = distxy( d[j,], d[i,] ) 
    }
  }

并得到（我可以平方结果，但这在这里并不重要）

     [,1] [,2] [,3]
[1,]    0   -5  -42
[2,]    5    0  -37
[3,]   42   37    0

我在做同样的事时遇到麻烦 rcpp

到目前为止，我取得的成就是

// [[Rcpp::export]]
NumericVector FunCpp(NumericMatrix x) {
  int nrow = x.nrow();
  NumericMatrix out(nrow);

  for (int i = 0; i < nrow; i++) {
    for (int j = 0; j < nrow; j++) { 
      out[i,j] = sum( x(i,_) - x(j,_) ) ; 
    }
  }
  return out;
}

但是距离矩阵是不正确的。任何想法？

d = rbind(c(24, 24, 22, 20), 
      c(21, 24, 30, 20), 
      c(44, 34, 41, 13))

努尔塞尔

您的Rcpp代码中有两个语法错误：

返回aNumericVector而不是aNumericMatrix
使用operator[]由两个维度指数（out[i,j]）

这是一个清理的版本：

#include <Rcpp.h>

inline double distxy(Rcpp::NumericVector x, Rcpp::NumericVector y) {
    return Rcpp::sum(x - y);
}

// [[Rcpp::export]]
Rcpp::NumericMatrix FunCpp(Rcpp::NumericMatrix x) {
    int nrow = x.nrow();
    Rcpp::NumericMatrix out(nrow);

    for (int i = 0; i < nrow; i++) {
        for (int j = 0; j < nrow; j++) { 
            out(j, i) = distxy(x.row(j), x.row(i)); 
        }
    }

    return out;
}

针对您的R函数进行测试，

m <- matrix(
    c(24, 24, 22, 20,
      21, 24, 30, 20,
      44, 34, 41, 13),
    nrow = 3, byrow = TRUE
)

all.equal(FunR(m), FunCpp(m))
#[1] TRUE

至于平方，你可以使用std::pow里面distxy，

return std::pow(Rcpp::sum(x - y), 2);

或内部FunCpp循环中的内部：

out(j, i) = std::pow(distxy(x.row(j), x.row(i)), 2);

distxy <- function(x,y) sum(x - y) 

FunR <- function(d) {
    len <- nrow(d)
    mat <- matrix(0, len, len)

    for(j in 1:len){
        for(i in 1:len){
            mat[j,i] <- distxy(d[j,], d[i,]) 
        }
    }
    mat
}

本文收集自互联网，转载请注明来源。

如有侵权，请联系[email protected] 删除。