9

考虑一个矩阵，每条线指定一个二维区域，另一个矩阵指定平面中的点：

xmin <- c(3, 14, 25, 61)
xmax <- c(5, 18, 27, 65)
ymin <- c(33, 12, 83, 2)
ymax <- c(35, 16, 90, 6)
regions <- cbind(xmin, xmax, ymin, ymax)

x <- c(7, 26, 4, 16)
y <- c(4, 85, 30, 13)
points <- cbind(x, y)

regions获取包含 中每个点的索引的最快方法是什么points？

我想要实现的一个例子是：

apply(points, 1, function(x){
    which(regions[,'xmin'] < x[1] & regions[,'xmax'] > x[1] & regions[,'ymin'] < x[2] & regions[,'ymax'] > x[2])
})

regions但是随着两者中的行数points接近 1E5 这变得相当慢，我正在寻找一种适当的矢量化方法......

提前致谢...

最好的托马斯

编辑：

对于任何感兴趣的人，我使用 Rcpp 在 C++ 中创建了一个函数，它提供了大约 50 倍的性能改进。我不精通C++，所以它可能会做得更好......

cppFunction('
    IntegerVector findInRegion(NumericVector x, NumericVector y, NumericVector xmin, NumericVector xmax, NumericVector ymin, NumericVector ymax){
        int pointSize = x.size();
        int regionSize = xmin.size();
        IntegerVector ans(pointSize);
        for(int i = 0; i < pointSize; i++){
            ans[i] = NA_INTEGER;
        }

        for(int i = 0; i < pointSize; i++){
            for(int j = 0; j < regionSize; j++){
                if(x[i] > xmin[j]){
                    if(x[i] < xmax[j]){
                        if(y[i] > ymin[j]){
                            if(y[i] < ymax[j]){
                                ans[i] = j+1;
                            };
                        };
                    };
                };
            };
        };
        return ans;
    }
')

findRegion <- function(points, regions){
    if(!all(c('x', 'y') %in% colnames(points))){
        stop('points must contain columns named \'x\' and \'y\'')
    }
    if(!all(c('xmin', 'xmax', 'ymin', 'ymax') %in% colnames(regions))){
        stop('regions must contain columns named \'xmin\', \'xmax\', \'ymin\' and \'ymax\'')
    }
    findInRegion(points[, 'x'], points[,'y'], regions[, 'xmin'], regions[, 'xmax'], regions[, 'ymin'], regions[, 'ymax'])
}

这个函数的一个缺点是它假设一个点只能属于一个区域......

4

2 回答 2

4

这是一个非常有趣的问题。我做了一些初步测试，这似乎可能会更快，但我真的不知道它的扩展性如何。如果您可以测试您的真实数据并报告一些时间，我会很感兴趣：

#  Are X coords greater than xmin
lx <- outer( points[,1] , regions[,1] , ">" )

#  Are X coords less than xmax
hx <- outer( points[,1] , regions[,2] , "<" )

#  Ditto for Y coords
ly <- outer( points[,2] , regions[,3] , ">" )
hy <- outer( points[,2] , regions[,4] , "<" )

#  These matrices for X and Y points have 1 if coords is in range, 0 otherwise
inx <- lx * hx
iny <- ly * hy

#  The final result matrix has 1 if both X and Y coords are in range and 0 if not
#  Rows are points, columns are regions
res <- inx * iny

对于 100000 个点和 100000 个区域的数据，除非您有大量RAM，否则这种方法将不起作用。但是，我认为如果您将区域数量分成大约 1000 个的块，它会非常有用。在我的桌面上，100,000 个点和 1,000 个区域需要 5 秒：

Unit: seconds
        expr      min      lq  median       uq      max neval
 eval(simon) 4.528942 4.55258 4.59848 4.607572 4.671511     5

作为我在您的方法和这个方法之间看到的时间差异幅度的粗略指南apply，有 10,000 个点和 1,000 个区域（基于 5 次运行）：

Unit: milliseconds
        expr       min        lq    median        uq       max neval
 eval(simon)  394.7165  402.0919  403.0491  404.6943  428.7077     5
    eval(OP) 1359.5889 1364.6308 1372.4980 1383.1327 1491.4628     5

并且有 100,000 个点和 1,000 个区域（基于一次运行）：

Unit: seconds
        expr       min        lq    median        uq       max neval
 eval(simon)  4.352857  4.352857  4.352857  4.352857  4.352857     1
    eval(OP) 14.027390 14.027390 14.027390 14.027390 14.027390     1

这是我用来生成示例数据和运行基准测试的代码：

set.seed(4862)
xmin <- sample(25,1000,repl=T)
xmax <- xmin + sample(15,100,repl=T)
ymin <- sample(25,1000,repl=T)
ymax <- ymin + sample(15,1000,repl=T)
regions <- cbind(xmin, xmax, ymin, ymax)

x <- sample(25,100000,repl=T)
y <- sample(25,100000,repl=T)
points <- cbind(x, y)


OP <- quote({ res <- apply(points, 1, function(x){
    which(regions[,'xmin'] < x[1] & regions[,'xmax'] > x[1] & regions[,'ymin'] < x[2] & regions[,'ymax'] > x[2])
}) })


simon <- quote({
lx <- outer( points[,1] , regions[,1] , ">" )
hx <- outer( points[,1] , regions[,2] , "<" )
ly <- outer( points[,2] , regions[,3] , ">" )
hy <- outer( points[,2] , regions[,4] , "<" )
inx <- lx * hx
iny <- ly * hy
res <- inx * iny })

require(microbenchmark)
microbenchmark( eval(simon) , eval(OP) , times = 1L )

我建议分块进行。HTH。

于 2013-06-10T08:56:35.797 回答
4

这是另一种解决方案，使用带有 SQLite 的R-tree索引（一种旨在存储边界框的数据库索引）。结果证明它比 Simon 的（7 秒）稍慢，可能是因为数据被复制到磁盘。

# Sample data: data.frames, rather than matrices
regions <- data.frame(id=1:length(xmin), xmin, xmax, ymin, ymax)
points  <- data.frame(x, y)

library(RSQLite)
con <- dbConnect("SQLite", dbname = "/tmp/a.sqlite") 
dbGetQuery( con, "CREATE VIRTUAL TABLE regions USING rtree (id, xmin, xmax, ymin, ymax)" )
dbWriteTable( con, "regions", regions, row.names = FALSE, append = TRUE )
dbWriteTable( con, "points",  points, row.names = TRUE )
res <- dbGetQuery( con, "
  SELECT points.row_names, regions.id
  FROM   points, regions
  WHERE  xmin <= x AND x <= xmax
  AND    ymin <= y AND y <= ymax
" )
于 2013-06-10T15:00:20.817 回答