R - Vectorize nested for loops to assign new values to a matrix

Question

I'm currently trying to vectorize this nested for loop to save time during execution but it doesnt seem to work. What I want is to go through every cell of my matrix and check if the value is 0 or 1 then change the value based on a condition. This is the algorithm for the Forest-fire model

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
      }
    }
    
  }
}

Here is the neighbours function :

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)
}

And some code to make it work :

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

I've started by trying with sapply but couldn't get results back into the matrix and for the past hour I've been trying to use nested foreach loops

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)
      }
    }
  }

But I'm only getting back the lines I need to change I'm not familiar with vectorization so perhaps I'm missing some basic steps !

user_11355066 · Answer

As you seem very comfortable in loops, you may want to recode in rcpp. Your code will translate very quickly.

Here's a draft that goes towards more efficient in R and is about 2.5 times more efficient on this small dataset.

## 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}
        )
    }
}

It's likely possible that the outer loop can be removed but additional thought would be needed for the neighbors algorithm to allow more than one i at a time.

R - Vectorize nested for loops to assign new values to a matrix

1 Answers