Fast Conway’s game of life in R

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

# The function ------------------
# Arguments:
# side - side of the game of life arena (matrix)
# steps - number of animation steps
# filename - name of the animated gif file <- 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)
}, steps=300, file="conway.gif")

4 Responses

  1. Rob Knell
    Rob Knell March 14, 2014 at 9:07 am |

    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,])

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


    Rob Knell

  2. Getting started with R | Conservation planning for migratory species

    […] going to show you how to get one of these animations using R based on some code on the blog Are you cereal? that I’ve altered. This is what it should look like at the […]

  3. Andy South
    Andy South February 5, 2015 at 6:41 pm |

    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,


Leave a Reply