data.table里面的rcpp函数加入

Rcpp function inside data.table join

问题

为什么我的 Rcpp 函数在 data.table 连接中产生的结果与在连接外使用时不同(并且不正确)?

例子

我有两个 data.table,我想找出两个表中每对坐标之间的欧氏距离。

为了进行距离计算,我定义了两个函数,一个在 base R 中,另一个使用 Rcpp

library(Rcpp)
library(data.table)

rEucDist <- function(x1, y1, x2, y2) return(sqrt((x2 - x1)^2 + (y2 - y1)^2))

cppFunction('NumericVector cppEucDistance(NumericVector x1, NumericVector y1,
                                          NumericVector x2, NumericVector y2){

            int n = x1.size();
            NumericVector distance(n);

            for(int i = 0; i < n; i++){
              distance[i] = sqrt(pow((x2[i] - x1[i]), 2) + pow((y2[i] - y1[i]), 2));
            }
            return distance;
}')


dt1 <- data.table(id = rep(1, 6),
                  seq1 = 1:6,
                  x = c(1:6),
                  y = c(1:6))

dt2 <- data.table(id = rep(1, 6),
                  seq2 = 7:12,
                  x = c(6:1),
                  y = c(6:1))

先做连接,再计算距离时,两个函数产生相同的结果

dt_cpp <- dt1[ dt2, on = "id", allow.cartesian = T]
dt_cpp[, dist := cppEucDistance(x, y, i.x, i.y)]

dt_r <- dt1[ dt2, on = "id", allow.cartesian = T]
dt_r[, dist := rEucDist(x, y, i.x, i.y)]

all.equal(dt_cpp$dist, dt_r$dist)
# [1] TRUE

但是,如果我在联接内进行计算,结果会有所不同; cpp 版本不正确。

dt_cppJoin <- dt1[
    dt2,
    { (cppEucDistance(x, y, i.x, i.y)) },
    on = "id",
    by = .EACHI
    ]

dt_rJoin <- dt1[
    dt2,
    { (rEucDist(x, y, i.x, i.y)) },
    on = "id",
    by = .EACHI
    ]

all.equal(dt_cppJoin$V1, dt_rJoin$V1)
# "Mean relative difference: 0.6173913"

## note that the R version of the join is correct
all.equal(dt_r$dist, dt_rJoin$V1)
# [1] TRUE

Rcpp 实现导致连接版本给出不同结果的原因是什么?

tl;博士

  1. 最终,这归结为我对 EACHI 的工作原理缺乏了解。

  2. 我的 C++ 循环依赖于 x1x2 的长度相同。使用 EACHI 时情况并非如此,因此我在 C++ 函数中的向量子集设置将不正确。因此我需要一个更复杂的 C++ 函数。


我相信我看到的 'issue' 归结为 EACHI 正在做的事情。看了一遍又一遍Arun's answer,我想我明白了。

特别是 EACHI "evaluates j-expression on the matching rows for each row in Y"。所以它采用 Y table (在我的例子中是 dt2 )并一次评估它的每一行。如果你在Rcpp函数里面使用几个prints就可以看出这一点

## I'm reducing the data sets so the 'prints' are readible
dt1 <- data.table(id = rep(1, 2),
                  seq1 = 1:2,
                  x = c(1:2),
                  y = c(1:2))

dt2 <- data.table(id = rep(1, 2),
                  seq2 = 7:8,
                  x = c(6:5),
                  y = c(6:5))

cppFunction('NumericVector cppEucDistance(NumericVector x1, NumericVector y1,
                                          NumericVector x2, NumericVector y2){

  int n = x1.size();
  NumericVector distance(n);

  Rcpp::Rcout << "x1 size: " << x1.size() << std::endl;
  Rcpp::Rcout << "x2 size: " << x2.size() << std::endl;
  Rcpp::Rcout << "n size: " << n << std::endl;

  for(int i = 0; i < n; i++){
    distance[i] = sqrt(pow((x2[i] - x1[i]), 2) + pow((y2[i] - y1[i]), 2));
  }
  return distance;
}')

dt_cppJoin <- dt1[
  dt2,
  {print(i.x); (cppEucDistance(x, y, i.x, i.y)) },
  on = "id",
  by = .EACHI
  ]
# [1] 6
# x1 size: 2
# x2 size: 1
# n size: 2
# [1] 5
# x1 size: 2
# x2 size: 1
# n size: 2

请注意打印语句的输出,x2 每次迭代的长度仅为 1。而 n 是 2。所以当 i 到达第二个元素时,我的 x2[i] 显然会 'fail'(记住 C++ 数组是从 0 开始索引的)

而去掉 EACHI 与全连接方法具有相同的效果

dt_cppJoin <- dt1[
  dt2,
  {print(i.x); (cppEucDistance(x, y, i.x, i.y)) },
  on = "id"
  ]
# [1] 6 6 5 5
# x1 size: 4
# x2 size: 4
# n size: 4