# Simulation of animal skins

This is a simulation of animal skins, like a cellular automat. I'd like it improved.

#Create Matrix with random numbers 0/1 with a 50% chance
createMatrix <- function(br,gen) {
mat <- matrix(rbinom(br*gen,1,0.5),br,gen)
return(mat)
}

nextMatrix <- function(mat,w) {

#Make that Matrix continuous
wideMatrix <- cbind(mat,mat,mat)
bigMatrix <- rbind(wideMatrix,wideMatrix,wideMatrix)

newMat <- mat

for(j in 1:nrow(mat)){
for(i in 1:ncol(mat)){
act <- sum(bigMatrix[(197+j):(203+j),(197+i):(203+i)]) #radius of 3
inh <- sum(bigMatrix[(195+j):(206+j),(195+i):(206+i)])*w #radius of 6
if(act>inh) { newMat[j,i] <- 1 }
if(act<inh) { newMat[j,i] <- 0 }
}
}
return(newMat)
}

x<- createMatrix(200,200)

for(i in 1:5) {
x <- nextMatrix(x,0.33)
}

image(x, axes=FALSE,col = c("black","darkgoldenrod"))


In particular, the part where I sort out the part with the radius within a matrix:

 act <- sum(bigMatrix[(197+j):(203+j),(197+i):(203+i)])
inh <- sum(bigMatrix[(195+j):(206+j),(195+i):(206+i)])*w


Is there any package I could use to get a "real radius" (it is a quadrat in my example)?

• What exactly needs to be improved: speed, flexibility, mathematical rigor, or something else entirely? – r2evans Dec 22 '14 at 7:51
• for example the 2 for loops, is there a way i could to it with apply instead of for? – Siklos Dec 23 '14 at 8:19
• With the specific nature of the loops, I don't know if that will improve speed much or clarity. It would benefit considerably from using Rcpp, but I don't usually go there until/unless speed is an issue. – r2evans Dec 23 '14 at 14:52

You can use roll_sum from the RcppRoll package to calculate rolling sums. That way I got an 80-fold speed increase.

Also from a memory perspective, your bigMatrix can be a lot smaller.

Below is my version of nextMatrix.

require(RcppRoll)
nextMatrix2 <- function(mat,w) {
# Make continuous matrix only as large as necessary
wideMatrix <- cbind(mat[, ncol(mat)-4:0], mat, mat[, 1:6])
bigMatrix <- rbind(wideMatrix[nrow(mat)-4:0, ] , wideMatrix, wideMatrix[1:6, ])
# use roll_sum from RcppRoll to get act/inh as matrices
actMat <- roll_sum(roll_sum(bigMatrix[3:208, 3:208], 7), 7, by.column=FALSE)
inhMat <- roll_sum(roll_sum(bigMatrix, 12), 12, by.column=FALSE)*w
# create the new matrix
newMat <- actMat > inhMat
newMat[actMat == inhMat] <- mat[actMat == inhMat]
# retrun newMat
return(newMat)
}


And the benchmarktest:

# load benchmarking package
require(microbenchmark)
# create data
x <- createMatrix(200,200)
y <- x
# benchmarktests
microbenchmark(
x <- nextMatrix(x, 0.33)
,
y <- nextMatrix2(y, 0.33)
)
## Unit: milliseconds
##                      expr        min         lq       mean    median         uq       max neval cld
## x <- nextMatrix(x, 0.33)  545.519907 566.455440 584.396621 577.26897 590.373952 687.29806   100   b
## y <- nextMatrix2(y, 0.33)   5.858406   6.131756   7.214752   6.24489   6.453976  62.63813   100  a
#
# checking for equality
all.equal(x, y)
##  TRUE