为什么我的 Rcpp 代码会给出意外的所有整数值输出？

Question

给定一个数据矩阵 X，我想计算任意两行 X 之间的成对距离矩阵。我有以下代码，它来自于稍微调整代码 here。

#include <Rcpp.h>
#include <cmath>
#include <algorithm>

using namespace Rcpp;
// generic function for l1_distance
template <typename InputIterator1, typename InputIterator2>
inline double l1_distance(InputIterator1 begin1, InputIterator1 end1, 
                            InputIterator2 begin2) {
   double rval = 0;
   InputIterator1 it1 = begin1;
   InputIterator2 it2 = begin2;
   while (it1 != end1) {
      double d1 = *it1++;
      double d2 = *it2++;
      rval += abs(d1 - d2);
   }
   return rval;  
}

// [[Rcpp::export]]
NumericMatrix rcpp_l1_distance(NumericMatrix mat) {

   // allocate the matrix we will return
   NumericMatrix rmat(mat.nrow(), mat.nrow());
   for (int i = 0; i < rmat.nrow(); i++) {
      for (int j = 0; j < i; j++) {
         NumericMatrix::Row row1 = mat.row(i);
         NumericMatrix::Row row2 = mat.row(j);
         double d = l1_distance(row1.begin(), row1.end(), row2.begin());
         rmat(i,j) = d;
         rmat(j,i) = d;
      }
   }
   return rmat;
}

问题是这段代码返回了一个包含所有整数值的矩阵。整数值似乎与我想要的距离值正相关，这使得它更加混乱。我还计算了一个成对的 l2 距离矩阵和成对的标准化 l1 距离（将两行之间的 l1 距离除以它们的 l1 范数之和）矩阵，它们的行为都符合预期。

谁能告诉我是哪一部分出错了？

您可以执行以下操作以获得奇怪的结果

library(Rcpp)
sourceCpp("distance.cpp") #the file containing the cpp code above
X = matrix(rnorm(16), 4, 4)
rcpp_l1_distance(X)

提前致谢！

Answer 1

编译你的代码会给我这些警告：

> Rcpp::sourceCpp('Desktop/mat.cpp')
mat.cpp:16:15: warning: using integer absolute value function 'abs' when argument is of floating point type [-Wabsolute-value]
      rval += abs(d1 - d2);
              ^
mat.cpp:16:15: note: use function 'std::abs' instead
      rval += abs(d1 - d2);
              ^~~
              std::abs
mat.cpp:16:15: warning: using integer absolute value function 'abs' when argument is of floating point type [-Wabsolute-value]
      rval += abs(d1 - d2);
              ^
mat.cpp:30:18: note: in instantiation of function template specialization 'l1_distance<Rcpp::MatrixRow<14>::iterator, Rcpp::MatrixRow<14>::iterator>' requested here
      double d = l1_distance(row1.begin(), row1.end(), row2.begin());
                 ^
mat.cpp:16:15: note: use function 'std::abs' instead
      rval += abs(d1 - d2);
              ^~~
              std::abs
2 warnings generated.

...暗示abs是整数，见this help page，你可以使用fabs，或者std::abs，或者你可以使用sugar的operator-，abs和sum：

#include <Rcpp.h>

using namespace Rcpp;

// [[Rcpp::export]]
NumericMatrix rcpp_l1_distance(NumericMatrix mat) {

  // allocate the matrix we will return
  NumericMatrix rmat(mat.nrow(), mat.nrow());
  for (int i = 0; i < rmat.nrow(); i++) {
    NumericMatrix::Row row1 = mat.row(i);

    for (int j = 0; j < i; j++) {
      rmat(j,i) = rmat(i,j) = sum( abs(row1 - mat.row(j) )) ;
    }
  }
  return rmat;
}