patternpythonMinor

Simulation of animal skins

Submitted by: @import:stackexchange-codereview·Mar 10, 2026·

Viewed 0 times

cellular-automata codereview r stackoverflow python matrix

animalsimulationskins

Problem

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

Solution

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

Code Snippets

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

# 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

Context

StackExchange Code Review Q#72099, answer score: 2

Revisions (0)

No revisions yet.