r语言 - 最接近路径的点

Question

我有两组点，分别是path和centers。对于path中的每个点，我想要一种有效的方法来查找centers中最近点的 ID。我想在 R 中执行此操作，下面是一个简单的可重现示例。

set.seed(1)
n <- 10000
x <- 100*cumprod(1 + rnorm(n, 0.0001, 0.002))
y <- 50*cumprod(1 + rnorm(n, 0.0001, 0.002))
path <- data.frame(cbind(x=x, y=y))
centers <- expand.grid(x=seq(0, 500,by=0.5) + rnorm(1001), 
                       y=seq(0, 500, by=0.2) + rnorm(2501))
centers$id <- seq(nrow(centers))

x和y是坐标。我想在 path data.frame 中添加一列，该列具有给定 x 和 y 坐标的最近中心的 id。然后我想获取所有的唯一 ID。

我目前的解决方案确实有效，但是当问题的规模增加时，速度会非常慢。我想要更高效的东西。

path$closest.id <- sapply(seq(nrow(path)), function(z){
   tmp <- ((centers$x - path[z, 'x'])^2) + ((centers$y - path[z, 'y'])^2)
   as.numeric(centers[tmp == min(tmp), 'id'])
})
output <- unique(path$closest.id)

任何关于加快速度的帮助将不胜感激。

我认为data.table可能会有所帮助，但理想情况下，我正在寻找的是一种在搜索方面可能更聪明的算法，即而不是计算到每个中心的距离，然后只选择最小中心......要获取 ID...

我也很乐意使用 Rcpp/Rcpp11 如果这有助于提高性能。

我执行这种计算的最小可接受时间是 10 秒，但显然越快越好。

Answer 1

您可以使用

RANN包中的nn2执行此操作。在我的系统上，这会在 2 秒内计算出最接近每个path点的center。

library(RANN)
system.time(closest <- nn2(centers[, 1:2], path, 1))
#   user  system elapsed 
#   1.41    0.14    1.55 

sapply(closest, head)
#      nn.idx   nn.dists
# [1,] 247451 0.20334929
# [2,] 250454 0.12326323
# [3,] 250454 0.28540127
# [4,] 253457 0.05178687
# [5,] 253457 0.13324137
# [6,] 253457 0.09009626

下面是另一个示例，其中包含 250 万个候选点，这些候选点都落在path点的范围内(在您的示例中，centers的x和y范围比path点大得多(。在这种情况下，它有点慢。

set.seed(1)
centers2 <- cbind(runif(2.5e6, min(x), max(x)), runif(2.5e6, min(y), max(y)))
system.time(closest2 <- nn2(centers2, path, 1))
#   user  system elapsed 
#   2.96    0.11    3.07 
sapply(closest2, head)
#       nn.idx    nn.dists
# [1,]  730127 0.025803703
# [2,]  375514 0.025999069
# [3,] 2443707 0.047259283
# [4,]   62780 0.022747930
# [5,] 1431847 0.002482623
# [6,] 2199405 0.028815865

这可以与使用 sp::spDistsN1 的输出进行比较(对于此问题来说，输出要慢得多(：

library(sp)
apply(head(path), 1, function(x) which.min(spDistsN1(centers, x)))
#       1       2       3       4       5       6 
#  730127  375514 2443707   62780 1431847 2199405 

将点 id 添加到 data.frame path并简化为唯一值是微不足道的：

path$closest.id <- closest$nn.idx
output <- unique(path$closest.id)

Answer 2

这是一个

Rcpp11的解决方案。类似的东西可能适用于Rcpp，但有一些变化。

#define RCPP11_PARALLEL_MINIMUM_SIZE 1000
#include <Rcpp11>
inline double square(double x){
    return x*x ;
}
// [[Rcpp::export]]
IntegerVector closest( DataFrame path, DataFrame centers ){
    NumericVector path_x = path["x"], path_y = path["y"] ;
    NumericVector centers_x = centers["x"], centers_y = centers["y"] ;
    int n_paths = path_x.size(), n_centers = centers_x.size() ; 

    IntegerVector ids = sapply( seq_len(n_paths), [&](int i){
            double px = path_x[i], py=path_y[i] ;
            auto get_distance = [&](int j){
                return  square(px - centers_x[j]) + square(py-centers_y[j]) ;
            } ;
            double distance = get_distance(0) ;
            int res=0;
            for( int j=1; j<n_centers; j++){
                double d = get_distance(j)  ;
                if(d < distance){
                    distance = d ;
                    res = j ;
                }
            }
            return res + 1 ;
    }) ;
    return unique(ids) ;
}

我得到 ：

> set.seed(1)
> n <- 10000
> x <- 100 * cumprod(1 + rnorm(n, 1e-04, 0.002))
> y <- 50 * cumprod(1 + rnorm(n, 1e-04, 0.002))
> path <- data.frame(cbind(x = x, y = y))
> centers <- expand.grid(x = seq(0, 500, by = 0.5) +
+     rnorm(1001), y = seq(0, 500, by = 0.2) + rnorm(2501))
> system.time(closest(path, centers))
   user  system elapsed
 84.740   0.141  21.392

这利用了糖的自动并行化，即 sapply并行运行。 #define RCPP11_PARALLEL_MINIMUM_SIZE 1000部分是强制并行，否则默认情况下仅从 10000 开始。但在这种情况下，由于内部计算非常耗时，因此值得。

请注意，您需要 Rcpp11 的开发版本，因为已发布版本中unique已损坏。

Answer 3

此解决方案将示例数据集的处理时间缩短了近 RANN 解决方案的一半。

可以使用devtools::install_github("thell/Rcppnanoflann")

安装

Rcppnanoflann 解决方案利用了 Rcpp、RcppEigen 和nanoflann EigenMatrixAdaptor 以及 c++11 以产生与原始问题相同的唯一索引。

library(Rcppnanoflann)
system.time(o.nano<-nnIndex(centers,path))
##    user  system elapsed 
##    0.62    0.05    0.67

^{_{* 使用原始问题中定义的路径和中心值}}

为了获得与原始样品相同的结果，RANN 解决方案需要我们在这里稍作修改...

library(RANN)
system.time(o.flann<-unique(as.numeric(nn2(centers,path,1)$nn.idx)))
##    user  system elapsed 
##    1.24    0.07    1.30

identical(o.flann,o.nano)
## [1] TRUE

Rcppnanoflann的工作功能利用了Eigen's Map能够从给定P数据帧。

测试是使用 RcppParallel 包完成的，但kd_tree对象没有一个复制构造函数，因此需要为每个线程创建树这消耗了并行查询处理中的任何收益。

RcppEigen 和 Rcpp11 目前不一起玩，所以想法使用 Rcpp11 的并行 sapply 进行查询并不容易测试。

---

// [[Rcpp::export]]
std::vector<double> nnIndex(const Rcpp::DataFrame & P, const Rcpp::DataFrame & Q )
{
  using namespace Eigen;
  using namespace Rcpp;
  using namespace nanoflann;
  // Matrix of points to be queried against.
  const NumericVector & Px(P[0]);
  const NumericVector & Py(P[1]);
  MatrixX2d M(Px.size(), 2);
  M.col(0) = VectorXd::Map(&Px[0],Px.size());
  M.col(1) = VectorXd::Map(&Py[0],Py.size());
  // The points to query.
  const NumericVector & Qx(Q[0]);
  const NumericVector & Qy(Q[1]);
  double query_pt[2];
  size_t query_count(Qx.size());
  // Populate a 2d tree.
  KD_Tree kd_tree( 2, M, 10 );
  kd_tree.index->buildIndex();
  std::set<size_t> nn;
  std::vector<double> out;
  out.reserve(query_count);
  size_t index(0);
  double quadrance;
  for( size_t i=0 ; i < query_count; ++i ) {
    query_pt[0] = Qx[i];
    query_pt[1] = Qy[i];
    kd_tree.index->knnSearch( &query_pt[0],1, &index, &quadrance);
    if( nn.emplace(index).second ) out.emplace_back(index+1);
  }
  return out;
}