Counting matrix elements that have "4" or "5" as a neighbor

Question

I have a matrix (500x500) of integers. For each entry, I need to look at its surrounding neighbours (so 8 elements) and determine what the integers are, and run a function on these integers.

Here what my function looks like right now (R code with the RCpp package).

fx2<- cppFunction('NumericMatrix getNeighbours(NumericMatrix x) {
            int nrow = x.nrow(), ncol = x.ncol();

            //create a result matrix of probabilities
            NumericMatrix outProb(nrow, ncol);

            // some counters needed to calculate probability 
            int kc = 0; int kpre = 0; int ks = 0;

            // i and j loop through the elements one by one. 
            //  - ignore the boundary as I dont know how to handle it for now
            for (int i = 1; i < (nrow-1); i++) {
              for (int j = 1; j < (ncol-1); j++) {

                // go through the neighbours
                for(int k = -1; k <= 1; k++){
                  for(int l = -1; l <= 1; l++){

                    // check what values the neighbours are

                    // increment the kc
                    if(x(i + k, j + l) == 4){
                      kc++;
                    }

                    // increment the ks counter
                    if(x(i + k, j + l) == 5){
                      ks++;
                    }

                    // more if statements removed for readibility

                  }
                } //end of loop for neighbours

                // calculate a probability for the resulting matrix
                outProb(i, j) = functionof(kc, ks, kpre);                
              }              
            }
            return outProb;
            }')

I am looking to make this code faster/more efficient.

Example of result: Suppose the functionof just doubles the value of the counter.

> health
     [,1] [,2] [,3] [,4] [,5]
[1,]    4    1    1    1    1
[2,]    1    1    1    1    1
[3,]    1    1    1    1    1
[4,]    1    1    1    1    1
[5,]    1    1    1    1    1

We see that (2, 1), (1, 2) and (2, 2) have 4 as their neighbours. So the output would be

     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    0    2    0    0    0    0
[2,]    2    2    0    0    0    0
[3,]    0    0    0    0    0    0
[4,]    0    0    0    0    0    0
[5,]    0    0    0    0    0    0
[6,]    0    0    0    0    0    0

Note that there are other elements in the matrix that have 4 as a neighbour. I've omitted that from the result matrix above.

I currently wrote the code in R but even after vectorizing/parallelizing it is still very slow. So now I've moved on to a C++ function that I can call using R's RCpp package.

Answer 1

This is how you could write your code efficiently in base R:

offset <- function(x, v_offset = 0L, h_offset = 0L) {
   stopifnot(h_offset %in% -1:1, v_offset %in% -1:1)
   nr <- nrow(x)
   nc <- ncol(x)
   if (h_offset == +1L) x <- cbind(x[, -1], -99L)
   if (h_offset == -1L) x <- cbind(-99L, x[, -nc])
   if (v_offset == +1L) x <- rbind(x[-1, ], -99L)
   if (v_offset == -1L) x <- rbind(-99L, x[-nr, ])
   x
}

getNeighbours <- function(x, functionof) {
   v_offsets <- c(-1, -1, -1,  0,  0, +1, +1, +1)
   h_offsets <- c(-1,  0, +1, -1, +1, -1,  0, +1)
   neigh_x   <- Map(offset, list(x), v_offsets, h_offsets)
   kc   <- Reduce(`+`, Map(function(x) x == 4L, neigh_x))
   ks   <- Reduce(`+`, Map(function(x) x == 5L, neigh_x))
   kpre <- 0
   functionof(kc, ks, kpre)
}

It builds a list (neigh_x) of eight matrices in memory, one for each neighboring position. Instead of four loops, you are down to a single loop (8 iterations) wherever you see Map or Reduce, everything else uses vectorized operations.

Testing the code:

health <- matrix(1, 5, 5)
health[1, 1] <- 4
test_fun <- function(kc, ks, kpre) 2 * (kc + ks + kpre)

getNeighbours(health, functionof = test_fun)
#      [,1] [,2] [,3] [,4] [,5]
# [1,]    0    2    0    0    0
# [2,]    2    2    0    0    0
# [3,]    0    0    0    0    0
# [4,]    0    0    0    0    0
# [5,]    0    0    0    0    0