查找多个点之间的最短距离

一种可能性是将识别最近元素(每组一个)的问题表述为混合整数规划。我们可以为是否选择每个点i定义决策变量y_i,以及为是否同时选择点i和j定义x_{ij}(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 

这远非瞬时的——100点7组数据集需要15秒,200点7组数据集需要323秒。不过,它比遍历100点数据集中的所有9200万个7元组或200点数据集中的所有138亿个7元组要快得多。您可以使用类似于Rglpk包中的解算器设置运行时限制,以获得在该限制内获得的最佳解。

所以我想你必须做一个 分支和边界 优化方法。

现在做一些简单的优化:对于每个标签,尝试是否有一些点可以用来代替当前点来改善结果。当你找不到任何进一步的改进时停止。

对于这个初始猜测,计算距离。这将给你一个上限,允许你提前停止搜索。您还可以计算下限,即所有最佳两个标签解的总和。

然后,您可以开始枚举解(可能首先从最小的标签开始),但只要当前解+剩余的下限大于最著名的解(分支和界限),就停止递归。

您也可以尝试对点进行排序,例如按与剩余标签的最小距离排序,以期快速找到更好的边界。

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

我开发了一个简单的算法来快速解决这个问题。第一步是在整个点区域上覆盖网格。第一步是将每组中的每个点分配到其所在的单元或单位正方形。接下来我们转到图的左下角,遍历一个单元格,然后向上遍历一个单元格。这是起始单元格。然后,我们定义一个感兴趣的区域,该区域由该单元及其所有8个相邻单元组成。然后进行测试,以确定每个组的至少一个点是否在这9个单元格区域内。如果是这样,则计算从该区域中表示的每个点到所有其他组中所有其他点的距离。换句话说,该9单元区域中的所有点组合用于获得总距离,其中用于距离计算的成对点从不来自同一组。从这些计算中,将涉及每组单个点的最小距离保存为可能的解。然后,通过右边的一个单元格重复整个过程。当中心单元向右移动时,计算每个9单元区域。这是从右端一个单元格停止的。当第一行完成时,该过程继续进行,向上移动一行,然后从左侧的一个单元格重新开始。因此,在完成顶行时,已考虑每个单元格。解决方案是根据每个9单元区域的所有测试计算出的最小距离。

我们考虑9单元区域而不仅仅是逐个单元的原因是,我们可能会错过位于单元角部的不同组的密集点。

选择正确的单元格或网格大小很重要。如果单元太小,则将找不到可能的解决方案,因为没有任何区域将包含每个组中的至少一个点。如果单元格太大,则每组将有许多点,计算时间将过长。幸运的是,通过反复试验可以很快找到这种最佳单元大小。