2
\$\begingroup\$

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.

\$\endgroup\$
4
\$\begingroup\$

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
\$\endgroup\$
  • 1
    \$\begingroup\$ This is pretty neat! Could you help explain what is happening? \$\endgroup\$ – masfenix Jun 26 '16 at 18:17
  • \$\begingroup\$ 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. \$\endgroup\$ – flodel Jun 26 '16 at 18:31
  • \$\begingroup\$ 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. \$\endgroup\$ – masfenix Jun 26 '16 at 18:33
  • \$\begingroup\$ 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++. \$\endgroup\$ – flodel Jun 26 '16 at 18:39
  • \$\begingroup\$ 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 \$\endgroup\$ – masfenix Jun 26 '16 at 19:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.