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 &lt;- function(br,gen) {
mat &lt;- matrix(rbinom(br*gen,1,0.5),br,gen)  
    return(mat)
}

nextMatrix &lt;- function(mat,w) {

    #Make that Matrix continuous 
    wideMatrix &lt;- cbind(mat,mat,mat)
    bigMatrix &lt;- rbind(wideMatrix,wideMatrix,wideMatrix)

    newMat &lt;- mat                    

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

x&lt;- createMatrix(200,200)

for(i in 1:5) {
    x &lt;- 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 &lt;- sum(bigMatrix[(197+j):(203+j),(197+i):(203+i)])
 inh &lt;- 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)?
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 &lt;- function(mat,w) {
  # Make continuous matrix only as large as necessary
  wideMatrix &lt;- cbind(mat[, ncol(mat)-4:0], mat, mat[, 1:6])
  bigMatrix &lt;- rbind(wideMatrix[nrow(mat)-4:0, ] , wideMatrix, wideMatrix[1:6, ])
  # use roll_sum from RcppRoll to get act/inh as matrices
  actMat &lt;- roll_sum(roll_sum(bigMatrix[3:208, 3:208], 7), 7, by.column=FALSE)
  inhMat &lt;- roll_sum(roll_sum(bigMatrix, 12), 12, by.column=FALSE)*w
  # create the new matrix 
  newMat &lt;- actMat &gt; inhMat
  newMat[actMat == inhMat] &lt;- mat[actMat == inhMat]
  # retrun newMat
  return(newMat)
}


And the benchmarktest:

# load benchmarking package
require(microbenchmark)
# create data
x &lt;- createMatrix(200,200)
y &lt;- x
# benchmarktests
microbenchmark(
  x &lt;- nextMatrix(x, 0.33)
  ,
  y &lt;- nextMatrix2(y, 0.33)
)
## Unit: milliseconds
##                      expr        min         lq       mean    median         uq       max neval cld
## x &lt;- nextMatrix(x, 0.33)  545.519907 566.455440 584.396621 577.26897 590.373952 687.29806   100   b
## y &lt;- 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