最接近路径的点

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。然后，我想获得所有唯一的ids。

我目前的解决办法确实有效，但当问题的规模增加时，我的解决办法非常缓慢。我想要更有效率的东西。

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可能会有所帮助，但理想情况下，我要寻找的是一种在搜索方面可能更聪明的算法，即不是计算到每个中心的距离，然后只选择最小的一个.为了拿到身份证..。

我也很高兴使用Rcpp/Rcpp11，如果这有助于提高性能的话。

我完成这种计算的最小可接受时间是10秒，但显然更快会更好。

Answer 1

您可以使用来自nn2包的RANN来完成这一任务。在我的系统中，这将在2秒内计算到每个path点的最近的path。

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添加到path data.frame并将其简化为唯一值是非常简单的：

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在发布的版本中被破坏了。