Ближайшая точка к пути

Question

У меня есть два набора точек, называемых path а также centers, Для каждой точки в pathЯ хотел бы эффективный метод для нахождения идентификатора ближайшей точки в 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 являются координатами. Я хотел бы добавить столбец в path data.frame, который имеет идентификатор ближайшего центра для данной координаты x и y. Затем я хочу получить все уникальные идентификаторы.

Мое решение на данный момент работает, но очень медленно, когда масштаб проблемы увеличивается. Я хотел бы что-то гораздо более эффективное.

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 пакет. В моей системе это вычисляет ближайший center каждому из ваших path указывает менее чем за 2 секунды.

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

Вот еще один пример с 2,5 миллионами баллов-кандидатов, которые все попадают в пределы 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 

Добавление идентификатора точки к 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 сломан в выпущенной версии.

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 для создания входных данных для собственной матрицы фиксированного типа из заданного P dataframe.

Тестирование было выполнено с помощью пакета RcppParallel, но у объекта kd_tree нет конструктора копирования, поэтому необходимо было создать дерево для каждого потока, который съел бы все выгоды при параллельной обработке запросов.

RcppEigen и Rcpp11 в настоящее время не работают вместе, поэтому идея использования параллельного sapply Rcpp11 для запроса нелегко проверить.

---

// [[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;
}