Counting matrix elements that have “4” or “5” as a neighbor

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.

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

• This is pretty neat! Could you help explain what is happening? – masfenix Jun 26 '16 at 18:17
• Maybe this will help understand. You probably know that x == 4 returns a matrix showing you what cells in x have a value of 4. Well, with the offset function I built, you can for example do offset(x, h_offset = 1) == 4 to find what cells have a 4 to their right (horizontal offset of 1) . neigh_x is a list of eight matrices computing the shifted versions of x. Map then helps you loop over that list to return a list of eight boolean matrices where the neighboors are (eg.) equal to 4. The Reduce(+, ...) turns these 8 matrices into their sum. I bet you'll figure out the rest. – flodel Jun 26 '16 at 18:31
• Ahh thanks, that makes sense! Two more quick questions: Is there a way to parallelize this code? I have a large cluster I am going to run this on. Second, are there libraries in C++ that can help achieve this functional method. I'd like to compare the speeds. – masfenix Jun 26 '16 at 18:33
• parallelize: do you have one gigantic matrix or a gigantic number of matrices? If the former, the single (Map/Reduce) loop has only 8 iterations, so I don't see a huge need to parallelize if your matrix already holds in memory. If the latter, sure, you could search the internet. Equivalent c++ code: isn't what you did already sufficient? It is probably already faster than my R solution. R can't really compete with c++. – flodel Jun 26 '16 at 18:39
• There is a slight modification I am trying to make but can't. I only want to run this "logic" of finding neighbours if the element in the matrix is 1 or 2. Is there a way to modify this somehow. I am failing to see where I can modify this – masfenix Jun 26 '16 at 19:11