Блог

Главная - Вопросы по программированию - эффективность вычисления евклидова расстояния между точками и опорными точками в r

эффективность вычисления евклидова расстояния между точками и опорными точками в r

#r #dplyr #euclidean-distance

Вопрос:

У меня есть список точек a1, a2, b1 и b2 и набор опорных точек. Я хотел бы вычислить расстояние между четырьмя точками и всеми контрольными точками. Я могу сделать это без проблем, используя следующий пример

 library(dplyr)

x <- rep(1:5, 5)
y <- rep(1:5, each = 5)  
a1.x <- c(4.5)
a1.y <- c(6)
a2.x <- c(0.8)
a2.y <- c(3.2)
b1.x <- c(2.5)
b1.y <- c(5)
b2.x <- c(3.8)
b2.y <- c(1.5)
time <- 1

time1 <- as.data.frame(cbind(time,x,y,a1.x,a1.y,a2.x,a2.y,b1.x,b1.y,b2.x,b2.y))

x <- rep(1:5, 5)
y <- rep(1:5, each = 5)  
a1.x <- c(4)
a1.y <- c(5)
a2.x <- c(1.5)
a2.y <- c(3.9)
b1.x <- c(1.4)
b1.y <- c(4.6)
b2.x <- c(6)
b2.y <- c(5.2)
time <- 2


time2 <- as.data.frame(cbind(time,x,y,a1.x,a1.y,a2.x,a2.y,b1.x,b1.y,b2.x,b2.y))

df <- rbind(time1,time2)

df <- df %>% 
  mutate(dista1 = sqrt((x-a1.x)^2   (y-a1.y)^2)) %>% 
  mutate(dista2 = sqrt((x-a2.x)^2   (y-a2.y)^2)) %>% 
  mutate(distb1 = sqrt((x-b1.x)^2   (y-b1.y)^2)) %>% 
  mutate(distb2 = sqrt((x-b2.x)^2   (y-b2.y)^2))

Это работает без проблем и довольно быстро. Однако в большем наборе данных это замедляется за счет ввода всех столбцов или необходимости полагаться на циклы. Каков наиболее эффективный способ выполнения вышеперечисленных действий?

Редактировать — кроме того, у меня также есть фактор времени

Ответ №1:

Использование rdist функции в пакете fields проще:

 library(fields)   #use install.packages("fields") first
pts <- cbind(x= c(a1.x, a2.x, b1.x, b2.x), y=c(a1.y, a2.y, b1.y, b2.y))
ref <- cbind(x, y)
distances <- rdist(ref, pts)

colnames(distances) <- c("dista1", "dista2", "distb1", "distb2")
head(distances)
#        dista1   dista2   distb1    distb2
# [1,] 6.103278 2.209072 4.272002 2.8442925
# [2,] 5.590170 2.505993 4.031129 1.8681542
# [3,] 5.220153 3.111270 4.031129 0.9433981
# [4,] 5.024938 3.883298 4.272002 0.5385165
# [5,] 5.024938 4.741308 4.716991 1.3000000
# [6,] 5.315073 1.216553 3.354102 2.8442925

Если вы хотите сопоставить df в своем примере:

 df <- cbind(ref, a1.x, a1.y, a2.x, a2.y, b1.x, b1.y, b2.x, b2.y, distances)
head(df)
 #     x y a1.x a1.y a2.x a2.y b1.x b1.y b2.x b2.y   dista1   dista2   distb1    distb2
# [1,] 1 1  4.5    6  0.8  3.2  2.5    5  3.8  1.5 6.103278 2.209072 4.272002 2.8442925
# [2,] 2 1  4.5    6  0.8  3.2  2.5    5  3.8  1.5 5.590170 2.505993 4.031129 1.8681542
# [3,] 3 1  4.5    6  0.8  3.2  2.5    5  3.8  1.5 5.220153 3.111270 4.031129 0.9433981
# [4,] 4 1  4.5    6  0.8  3.2  2.5    5  3.8  1.5 5.024938 3.883298 4.272002 0.5385165
# [5,] 5 1  4.5    6  0.8  3.2  2.5    5  3.8  1.5 5.024938 4.741308 4.716991 1.3000000
# [6,] 1 2  4.5    6  0.8  3.2  2.5    5  3.8  1.5 5.315073 1.216553 3.354102 2.8442925

Если есть несколько раз, этот подход может быть расширен. Сначала извлеките из ваших time1 time2 объектов и, чтобы создать несколько точек и опорных матриц:

 time1.pts <- matrix(unlist(time1[1, 4:11]), 4, 2, byrow=TRUE)
time2.pts <- matrix(unlist(time2[1, 4:11]), 4, 2, byrow=TRUE)
ref1 <- matrix(unlist(time1[1, 2:3]), 4, 2, byrow=TRUE)
ref2 <- matrix(unlist(time2[1, 2:3]), 4, 2, byrow=TRUE)
ref <- list(ref1=ref1, ref2=ref2)
pts <- list(time1.pts=time1.pts, time2.pts=time2.pts)

Матрицы обрабатываются быстрее, чем фреймы данных, поэтому это должно быть быстрее, чем работа с фреймами данных. Теперь анализ:

 results <- lapply(seq(ntimes), function(i) rdist(ref[[i]], pts[[i]]))
distances <- do.call(rbind, results)
colnames(distances) <- c("dista1", "dista2", "distb1", "distb2")

distances Матрица содержит все расстояния. Теперь мы просто объединяем их с вашими df :

 df <- data.frame(df, distances)
options(digits=4)
head(df, 5); cat(". . . . .n"); tail(df, 5)
#   time x y a1.x a1.y a2.x a2.y b1.x b1.y b2.x b2.y dista1 dista2 distb1 distb2
# 1    1 1 1  4.5    6  0.8  3.2  2.5    5  3.8  1.5  6.103  2.209  4.272 2.8443
# 2    1 2 1  4.5    6  0.8  3.2  2.5    5  3.8  1.5  5.590  2.506  4.031 1.8682
# 3    1 3 1  4.5    6  0.8  3.2  2.5    5  3.8  1.5  5.220  3.111  4.031 0.9434
# 4    1 4 1  4.5    6  0.8  3.2  2.5    5  3.8  1.5  5.025  3.883  4.272 0.5385
# 5    1 5 1  4.5    6  0.8  3.2  2.5    5  3.8  1.5  5.025  4.741  4.717 1.3000
# . . . . .
#    time x y a1.x a1.y a2.x a2.y b1.x b1.y b2.x b2.y dista1 dista2 distb1 distb2
# 46    2 1 5    4    5  1.5  3.9  1.4  4.6    6  5.2      3  1.208 0.5657  5.004
# 47    2 2 5    4    5  1.5  3.9  1.4  4.6    6  5.2      2  1.208 0.7211  4.005
# 48    2 3 5    4    5  1.5  3.9  1.4  4.6    6  5.2      1  1.860 1.6492  3.007
# 49    2 4 5    4    5  1.5  3.9  1.4  4.6    6  5.2      0  2.731 2.6306  2.010
# 50    2 5 5    4    5  1.5  3.9  1.4  4.6    6  5.2      1  3.669 3.6222  1.020

Комментарии:

1. Это, безусловно, самый быстрый. Я собирался опубликовать решение в 2,5 раза быстрее, чем у Джонаса, но ваше в 14 раз быстрее моего.

2. Спасибо, что проверили Руи.

3. Мне очень нравится это решение, раньше я не сталкивался с пакетом fields. Но у меня также есть фактор времени. Есть ли способ сказать ей, чтобы она сложила расстояния друг на друга. таким образом, у вас будет df всех расстояний между точкой отсчета и позициями для времени 1, затем ниже этого времени 2 и так далее. отредактировали код вверху

Ответ №2:

Я бы предложил base-R использовать

 referencePointList <- list(a1 = c(4.5,6), a2 = c(0.8,3.2), b1 = c(2.5,5),b2 = c(3.8,1.5))

distanceDfToReferencePoints <- function(x,y,referencePointList) {
  distDf <- setNames(data.frame(do.call("cbind", lapply(referencePointList, function(rp) {
    sqrt((x-rp[1])^2 (y-rp[2])^2)
  }))), paste0("dist_",names(referencePointList)))
  cbind(data.frame(x=x,y=y),distDf)
}

Давайте поместим ваш метод в функцию, скажем

 f0 <- function() {
  df <- as.data.frame(cbind(x,y,a1.x,a1.y,a2.x,a2.y,b1.x,b1.y,b2.x,b2.y))
  df %>% 
    mutate(dista1 = sqrt((x-a1.x)^2   (y-a1.y)^2)) %>% 
    mutate(dista2 = sqrt((x-a2.x)^2   (y-a2.y)^2)) %>% 
    mutate(distb1 = sqrt((x-b1.x)^2   (y-b1.y)^2)) %>% 
    mutate(distb2 = sqrt((x-b2.x)^2   (y-b2.y)^2)) 
}

и сделайте контрольный показатель:

 microbenchmark::microbenchmark(distanceDfToReferencePoints(x,y,referencePointList),f0())

даю на своем автоответчике

 Unit: microseconds
                                                  expr    min      lq     mean  median     uq     max neval
 distanceDfToReferencePoints(x, y, referencePointList)  302.3  340.90  543.882  411.15  479.8 11869.5   100
                                                  f0() 6074.3 6557.85 7497.018 6994.55 7591.5 28291.7   100

Следовательно base-R , решение, по-видимому, намного быстрее. Может быть, это уже помогает вам.

Ответ №3:

подход data.table 

 library(data.table)
setDT(df)
#get names of points
refpoints <- unique( sub( "(^.*)\.[xy]", "\1", names( df[, -c(1,2)] ) ) )
#melt
DT <- melt(df, id.vars = c("x","y"), measure.vars = patterns(x2 = ".*\.x", y2 = ".*\.y"))
#set points' names
setattr(DT$variable, "levels", refpoints )
#calculate distance
DT[, distance := sqrt((x-x2)^2   (y-y2)^2)]
#cast to wide again
dcast(DT, x   y ~ paste0("dist_", variable), value.var = "distance")

#    x y  dist_a1   dist_a2  dist_b1   dist_b2
# 1: 1 1 6.103278 2.2090722 4.272002 2.8442925
# 2: 1 2 5.315073 1.2165525 3.354102 2.8442925
# 3: 1 3 4.609772 0.2828427 2.500000 3.1764760
# 4: 1 4 4.031129 0.8246211 1.802776 3.7536649
# 5: 1 5 3.640055 1.8110770 1.500000 4.4821870
# 6: 2 1 5.590170 2.5059928 4.031129 1.8681542
# 7: 2 2 4.716991 1.6970563 3.041381 1.8681542
# 8: 2 3 3.905125 1.2165525 2.061553 2.3430749
# 9: 2 4 3.201562 1.4422205 1.118034 3.0805844
#10: 2 5 2.692582 2.1633308 0.500000 3.9357337
#11: 3 1 5.220153 3.1112698 4.031129 0.9433981
#12: 3 2 4.272002 2.5059928 3.041381 0.9433981
#13: 3 3 3.354102 2.2090722 2.061553 1.7000000
#14: 3 4 2.500000 2.3409400 1.118034 2.6248809
#15: 3 5 1.802776 2.8425341 0.500000 3.5902646
#16: 4 1 5.024938 3.8832976 4.272002 0.5385165
#17: 4 2 4.031129 3.4176015 3.354102 0.5385165
#18: 4 3 3.041381 3.2062439 2.500000 1.5132746
#19: 4 4 2.061553 3.2984845 1.802776 2.5079872
#20: 4 5 1.118034 3.6715120 1.500000 3.5057096
#21: 5 1 5.024938 4.7413078 4.716991 1.3000000
#22: 5 2 4.031129 4.3680659 3.905125 1.3000000
#23: 5 3 3.041381 4.2047592 3.201562 1.9209373
#24: 5 4 2.061553 4.2755117 2.692582 2.7730849
#25: 5 5 1.118034 4.5694639 2.500000 3.7000000
#    x y  dist_a1   dist_a2  dist_b1   dist_b2

Ответ №4:

Мое предложение состоит в том, чтобы использовать Rfast::dista() и посмотреть время самостоятельно.

Ответ №5:

использование tidyverse:

 df %>%
  rename(x_new = x, y_new = y)%>%
  pivot_longer(3:ncol(df), names_pattern ="(\w )\.(\w )", 
               names_to = c('var', '.value')) %>%
  mutate(value = sqrt((x-x_new)^2   (y-y_new)^2)) %>%
  pivot_wider(c(x_new,y_new),var, values_from = value) 

 x_new y_new    a1    a2    b1    b2
   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
 1     1     1  6.10  2.21  4.27 2.84 
 2     2     1  5.59  2.51  4.03 1.87 
 3     3     1  5.22  3.11  4.03 0.943
 4     4     1  5.02  3.88  4.27 0.539
 5     5     1  5.02  4.74  4.72 1.3  
 6     1     2  5.32  1.22  3.35 2.84 
 7     2     2  4.72  1.70  3.04 1.87 
 8     3     2  4.27  2.51  3.04 0.943
 9     4     2  4.03  3.42  3.35 0.539
10     5     2  4.03  4.37  3.91 1.3  
# ... with 15 more rows

Если бы вы вообще могли извлекать столбцы вручную, вы могли бы даже сделать это быстрее:

 pts <- cbind(x= c(a1.x, a2.x, b1.x, b2.x), y=c(a1.y, a2.y, b1.y, b2.y))
ref <- cbind(x, y)
sqrt(laGP::distance(ref, pts))  
          [,1]      [,2]     [,3]      [,4]
 [1,] 6.103278 2.2090722 4.272002 2.8442925
 [2,] 5.590170 2.5059928 4.031129 1.8681542
 [3,] 5.220153 3.1112698 4.031129 0.9433981
 [4,] 5.024938 3.8832976 4.272002 0.5385165
 [5,] 5.024938 4.7413078 4.716991 1.3000000
 [6,] 5.315073 1.2165525 3.354102 2.8442925
 [7,] 4.716991 1.6970563 3.041381 1.8681542
 [8,] 4.272002 2.5059928 3.041381 0.9433981

или даже:

 apply(pts, 1, function(x)sqrt(rowSums((x-ref)^2)))
          [,1]      [,2]     [,3]      [,4]
 [1,] 6.103278 2.2090722 4.272002 2.8442925
 [2,] 5.315073 1.2165525 3.354102 2.8442925
 [3,] 5.220153 3.1112698 4.031129 0.9433981
 [4,] 4.031129 0.8246211 1.802776 3.7536649
 [5,] 5.024938 4.7413078 4.716991 1.3000000
 [6,] 5.590170 2.5059928 4.031129 1.8681542
 [7,] 4.716991 1.6970563 3.041381 1.8681542
 [8,] 3.905125 1.2165525 2.061553 2.3430749

Если вы используете R >= 4.1

 df |>
  reshape(matrix(3:ncol(df), 2), dir='long') |>
  transform(new_pt = sqrt((x-a1.x)^2   (y - a1.y)^2)) |>
  reshape(v.names =  c('a1.x','a1.y','new_pt'),dir='wide')

Если нет, замените |> трубу на %>%

Метки: