條件矩陣鄰接計算

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值爲零。是否有任何可能的矩陣方式來實現這一點，沒有一些可怕的行方式和col-wise循環算法？

---

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 

---

該解決方案的工作原理如下。

我們構造了一個邏輯索引矩陣，其中元素是NA AND與至少一個非NA元素（的上面，下面，左邊或右邊）相鄰。然後，我們可以用邏輯索引矩陣來下標m，並分配所需的替換值。

邏輯連接的LHS很容易;它只是is.na(m)。

邏輯連接的RHS是最棘手的部分。我們需要執行4個測試，每個方向的鄰接關係。一般算法如下：

1：索引關於鄰接方向的奇異索引，該奇異索引與任何其他索引相對於該鄰接方向不相鄰。例如，對於「正確的方向」，我們索引最左邊的列，因爲它不在任何其他索引的右側。換句話說，沒有列的右邊最左邊的列，所以我們可以忽略它（並且必須移除它）以進行「正確的方向」計算。

2：使用is.na()測試NA的子矩陣。

3：然後，我們必須綁定（cbind()用於水平鄰接方向，rbind()爲垂直）TRUE在相反側（即，在步驟1中除去索引的對面）所得的邏輯子矩陣的。這有效地導致鄰接方向上的最後一個索引總是在它的鄰接方向上具有（僞）NA，所以它將永遠不會被替換，因爲該鄰接方向。

4：邏輯AND的4個測試。對於中的每個鄰近單元的元素，結果將爲TRUE的邏輯矩陣。

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 

Answer 2

另一種方法：

require(raster) 
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 
     [,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