找到多個點之間的最短距離

Question

想象一下xy座標的小數據集。這些點由一個名爲indexR的變量組成，共有3組。所有的xy座標都是相同的單位。數據看起來大致像這樣：

# A tibble: 61 x 3 
    indexR  x  y 
    <dbl> <dbl> <dbl> 
1  1 837 924 
2  1 464 661 
3  1 838 132 
4  1 245 882 
5  1 1161 604 
6  1 1185 504 
7  1 853 870 
8  1 1048 859 
9  1 1044 514 
10  1 141 938 
# ... with 51 more rows 

的目標是確定哪些3點，從每個組，是彼此最接近，在最小化所選擇點之間的成對距離的總和的感覺。

我試圖通過考慮歐幾里德距離，如下所示。 （幸得@Mouad_S，在這個線程，和https://gis.stackexchange.com/questions/233373/distance-between-coordinates-in-r）

#dput provided at bottom of this post 
> df$dummy = 1 
> df %>% 
+ full_join(df, c("dummy" = "dummy")) %>% 
+ full_join(df, c("dummy" = "dummy")) %>% 
+ filter(indexR.x != indexR.y & indexR.x != indexR & indexR.y != indexR) %>% 
+ mutate(dist = 
+   ((.$x - .$x.x)^2 + (.$y- .$y.x)^2)^.5 + 
+   ((.$x - .$x.y)^2 + (.$y- .$y.y)^2)^.5 + 
+   ((.$x.x - .$x.y)^2 + (.$y.x- .$y.y)^2)^.5, 
+   dist = round(dist, digits = 0)) %>% 
+ arrange(dist) %>% 
+ filter(dist == min(dist)) 
# A tibble: 6 x 11 
    indexR.x x.x y.x dummy indexR.y x.y y.y indexR  x  y dist 
    <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> 
1  1 638 324  1  2 592 250  3 442 513 664 
2  1 638 324  1  3 442 513  2 592 250 664 
3  2 592 250  1  1 638 324  3 442 513 664 
4  2 592 250  1  3 442 513  1 638 324 664 
5  3 442 513  1  1 638 324  2 592 250 664 
6  3 442 513  1  2 592 250  1 638 324 664 

由此我們可以識別最接近的三個點一起（最小距離隔開;放大在下面的圖）。然而，當擴展這個指數R有4,5 ... n個組時，挑戰就來了。問題在於找到一個更實用或最優化的方法來進行此計算。

structure(list(indexR = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 
2, 2, 2, 2, 2, 3, 3), x = c(836.65, 464.43, 838.12, 244.68, 1160.86, 
1184.52, 853.4, 1047.96, 1044.2, 141.06, 561.01, 1110.74, 123.4, 
1087.24, 827.83, 100.86, 140.07, 306.5, 267.83, 1118.61, 155.04, 
299.52, 543.5, 782.25, 737.1, 1132.14, 659.48, 871.78, 1035.33, 
867.81, 192.94, 1167.8, 1099.59, 1097.3, 1089.78, 1166.59, 703.33, 
671.64, 346.49, 440.89, 126.38, 638.24, 972.32, 1066.8, 775.68, 
591.86, 818.75, 953.63, 1104.98, 1050.47, 722.43, 1022.17, 986.38, 
1133.01, 914.27, 725.15, 1151.52, 786.08, 1024.83, 246.52, 441.53 
), y = c(923.68, 660.97, 131.61, 882.23, 604.09, 504.05, 870.35, 
858.51, 513.5, 937.7, 838.47, 482.69, 473.48, 171.78, 774.99, 
792.46, 251.26, 757.95, 317.71, 401.93, 326.32, 725.89, 98.43, 
414.01, 510.16, 973.61, 445.33, 504.54, 669.87, 598.75, 225.27, 
789.45, 135.31, 935.51, 270.38, 241.19, 595.05, 401.25, 160.98, 
778.86, 192.17, 323.76, 361.08, 444.92, 354, 249.57, 301.64, 
375.75, 440.03, 428.79, 276.5, 408.84, 381.14, 459.14, 370.26, 
304.05, 439.14, 339.91, 435.85, 759.42, 513.37)), class = c("tbl_df", 
"tbl", "data.frame"), row.names = c(NA, -61L), .Names = c("indexR", 
"x", "y")) 

Answer 1

一種可能性是將每個組中最近的元素標識爲混合整數程序。我們可以定義決策變量y_i是否選擇每個點i，以及x_ {ij}是否選擇了點i和j（x_ {ij} = y_iy_j）。我們需要從每個組中選擇一個元素。

實際上，您可以使用lpSolve包（或其他R優化包之一）實現此混合整數程序。

opt.closest <- function(df) { 
    # Compute every pair of indices 
    library(dplyr) 
    pairs <- as.data.frame(t(combn(nrow(df), 2))) %>% 
    mutate(G1=df$indexR[V1], G2=df$indexR[V2]) %>% 
    filter(G1 != G2) %>% 
    mutate(dist = sqrt((df$x[V1]-df$x[V2])^2+(df$y[V1]-df$y[V2])^2)) 

    # Compute a few convenience values 
    n <- nrow(df) 
    nP <- nrow(pairs) 
    groups <- sort(unique(df$indexR)) 
    nG <- length(groups) 
    gpairs <- combn(groups, 2) 
    nGP <- ncol(gpairs) 

    # Solve the optimization problem 
    obj <- c(pairs$dist, rep(0, n)) 
    constr <- rbind(cbind(diag(nP), -outer(pairs$V1, seq_len(n), "==")), 
        cbind(diag(nP), -outer(pairs$V2, seq_len(n), "==")), 
        cbind(diag(nP), -outer(pairs$V1, seq_len(n), "==") - outer(pairs$V2, seq_len(n), "==")), 
        cbind(matrix(0, nG, nP), outer(groups, df$indexR, "==")), 
        cbind((outer(gpairs[1,], pairs$G1, "==") & 
         outer(gpairs[2,], pairs$G2, "==")) | 
         (outer(gpairs[2,], pairs$G1, "==") & 
         outer(gpairs[1,], pairs$G2, "==")), matrix(0, nGP, n))) 
    dir <- rep(c("<=", ">=", "="), c(2*nP, nP, nG+nGP)) 
    rhs <- rep(c(0, -1, 1), c(2*nP, nP, nG+nGP)) 
    library(lpSolve) 
    mod <- lp("min", obj, constr, dir, rhs, all.bin=TRUE) 
    which(tail(mod$solution, n) == 1) 
} 

這可以計算最接近的3個點，從每個集羣，在您的示例數據集：

df[opt.closest(df),] 
# A tibble: 3 x 3 
# indexR  x  y 
# <dbl> <dbl> <dbl> 
# 1  1 638.24 323.76 
# 2  2 591.86 249.57 
# 3  3 441.53 513.37 

它也可以計算與更多的積分和組數據集最佳的解決方案。以下是運行時用於與每個7組，100和200點的數據集：

make.dataset <- function(n, nG) { 
    set.seed(144) 
    data.frame(indexR = sample(seq_len(nG), n, replace=T), x = rnorm(n), y=rnorm(n)) 
} 
df100 <- make.dataset(100, 7) 
system.time(opt.closest(df100)) 
# user system elapsed 
# 11.536 2.656 15.407 
df200 <- make.dataset(200, 7) 
system.time(opt.closest(df200)) 
# user system elapsed 
# 187.363 86.454 323.167 

這遠遠瞬時 - 它需要15秒的100點，7組數據集和323秒鐘的200點，7組數據集。儘管如此，它比遍歷100點數據集中的所有9200萬個7元組或130點數據集中的所有138億7元組要快得多。您可以使用Rglpk軟件包中的求解器來設置運行時限制，以獲得在此限制內獲得的最佳解決方案。

Answer 2

您可以使用交叉連接一起把所有的點的組合，計算所有的三個點之間的總距離，然後取最小的那個。

df$id <- row.names(df) # to create ID's for the points 

df2 <- merge(df, df, by = NULL) # the first cross join 

df3 <- merge(df2, df, by = NULL) # the second cross join 



# eliminating rows where the points are of the same indexR 

df3 <- df3[df3$indexR.x != df3$indexR.y & df3$indexR.x != df3$indexR 
      & df3$indexR.y != df3$indexR,] 


## calculating the total distance 

df3$total_distance <- ((df3$x - df3$x.x)^2 + (df3$y- df3$y.x)^2)^.5 + 
    ((df3$x - df3$x.y)^2 + (df3$y- df3$y.y)^2)^.5 + 
    ((df3$x.x - df3$x.y)^2 + (df3$y.x- df3$y.y)^2)^.5 

## minimum distance 

df3[which.min(df3$total_distance),] 

indexR.x x.x y.x id.x indexR.y x.y y.y id.y indexR  x  y id total_distance 
155367  3 441.53 513.37 61  2 591.86 249.57 46  1 638.24 323.76 42  664.3373 

Answer 3

您不能列舉所有可能的解決方案，而且我也沒有看到任何明顯的捷徑。

所以，我想你將不得不做一個分支和約束優化的方法。

首先猜測一個相當好的解決方案。就像最近的兩個不同標籤的點一樣。然後添加最近的不同標籤，直到您覆蓋所有標籤。

現在做一些簡單的優化：對於每個標籤，嘗試是否有某個點可以用來代替當前點來改善結果。當你找不到任何進一步的改進時停下來。

對於這個初始猜測，計算距離。這會給你一個上限，這可以讓你儘早停止搜索。您還可以計算一個下限，即所有最佳雙標籤解決方案的總和。

現在您可以嘗試刪除點，其中每個標籤的最近鄰居+所有其他標籤的下界已經比您的初始解決方案更差。這將有希望消除很多點。

然後，您可以開始枚舉解決方案（可能首先從最小的標籤開始），但是隻要當前解決方案+剩餘下限大於您最知名的解決方案（分支和界限）就停止遞歸。

你也可以嘗試排序點，例如通過與其餘標籤的最小距離，希望快速找到更好的邊界。

我肯定不會選擇R實現這個...