r - 条件矩阵邻接计算

Question

我有一个矩阵，例如:

set.seed(1)
m = matrix(rep(NA,100), nrow=10)
m[sample(1:100,10)] = 1
m
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
 [1,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA
 [2,]   NA   NA   NA   NA   NA   NA    1   NA   NA    NA
 [3,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA
 [4,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA
 [5,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA
 [6,]    1   NA   NA   NA   NA   NA   NA   NA    1    NA
 [7,]   NA   NA    1    1   NA    1   NA   NA   NA     1
 [8,]   NA   NA   NA   NA   NA    1   NA   NA   NA    NA
 [9,]   NA   NA   NA   NA   NA   NA   NA   NA    1    NA
[10,]   NA    1   NA   NA   NA   NA   NA   NA   NA    NA

我想将下一个(相邻)的所有 NA 值转换为非 NA 值为零。没有一些可怕的逐行和逐列循环算法，是否有任何 swishy 矩阵方式来实现这一点？

注意。我已经重新设计了这个例子以减少歧义。我需要非 NA 值的上方、下方、左侧和右侧的所有 NA 值都变为零!

Answer 1

m[is.na(m) & !(cbind(is.na(m[,-1L]),T) & cbind(T,is.na(m[,-ncol(m)])) & rbind(is.na(m[-1L,]),T) & rbind(T,is.na(m[-nrow(m),])))] <- 0;
m;
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
##  [1,]   NA   NA   NA   NA   NA   NA    0   NA   NA    NA
##  [2,]   NA   NA   NA   NA   NA    0    1    0   NA    NA
##  [3,]   NA   NA   NA   NA   NA   NA    0   NA   NA    NA
##  [4,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA
##  [5,]    0   NA   NA   NA   NA   NA   NA   NA    0    NA
##  [6,]    1    0    0    0   NA    0   NA    0    1     0
##  [7,]    0    0    1    1    0    1    0   NA    0     1
##  [8,]   NA   NA    0    0    0    1    0   NA    0     0
##  [9,]   NA    0   NA   NA   NA    0   NA    0    1     0
## [10,]    0    1    0   NA   NA   NA   NA   NA    0    NA

该解决方案的工作原理如下。

我们用 TRUE 构造一个逻辑索引矩阵其中一个元素是 NA 和 与至少一个非NA元素相邻(上方、下方、左侧或右侧)。然后我们可以下标m使用逻辑索引矩阵并分配所需的替换值。

逻辑合取的 LHS 很容易；简直是is.na(m) .

逻辑合取的 RHS 是最棘手的部分。我们需要执行 4 次测试，每个邻接方向一个。一般算法是:

1:索引邻接方向的奇异索引，相对于该邻接方向不与任何其他索引相邻。例如，对于“右方向”，我们从最左边的列开始索引，因为它不在任何其他索引的右侧。换句话说，没有最左边的列在右边，所以我们可以忽略它(并且必须删除它)以进行“右方向”计算。

2:使用 is.na() 测试 NA 的子矩阵.

3:然后我们必须绑定(bind)(cbind() 用于水平邻接方向，rbind() 用于垂直方向)TRUE在结果逻辑子矩阵的另一侧(即与步骤 1 中删除的索引相反)。这有效地导致邻接方向上的最后一个索引在其邻接方向上始终具有(伪)NA，因此它永远不会由于该邻接方向而被替换。

4:逻辑 和 4个测试。结果将是一个带有 TRUE 的逻辑矩阵。对于在每个相邻单元格中都有 NA 的元素。

5:否定步骤 4 的结果。这将产生一个带有 TRUE 的逻辑矩阵对于在任何相邻单元格中至少有一个非 NA 的元素。

请注意，还有另一种方法可以做到这一点，这可能更直观一些。我们可以编写 4 个测试中的每一个来测试非 NA，而不是 NA，然后是逻辑 或 他们在一起。这也需要绑定(bind) FALSE而不是 TRUE对于最后一个索引。它看起来像这样:

m[is.na(m) & (cbind(!is.na(m[,-1L]),F) | cbind(F,!is.na(m[,-ncol(m)])) | rbind(!is.na(m[-1L,]),F) | rbind(F,!is.na(m[-nrow(m),])))] <- 0;
m;
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
##  [1,]   NA   NA   NA   NA   NA   NA    0   NA   NA    NA
##  [2,]   NA   NA   NA   NA   NA    0    1    0   NA    NA
##  [3,]   NA   NA   NA   NA   NA   NA    0   NA   NA    NA
##  [4,]   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA
##  [5,]    0   NA   NA   NA   NA   NA   NA   NA    0    NA
##  [6,]    1    0    0    0   NA    0   NA    0    1     0
##  [7,]    0    0    1    1    0    1    0   NA    0     1
##  [8,]   NA   NA    0    0    0    1    0   NA    0     0
##  [9,]   NA    0   NA   NA   NA    0   NA    0    1     0
## [10,]    0    1    0   NA   NA   NA   NA   NA    0    NA

第一种方法更可取，因为它只需要一个否定，而第二种方法需要 4 个否定。

基准测试

library(raster);
library(microbenchmark);

bgoldst1 <- function(m) { m[is.na(m) & !(cbind(is.na(m[,-1L]),T) & cbind(T,is.na(m[,-ncol(m)])) & rbind(is.na(m[-1L,]),T) & rbind(T,is.na(m[-nrow(m),])))] <- 0; m; };
bgoldst2 <- function(m) { m[is.na(m) & (cbind(!is.na(m[,-1L]),F) | cbind(F,!is.na(m[,-ncol(m)])) | rbind(!is.na(m[-1L,]),F) | rbind(F,!is.na(m[-nrow(m),])))] <- 0; m; };
geotheory <- function(m) { r <- raster(m,crs='+init=epsg:27700'); extent(r) <- extent(1,ncol(m),1,nrow(m)); b <- as.matrix(buffer(r,1)); m[is.na(m) & !is.na(b)] <- 0; m; };

set.seed(1L); m <- matrix(rep(NA,100),nrow=10L); m[sample(1:100,10L)] <- 1;

expected <- bgoldst1(m);
identical(expected,bgoldst2(m));
## [1] TRUE
identical(expected,geotheory(m));
## [1] TRUE

microbenchmark(bgoldst1(m),bgoldst2(m),geotheory(m));
## Unit: microseconds
##          expr      min        lq      mean    median        uq      max neval
##   bgoldst1(m)   89.380   96.0085  110.0142  107.9825  119.1015  197.149   100
##   bgoldst2(m)   87.242   97.5055  111.4725  107.3410  121.2410  176.194   100
##  geotheory(m) 5010.376 5519.7095 6017.3685 5824.4115 6289.9115 9013.201   100

set.seed(1L); NR <- 100L; NC <- 100L; probNA <- 0.9; m <- matrix(sample(c(1,NA),NR*NC,T,c(1-probNA,probNA)),NR);

expected <- bgoldst1(m);
identical(expected,bgoldst2(m));
## [1] TRUE
identical(expected,geotheory(m));
## [1] TRUE

microbenchmark(bgoldst1(m),bgoldst2(m),geotheory(m));
## Unit: milliseconds
##          expr       min        lq      mean    median        uq        max neval
##   bgoldst1(m)  6.815069  7.053484  7.265562  7.100954  7.220269   8.930236   100
##   bgoldst2(m)  6.920270  7.071018  7.381712  7.127683  7.217275  16.034825   100
##  geotheory(m) 56.505277 57.989872 66.803291 58.494288 59.451588 571.142534   100