5

我有两组点，称为path和centers。对于 中的每个点path，我想要一种有效的方法来查找中最近点的 ID centers。我想在 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是坐标。我想向pathdata.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 秒，但显然更快会更好。

4

3 回答 3

13

您可以nn2从RANN包中执行此操作。在我的系统上，这会在 2 秒内计算出离center您的每个点最近的点。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大得多）。在这种情况下它会慢一点。xypath

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 添加到pathdata.frame 并减少到唯一值是微不足道的：

path$closest.id <- closest$nn.idx
output <- unique(path$closest.id)
于 2014-12-08T14:47:45.603 回答
6

这是一个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在发布版本中已损坏。

于 2014-12-08T16:50:35.657 回答
1

该解决方案将样本数据集的处理时间几乎减少了 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 的 Map 功能从给定的数据帧创建固定类型 Eigen 矩阵的输入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;
}
于 2014-12-08T01:20:08.663 回答