查找多个点之间的最短距离
Find shortest distance between multiple points
想象一个 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
由此我们可以找出距离最近的三个点（最小距离；下图放大）。然而，当扩展它使得 indexR 有 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"))
您可以使用交叉连接获得所有点组合，计算所有三个点之间的总距离，然后取最小值。
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
一种可能性是将识别最接近的元素（每组一个）的问题表述为混合整数程序。我们可以为每个点 i 是否被 selected 定义决策变量 y_i，以及为点 i 和 j 是否都被 selected 定义 x_{ij}（x_{ij} = y_iy_j)。我们需要 select 每组中的一个元素。
实际上，您可以使用 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
这远非瞬时——100 点、7 组数据集需要 15 秒，200 点、7 组数据集需要 323 秒。不过，它比遍历 100 点数据集中的所有 9200 万个 7 元组或 200 点数据集中的所有 138 亿个 7 元组要快得多。您可以使用 Rglpk 包中的求解器设置运行时限制，以获得在该限制内获得的最佳解决方案。
您无力列举所有可能的解决方案，而且我看不到任何明显的捷径。
所以我猜你必须做一个分支定界优化方法。
先猜一个比较好的解法。就像具有不同标签的最近的两个点。然后添加最近的不同标签，直到覆盖所有标签。
现在做一些微不足道的优化：对于每个标签，尝试是否有一些点可以用来代替当前点来改善结果。当你找不到任何进一步的改进时就停止。
对于这个初始猜测，计算距离。这会给你一个上限，让你提前停止搜索。您还可以计算下限，即所有最佳双标签解决方案的总和。
现在你可以尝试移除点，其中每个标签的最近邻居+所有其他标签的下界已经比你的初始解决方案差。这将有望消除很多点。
然后您可以开始枚举解决方案（可能首先从最小的标签开始），但是只要当前解决方案 + 剩余的下限大于您最知名的解决方案（分支定界法）就停止递归。
您也可以尝试对点进行排序，例如通过与剩余标签的最小距离，希望能快速找到更好的界限。
我当然不会选择 R 来实现这个...
我开发了一个简单的算法来快速解决这个问题。第一步是在整个点区域上覆盖一个网格。第一步是将每个组中的每个点分配给它所在的单元格或单元格。接下来我们转到图表的左下角，越过一个单元格并向上移动一个单元格。这是起始单元格。然后我们定义一个由该单元格及其所有 8 个邻居组成的感兴趣区域。然后进行测试以确定每个组中的至少一个点是否在这 9 个单元格区域内。如果是这样，则计算从每个点组到该区域中表示的每个点到所有其他组的所有其他点的距离。换句话说，这个 9 单元格区域中的所有点组合都用于获得总距离，其中用于距离计算的成对点永远不会来自同一组。从这些计算中，涉及每个组中单个点的最小距离的计算被保存为可能的解决方案。然后通过向右移动一个单元格来重复整个过程。当中央单元格向右移动时，将计算每个 9 单元格区域。这是从右端停止一个单元格。当第一行完成时，该过程通过上升一行并在左侧再次开始但再次开始一个单元格来继续。因此，当顶行完成时，每个单元格都已被考虑。解决方案将是根据为每个 9 单元格区域完成的所有测试计算出的最小距离。
我们考虑 9 单元格区域而不是逐个单元格的原因是我们可能会错过来自位于单元格角落的不同组的紧密间隔的点。
选择正确的单元格或网格大小很重要。如果单元格太小，则无法找到可能的解决方案，因为 none 个区域将包含每个组中的至少一个点。如果单元格太大，那么每组的点就会很多，计算时间也会过长。幸运的是，可以通过反复试验快速找到最佳单元格大小。
我已经运行多次使用不同数量的组和组中的点数来使用此算法。对于所有组中随机分散的点，我发现 15 x 15 的网格大小适用于 10 组 - 400 点（每组 40 点）的情况。该示例 运行 不到一秒钟。
想象一个 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
由此我们可以找出距离最近的三个点（最小距离；下图放大）。然而，当扩展它使得 indexR 有 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"))
您可以使用交叉连接获得所有点组合，计算所有三个点之间的总距离，然后取最小值。
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
一种可能性是将识别最接近的元素（每组一个）的问题表述为混合整数程序。我们可以为每个点 i 是否被 selected 定义决策变量 y_i，以及为点 i 和 j 是否都被 selected 定义 x_{ij}（x_{ij} = y_iy_j)。我们需要 select 每组中的一个元素。
实际上，您可以使用 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
这远非瞬时——100 点、7 组数据集需要 15 秒，200 点、7 组数据集需要 323 秒。不过，它比遍历 100 点数据集中的所有 9200 万个 7 元组或 200 点数据集中的所有 138 亿个 7 元组要快得多。您可以使用 Rglpk 包中的求解器设置运行时限制，以获得在该限制内获得的最佳解决方案。
您无力列举所有可能的解决方案，而且我看不到任何明显的捷径。
所以我猜你必须做一个分支定界优化方法。
先猜一个比较好的解法。就像具有不同标签的最近的两个点。然后添加最近的不同标签，直到覆盖所有标签。
现在做一些微不足道的优化：对于每个标签，尝试是否有一些点可以用来代替当前点来改善结果。当你找不到任何进一步的改进时就停止。
对于这个初始猜测，计算距离。这会给你一个上限，让你提前停止搜索。您还可以计算下限，即所有最佳双标签解决方案的总和。
现在你可以尝试移除点，其中每个标签的最近邻居+所有其他标签的下界已经比你的初始解决方案差。这将有望消除很多点。
然后您可以开始枚举解决方案（可能首先从最小的标签开始），但是只要当前解决方案 + 剩余的下限大于您最知名的解决方案（分支定界法）就停止递归。
您也可以尝试对点进行排序，例如通过与剩余标签的最小距离，希望能快速找到更好的界限。
我当然不会选择 R 来实现这个...
我开发了一个简单的算法来快速解决这个问题。第一步是在整个点区域上覆盖一个网格。第一步是将每个组中的每个点分配给它所在的单元格或单元格。接下来我们转到图表的左下角，越过一个单元格并向上移动一个单元格。这是起始单元格。然后我们定义一个由该单元格及其所有 8 个邻居组成的感兴趣区域。然后进行测试以确定每个组中的至少一个点是否在这 9 个单元格区域内。如果是这样，则计算从每个点组到该区域中表示的每个点到所有其他组的所有其他点的距离。换句话说，这个 9 单元格区域中的所有点组合都用于获得总距离，其中用于距离计算的成对点永远不会来自同一组。从这些计算中，涉及每个组中单个点的最小距离的计算被保存为可能的解决方案。然后通过向右移动一个单元格来重复整个过程。当中央单元格向右移动时，将计算每个 9 单元格区域。这是从右端停止一个单元格。当第一行完成时，该过程通过上升一行并在左侧再次开始但再次开始一个单元格来继续。因此，当顶行完成时，每个单元格都已被考虑。解决方案将是根据为每个 9 单元格区域完成的所有测试计算出的最小距离。
我们考虑 9 单元格区域而不是逐个单元格的原因是我们可能会错过来自位于单元格角落的不同组的紧密间隔的点。
选择正确的单元格或网格大小很重要。如果单元格太小，则无法找到可能的解决方案，因为 none 个区域将包含每个组中的至少一个点。如果单元格太大，那么每组的点就会很多，计算时间也会过长。幸运的是，可以通过反复试验快速找到最佳单元格大小。
我已经运行多次使用不同数量的组和组中的点数来使用此算法。对于所有组中随机分散的点，我发现 15 x 15 的网格大小适用于 10 组 - 400 点（每组 40 点）的情况。该示例 运行 不到一秒钟。