R - 向量化嵌套 for 循环以将新值分配给矩阵

Question

我目前正在尝试对这个嵌套的 for 循环进行矢量化以在执行期间节省时间，但它似乎不起作用。我想要的是遍历矩阵的每个单元格并检查值是 0 还是 1，然后根据条件更改值。这是森林火灾模型的算法

for (i in 1:nrow(X)) {
  for (j in 1:ncol(X)) {
    
    if (X[i, j] == 2) {
      if (runif(1) > (1 - a)^neighbours(X, i, j)) {
        B[i, j] <- 1
      }
    } 
    else if (X[i, j] == 1) {
      burning <- TRUE
      if (runif(1) < b) {
        B[i, j] <- 0
      }
    }
    
  }
}

这里是邻居函数：

neighbours <- function(A, i, j) {
  # calculate number of neighbours of A[i,j] that are infected
  # we have to check for the edge of the grid
  nbrs <- 0
  # sum across row i - 1
  if (i > 1) {
    if (j > 1) nbrs <- nbrs + (A[i-1, j-1] == 1)
    nbrs <- nbrs + (A[i-1, j] == 1)
    if (j < ncol(A)) nbrs <- nbrs + (A[i-1, j+1] == 1)
  }
  # sum across row i
  if (j > 1) nbrs <- nbrs + (A[i, j-1] == 1)
  nbrs <- nbrs + (A[i, j] == 1)
  if (j < ncol(A)) nbrs <- nbrs + (A[i, j+1] == 1)
  # sum across row i + 1
  if (i < nrow(A)) {
    if (j > 1) nbrs <- nbrs + (A[i+1, j-1] == 1)
    nbrs <- nbrs + (A[i+1, j] == 1)
    if (j < ncol(A)) nbrs <- nbrs + (A[i+1, j+1] == 1)
  }
  return(nbrs)
}

还有一些让它工作的代码：

set.seed(3)
X <- matrix(2, 21, 21)
X[11, 11:13] <- 1
burning <- FALSE
a= 0.2
b = 0.4
B <- X

我开始尝试使用 sapply，但无法将结果返回到矩阵中，过去一个小时我一直在尝试使用嵌套的 foreach 循环

library(foreach)
B <-
foreach(i=1:nrow(X), .combine='cbind') %:%
  foreach(j=1:ncol(X), .combine='c') %do% {
    if (X[i, j] == 2) {
      if (runif(1) > (1 - a)^neighbours(X, i, j)) {
        1
      }
    } 
    else if (X[i, j] == 1) {
      burning <- TRUE
      if (runif(1) < b) {
        0
        print(i)
        print(j)
      }
    }
  }

但我只是恢复了我需要更改的线路 我不熟悉矢量化，所以也许我错过了一些基本步骤！

Answer 1

由于您似乎对循环非常熟悉，因此您可能需要在 rcpp 中重新编码。您的代码将很快翻译。

这是一个在 R 中提高效率的草稿，在这个小数据集上效率提高了大约 2.5 倍。

## get constants out of the away above the loop
A = X == 1L
nr = nrow(X)
nc = ncol(X)

for (i in 1:nr) {
    i_start = i - (i > 1L)
    i_stop = i + (i < nr)
    for (j in 1:nc) {
        j_start = j - (j > 1L)
        j_stop = j + (j < nc)
        switch(X[i, j],
             ##refactoring of neighbours function which is a partial rolling sum of sub-matrixes equal to 1.
               2, if (runif(1) > (1 - a)^sum(A[i_start:i_stop, j_start:j_stop]))  B[i, j] <- 1,
               1, {burning = TRUE
                   if (runif(1) < b) B[i, j] = 0}
        )
    }
}

可能可以移除外循环，但需要额外考虑邻居算法一次允许多个i。