条件矩阵邻接计算

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矩阵方法来实现这一点,没有一些可怕的行和逐行循环算法？

---

NB.我重新修改了这个例子,不那么模棱两可.我需要非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)。

逻辑连接的右侧是最棘手的部分。我们需要执行 4 项测试，每个邻接方向一项。通用算法是：

1：从与该邻接方向的任何其他索引不相邻的邻接方向的奇异索引中索引。例如，对于“右方向”，我们从最左边的列建立索引，因为它不在任何其他索引的右侧。换句话说，不存在最左边的列在其右侧的列，因此我们可以忽略它（并且必须将其删除）以进行“正确方向”计算。

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

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

4： 4 个测试的逻辑与。结果将是一个逻辑矩阵，其中每个TRUE相邻单元格中都有 NA 的元素。

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

请注意，还有另一种方法可以做到这一点，这可能稍微更直观。我们可以编写 4 个测试中的每一个来测试非 NA（而不是 NA），然后将它们逻辑或在一起。这也需要绑定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