Here I demonstrate a simple way to code Conway's game of life (GoL) in R and to produce the animation above. Cellular automata in R are usually painfully slow if you iterate through all grid cells in an array. A couple of years ago my friend **Martin Weiser** came with an idea to avoid the individual iterations. He suggested to make GoL in R fast by taking the advantage of the somewhat optimized matrix operations. This is how it works:

**1.** You have an rectangular array (matrix) of living (1) and dead (0) cells. You need to figure out what is the number **Ni** of living cells in the neighborhood of each focal cell **i**.

**2.** Once you know **Ni**, you change the state of the focal cell **i** accordingly (to 0 or to 1; using the rules of GoL). However, iterating through all cells **i** would take ages in R. What you first need to do is to **make 8 copies of the original array**. Then:

- In copy 1 you delete the leftmost column and add column of zeroes to the very right.

- In copy 2 you do the opposite.

- In copy 3 you delete the uppermost row and add row of zeroes to the bottom.

- In copy 4 you delete the row at the bottom and add a row of zeroes to the very top.

So you are essentially shifting the copies of the original array by 1 position to the left, right, up and down.

**3.** Then you do similar thing with copies 5-8 but now you shift these by 1 position diagonally. I am sure you can figure it out.

**4.** Next, you stack all of the copies ("sum them up" using the + operator in R) and this will give you a new array of **Ni** values.

**5.** In the next step you use logical subscripting to apply the rules of GoL (conditional on the array of **Ni** values) to the original array.

**6.** You go to step 1.

In R this implementation works MUCH faster than any looping solution. In MatLab you would do this by using the the **shift** function and so your code would be even simpler. Also note that in reality the processor does actually iterate through all grid cells; the solution presented here does not make all of the iterations happen at the same time. My solution only hides the iteration process by using the R-optimized form of matrix operations and logical subscripting.

See the code below for how it is actually done. Note that you may run out of memory if you specify a very large array. Also, make sure to install the **caTools** library. Have fun!

# library that allows the animated .gif export library(caTools) # The game.of.life() function ------------------ # Arguments: # side - side of the game of life arena (matrix) # steps - number of animation steps # filename - name of the animated gif file game.of.life <- function(side, steps, filename){ # the sideXside matrix, filled up with binomially # distributed individuals X <- matrix(nrow=side, ncol=side) X[] <- rbinom(side^2,1,0.4) # array that stores all of the simulation steps # (so that it can be exported as a gif) storage <- array(0, c(side, side, steps)) # the simulation for (i in 1:steps) { # make the shifted copies of the original array allW = cbind( rep(0,side) , X[,-side] ) allNW = rbind(rep(0,side),cbind(rep(0,side-1),X[-side,-side])) allN = rbind(rep(0,side),X[-side,]) allNE = rbind(rep(0,side),cbind(X[-side,-1],rep(0,side-1))) allE = cbind(X[,-1],rep(0,side)) allSE = rbind(cbind(X[-1,-1],rep(0,side-1)),rep(0,side)) allS = rbind(X[-1,],rep(0,side)) allSW = rbind(cbind(rep(0,side-1),X[-1,-side]),rep(0,side)) # summation of the matrices X2 <- allW + allNW + allN + allNE + allE + allSE + allS + allSW # the rules of GoL are applied using logical subscripting X3 <- X X3[X==0 & X2==3] <- 1 X3[X==1 & X2<2] <- 0 X3[X==1 & X2>3] <- 0 X <- X3 # each simulation step is stored storage[,,i] <- X2 # note that I am storing the array of Ni values - # - this is in order to make the animation prettier } storage <- storage/max(storage) # scaling the results # to a 0-1 scale # writing the results into an animated gif write.gif(storage, filename, col="jet", delay=5) } game.of.life(side=150, steps=300, file="conway.gif")

This version doesn't wrap around the boundaries of the grid. You can replace the columns or rows of zeros with the respective columns that are currently deleted in your version to run the GOL on a torus rather than on a plane, e.g.

allW <- cbind(grid[,-1],grid[,1])

allN <- rbind(grid[-1,],grid[1,])

allE <- cbind(grid[,grid.size],grid[,-grid.size])

allS <- rbind(grid[grid.size,],grid[-grid.size,])

allNW<-rbind(allW[-1,],allW[1,])

allNE<-rbind(allE[-1,],allE[1,])

allSW<-rbind(allW[grid.size,],allW[-grid.size,])

allSE<-rbind(allE[grid.size,],allE[-grid.size,])

where grid.size is the size of one axis of the grid and grid is the matrix.

Cheers

Rob Knell

Pingback: Getting started with R | Conservation planning for migratory species

Thanks ! I found this code very useful in recent work I've been doing developing a simulation of tsetse fly populations. In a recent blog post I describe how I modify the code to create movement functions that accept 1) a matrix of the spatial distribution of a population, and 2) the proportion that move.

Best wishes,

Andy

Thanks! I am glad that it's been useful