Как я могу ускорить этот двойной цикл путем векторизации в R? - PullRequest
Новая
       34

Как я могу ускорить этот двойной цикл путем векторизации в R?

0 голосов
/ 07 июня 2019

У меня есть две таблицы данных. Один с информацией о домах, включая две колонки с координатами (один для Востока и один для Норда). И один с информацией о станциях, включая две колонки с координатами. Координаты в швейцарском формате LV95. 

    > head(d.small)
         id   GKODE   GKODN
1:      100 2608850 1262583
2: 10000012 2612952 1267232
3: 11776807 2612952 1267232
4: 10000061 2623360 1247413
5: 10000223 2622938 1259411
6:  9997602 2609194 1262383

> head(haltestelle.small)
   y_Koord_Ost x_Koord_Nord
1:     2633061      1257737
2:     2630544      1252831
3:     2628497      1256841
4:     2629649      1255767
5:     2619156      1264531
6:     2619741      1247012

Теперь мне нужно расстояние до ближайшей станции.

Мой код ниже вычисляет расстояние между одним домом и всеми станциями, берет минимальное расстояние и добавляет соответствующий индекс.

К сожалению, мой код слишком медленный. Как мне векторизовать мой цикл? 

dist.oev <- data.table(dist.oev=rep(1, nrow(d.small)), dist.oev.index=rep(1, nrow(d.small)))

for (i in 1:nrow(d.small)) {
  cat(i, " ")
  for ( j in 1:nrow(haltestelle.small)) {
    diff.ost <- d.small[i, .(GKODE)] - haltestelle.small[j, .(y_Koord_Ost)]
    diff.nord <- d.small[i, .(GKODN)] - haltestelle.small[j, .(x_Koord_Nord)]
    dist.oev[i,1] <- min(sqrt(diff.ost^2 + diff.nord^2))
    dist.oev[i,2] <- which.min(sqrt(diff.ost^2 + diff.nord^2))
  }
}

Ответы [ 3 ]

1 голос
/ 07 июня 2019

Это может быть то, что вы ищете: 

d.small[,
        mindist := sqrt(min((GKODE - haltestelle.small[["y_Koord_Ost"]])^2 +
                            (GKODN - haltestelle.small[["x_Koord_Nord"]])^2)),
        by = id]

#          id   GKODE   GKODN   mindist
# 1:      100 2608850 1262583 10488.486
# 2: 10000012 2612952 1267232  6766.463
# 3: 11776807 2612952 1267232  6766.463
# 4: 10000061 2623360 1247413  3641.148
# 5: 10000223 2622938 1259411  6124.327
# 6:  9997602 2609194 1262383 10190.944

Данные (в воспроизводимом формате): 

d.small <- fread("id   GKODE   GKODN
     100 2608850 1262583
10000012 2612952 1267232
11776807 2612952 1267232
10000061 2623360 1247413
10000223 2622938 1259411
 9997602 2609194 1262383")

haltestelle.small <- fread("y_Koord_Ost x_Koord_Nord
2633061      1257737
2630544      1252831
2628497      1256841
2629649      1255767
2619156      1264531
2619741      1247012")

С небольшими правками ваш код не содержит ошибок: 

for (i in 1:nrow(d.small)) {
  diff.ost  <- d.small[i, GKODE] - haltestelle.small[, y_Koord_Ost]
  diff.nord <- d.small[i, GKODN] - haltestelle.small[, x_Koord_Nord]
  dist.oev[i,1] <- sqrt(min(diff.ost^2 + diff.nord^2)) # take sqrt outside min for efficiency
  dist.oev[i,2] <- which.min(diff.ost^2 + diff.nord^2) # sqrt unnecessary (monotonic transformation)
}

Получить индексы также можно примерно так: 

d.small[,
        c("mindist", "mindist_index") := {
          dist = (GKODE - haltestelle.small[["y_Koord_Ost"]])^2 + (GKODN - haltestelle.small[["x_Koord_Nord"]])^2
          .(sqrt(min(dist)), which.min(dist))
        },
        by = id]
1 голос
/ 07 июня 2019

Вот решение, которое должно быть относительно быстрым, если оно умещается в вашей оперативной памяти.Я привожу воспроизводимый пример, хотя это не ваши точные данные.h содержит координаты вашей станции, а d содержит координаты вашего дома: 

h <- data.frame(x=rnorm(1000),y=rnorm(1000))
d <- data.frame(x=rnorm(500),y=rnorm(500))
xdiff <- sapply(h$x,function(x1)sapply(d$x,function(x2)(x1-x2)^2))
ydiff <- sapply(h$y,function(y1)sapply(d$y,function(y2)(y1-y2)^2))
dist2 <-  xdiff+ydiff
closest <- apply(dist2,1,which.min)
min.dist <- sqrt(dist2[cbind(seq(500),closest)])

Код сначала вычисляет все квадратные различия по x, а затем все квадратные различия по y и добавляет их для получения евклидовых квадратоврасстояние в матрице.

0 голосов
/ 08 июня 2019

Решение, использующее пакет sf:

(я не уверен, правильно ли настроены широта и долгота, поскольку я не знаком с форматированием LV95.)

Пример данных 

# Sample data
library(data.table)
d.small <- fread("id   GKODE   GKODN
100 2608850 1262583
10000012 2612952 1267232
11776807 2612952 1267232
10000061 2623360 1247413
10000223 2622938 1259411
9997602 2609194 1262383")

haltestelle.small <- fread("y_Koord_Ost x_Koord_Nord
2633061      1257737
2630544      1252831
2628497      1256841
2629649      1255767
2619156      1264531
2619741      1247012")

Код 

# Create spatial objects
library(sf)
d.sf <- d.small %>% sf::st_as_sf( coords = c("GKODE", "GKODN"), crs = 2056)     
haltestelle.sf <- haltestelle.small %>% sf::st_as_sf( coords = c("y_Koord_Ost", "x_Koord_Nord"), crs = 2056) 

# Calculate nearest haltestelle for each d.small
d.sf %>% 
  dplyr::group_by( id ) %>%
  dplyr::mutate( np = sf::st_nearest_feature( geometry, haltestelle.sf ),
                 dist_np = as.numeric( sf::st_distance( geometry, d.sf[np,] ) ) )

Выход 

        id    np dist_np          geometry
     <int> <int>   <dbl>       <POINT [m]>
1      100     5  14441. (2608850 1262583)
2 10000012     5  12684. (2612952 1267232)
3 11776807     5  12684. (2612952 1267232)
4 10000061     6  20610. (2623360 1247413)
5 10000223     3  12684. (2622938 1259411)
6  9997602     5  14062. (2609194 1262383)
...