【发布时间】：2016-01-22 17:53:27
【问题描述】：
我一直试图在 Internet 上找到解决方案，但未能解决我的具体问题。 我有一个矩阵 D 告诉我一组矩形（列）的最大和最小 x 和 y 坐标。现在我想知道它们中的一些相交的频率。我想评估其中两个是否相交的方法需要遍历所有可能的列组合（也“向后”）。现在对我来说第一个直观的解决方案是嵌套一个 for 循环，通过列进入另一个执行相同操作的循环。这是一个例子：
n <- 5
D <- matrix(c(1,2,1,2,1,3,1,3,2,4,2,4,3,5,1,2,3,4,2,4),nrow=4,ncol=n)
E <- mat.or.vec(nr=n, nc=n)
for (i in 1:n){
for (j in 1:n){
if (i != j &&
(D[1,i] <= D[1,j] & D[1,j] < D[2,i]) &&
((D[3,i] <= D[3,j] & D[3,j] < D[4,i]) | (D[3,i] < D[4,j] & D[4,j] <= D[4,i]))) {
E[i,j] <- 1
}
}
}
[,1] [,2] [,3] [,4] [,5]
[1,] 0 1 0 0 0
[2,] 1 0 1 0 0
[3,] 0 0 0 0 1
[4,] 0 0 0 0 0
[5,] 0 0 0 0 0
我一直在尝试各种替代解决方案，包括使用foreach、outer、apply 的不同版本以及它们的组合。我发现很多版本都给了我相同的正确结果，但是由于某种原因，双 for 循环的执行速度比任何其他方法都快。
问题是，矩阵中有大约 4000 列，我必须运行整个结果函数大约 50 次。所以循环花费的时间仍然太长。谁能告诉我为什么其他解决方案的性能更差（因为我认为循环是最糟糕的事情）并且可能会建议一些可以快速处理大量数据的方法？
编辑： 我尝试了使用 918 个矩形建议的两种解决方案并对其进行了基准测试。
D1 <- t(D)
D1 <- D1[,c(1,3,2,3,2,4,1,4)]
ID <- paste0('rect',seq_len(nrow(D1)))
polys <- SpatialPolygons(mapply(function(poly, id) {
xy <- matrix(poly, ncol=2, byrow=TRUE)
Polygons(list(Polygon(xy)), ID=id)
}, split(D1, row(D1)), ID))
polys.df <- SpatialPolygonsDataFrame(polys,data.frame(id=ID,row.names=ID))
res <- gOverlaps(polys.df, byid = TRUE) |
gContains(polys.df, byid = TRUE) |
gWithin(polys.df, byid = TRUE)
diag(res) <-0
intersecting_rects <- function(edge1, edge2){
#* index 1 is left edge, index 2 is right edge
#* index 3 is lower edge, index 4 is upper edge
(edge1[1] < edge2[2] && edge1[2] > edge2[1] &&
edge1[3] < edge2[4] && edge1[4] > edge2[3])
}
microbenchmark(
poly = {
res <- gOverlaps(polys.df, byid = TRUE) |
gContains(polys.df, byid = TRUE) |
gWithin(polys.df, byid = TRUE) |
gCovers(polys.df, byid = TRUE)
diag(res) <- 0
},
full_loop = {
E <- mat.or.vec(nr=n, nc=n)
for (i in 1:n){
for (j in 1:n){
E[i, j] <- as.numeric(intersecting_rects(D[, i], D[, j]))
}
}
},
partial_loop = {
E <- mat.or.vec(nr=n, nc=n)
for (i in 1:n){
for (j in (i:n)[-1]){
E[i, j] <- as.numeric(intersecting_rects(D[, i], D[, j]))
}
}
E[lower.tri(E)] <- t(E)[lower.tri(E)]
}
)
Unit: seconds
expr min lq mean median uq max neval
poly 5.280785 5.985167 6.277424 6.288105 6.567685 7.141447 100
full_loop 8.328631 9.700908 10.091188 10.112840 10.490162 11.450817 100
partial_loop 4.335070 4.921649 5.203672 5.188610 5.503550 6.000614 100
我们看到，部分循环已经比完整循环有了显着改进，但仅比空间对象的解决方案略好。因此，我选择了 Rolands 解决方案。再次感谢！
【问题讨论】：
-
请提供一个可重现的最小示例。 stackoverflow.com/questions/5963269/…
-
如果您的矩形数量相对较少，我会考虑为每个矩形分配一个位域，然后用它们“绘制”（或）矩阵。
-
嗯，我有大约 4000 个，所以这是一个相对较小的数字吗？
-
如果我错了，请纠正我，但
E应该是一个对称矩阵，不是吗？也就是说，既然E[2, 3] == 1，不应该E[3,2] == 1吗？换句话说，当矩形 2 与矩形 3 相交时，矩形 3 不也与矩形 2 相交吗？我问的原因是，如果这是真的，你不需要遍历矩阵的整个维度。你只需要做上矩阵或三角矩阵，它可以把你的循环减半。 -
另外，
for循环的名声不好。如果使用得当，它真的不是那么糟糕的工具。for循环在必须重新评估它们正在创建的对象类型或用于制作不断扩展的对象时往往表现最差。您已经预先定义了E的尺寸并且没有更改其大小或类型这一事实应该可以最大限度地减少与for循环相关的问题。我认为你必须非常有创意才能让某些东西运行得更快。
标签： r