0

我正在尝试在 Rcpp 中实现一个简单的代码来计算和填充距离矩阵的条目。问题是 Rcpp 代码（如下）返回一个矩阵 D ，其中所有元素的值都为零。这个问题似乎没有在论坛的任何地方得到解决 - 我很感激一些建议！

  src_d_err_c <- '
  using namespace Rcpp;
  double d_err_c(NumericVector cx, NumericVector csx, NumericVector cy, NumericVector csy) {
    using namespace Rcpp;    
    NumericVector d = (cx - cy)*(cx - cy) / csx;
    double s = std::accumulate(d.begin(), d.end(), 0.0);

    return s;
  }'

  src_d_mat = '   
  using namespace Rcpp;
  // input

  Rcpp::NumericMatrix cX(X);
  Rcpp::NumericMatrix cY(Y);
  Rcpp::NumericMatrix cSX(SX);
  Rcpp::NumericMatrix cSY(SY);  
  int N1 = cX.nrow();
  int N2 = cY.nrow();
  NumericMatrix D(N1, N2);
  NumericVector v(N1);

  for (int x = 0; x++; x<N1){
    v[x] = x;
    for (int y = 0; y++; y<N2) {
      D(x,y) = d_err_c(cX(x,_), cSX(x,_), cY(y,_), cSY(y,_));
    };
  };
  return wrap(v);'  
  fun <- cxxfunction(signature(X = "numeric",  SX = "numeric", 
                               Y = "numeric",  SY = "numeric"),
                     body = src_d_mat, includes = src_d_err_c, 
                     plugin = "Rcpp")
4

2 回答 2

4

循环的参数for顺序错误：条件应该在中间，增量在末尾。

for (int x = 0; x < N1; x++)
于 2013-09-15T11:24:41.150 回答
1

@Vincent 正确地指出了一个主要错误（实际上不是循环），但还有另一个主要错误：v当您的计算进入D您想要返回的内容时，您返回。（实际上，您也根本不需要v。）

这是一个修复版本，它使用正确的循环索引，返回正确的对象，并省略未使用的代码。它还切换到 Rcpp 属性。

将其保存在一个文件中，例如“distmat.cpp”，然后使用sourceCpp("distmat.cpp")，之后d_mat您可以调用一个新函数。

#include <Rcpp.h>

using namespace Rcpp;

double d_err_c(NumericVector cx, NumericVector csx, 
               NumericVector cy, NumericVector csy) {
    NumericVector d = (cx - cy)*(cx - cy) / csx;
    return std::accumulate(d.begin(), d.end(), 0.0);
}

// [[Rcpp::export]]
NumericMatrix d_mat(NumericMatrix cX, NumericMatrix cY, 
                    NumericMatrix cSX, NumericMatrix cSY) {
    int N1 = cX.nrow();
    int N2 = cY.nrow();
    NumericMatrix D(N1, N2);

    for (int x = 0; x<N1; x++) {
        for (int y = 0; y<N2; y++) {
            D(x,y) = d_err_c(cX(x,_), cSX(x,_), cY(y,_), cSY(y,_));
        }
    }
    return D;
}
于 2013-09-15T13:48:44.893 回答