Матрица расстояний с геосферой: избежать повторного исчисления
Если проблема заключается в дополнительном времени компиляции и раздувании бинарного размера, полученном путем компиляции .h как части всех используемых им .cpp-модулей, во многих случаях то, что вы можете сделать, - это сделать класс шаблона опустившимся из не-шаблонного базовый класс для не зависящих от типа частей интерфейса, и что базовый класс может иметь свою реализацию в файле .cpp.
3 ответа
Вы можете подготовить фрейм данных возможных комбинаций без повторений (с gtools пакетами). Затем рассчитать расстояния для этих пар. Вот код:
library(gtools)
library(geosphere)
library(data.table)
coords <- data.table(coordX = c(1, 2, 5, 9), coordY = c(2, 2, 0, 1))
pairs <- combinations(n = nrow(coords), r = 2, repeats.allowed = F, v = c(1:nrow(coords)))
distances <- apply(pairs, 1, function(x) {
distm(coords[x[1], ], coords[x[2], ], fun = distGeo)
})
# Construct distances matrix
dist_mat <- matrix(NA, nrow = nrow(coords), ncol = nrow(coords))
dist_mat[upper.tri(dist_mat)] <- distances
dist_mat[lower.tri(dist_mat)] <- distances
dist_mat[is.na(dist_mat)] <- 0
print(dist_mat)
Результаты:
[,1] [,2] [,3] [,4]
[1,] 0.0 111252.1 497091.1 400487.6
[2,] 111252.1 0.0 897081.9 786770.1
[3,] 497091.1 400487.6 0.0 458780.1
[4,] 897081.9 786770.1 458780.1 0.0
Если вы хотите вычислить все попарные расстояния для точек x, лучше использовать distm(x), а не distm(x,x). Функция distm возвращает одну и ту же симметричную матрицу в обоих случаях, но когда вы передаете ей один аргумент, она знает, что матрица симметрична, поэтому она не будет выполнять ненужные вычисления.
Вы можете рассчитать время.
library("geosphere")
n <- 500
xy <- matrix(runif(n*2, -90, 90), n, 2)
system.time( replicate(100, distm(xy, xy) ) )
# user system elapsed
# 61.44 0.23 62.79
system.time( replicate(100, distm(xy) ) )
# user system elapsed
# 36.27 0.39 38.05
Вы также можете посмотреть код R для geosphere::distm, чтобы проверить, что он обрабатывает два случая по-разному.
В стороне: Быстрый поиск в Google находит parallelDist: Матрица параллельного расстояния в CRAN. Геодезическое расстояние является опцией.
Использование combn() из базы R может быть немного проще и, вероятно, быстрее, чем загрузка дополнительных пакетов. Затем distm() использует distGeo() в качестве источника, поэтому использование последнего должно быть еще быстрее.
coords <- as.data.frame(coords) # this won't work with data.tables though
cbind(t(combn(1:4, 2)), unique(geosphere::distGeo(coords[combn(1:4, 2), ])))
# [,1] [,2] [,3]
# [1,] 1 2 111252.1
# [2,] 1 3 497091.1
# [3,] 1 4 897081.9
# [4,] 2 3 786770.1
# [5,] 2 4 400487.6
# [6,] 3 4 458780.1
Мы могли бы проверить это с помощью эталона.
Unit: microseconds
expr min lq mean median uq max neval cld
distm 555.690 575.846 597.7672 582.352 596.1295 904.718 100 b
distGeo 426.335 434.372 450.0196 441.516 451.8490 609.524 100 a
Выглядит хорошо.