3
\$\begingroup\$

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

\$\endgroup\$
  • \$\begingroup\$ What exactly needs to be improved: speed, flexibility, mathematical rigor, or something else entirely? \$\endgroup\$ – r2evans Dec 22 '14 at 7:51
  • \$\begingroup\$ for example the 2 for loops, is there a way i could to it with apply instead of for? \$\endgroup\$ – Siklos Dec 23 '14 at 8:19
  • \$\begingroup\$ 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. \$\endgroup\$ – r2evans Dec 23 '14 at 14:52
2
\$\begingroup\$

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)
## [1] TRUE
\$\endgroup\$

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.