Animating simple cellular automata with Python
Recently, I had a talk on Prim’s algorithm at a meeting of our local algorithmic club. While preparing for the talk, I created several animations of the maze generating process using matplotlib. I can’t say it was a fun time. Animating several thousands of frames took ages, not to mention that the API is rather lowlevel.
Agreed, matplotlib is well suited for producing short animations of plots, but for combining a large number of images into a movie or animated GIF I needed to find another tool. I looked around and found an excellent MoviePy library.
At the next meeting of our club, we will be delving deeper into generating mazes, this time using cellular automata. This makes it a perfect opportunity to describe here how one can implement such automata in Python and visualize their evolution with beautiful video clips produced by MoviePy.
Here’s what you’ll find in today’s post.
 Introduction to cellular automata.
 Counting neighbours as a 2D convolution
 Implementing evolution of cellular automata using numpy stack
 Drawing state of an automaton using matplotlib
 Creating animation of evolving automata using MoviePy
 Gallery of nice examples produced by our scripts
Prerequisites
To run the examples from this post you’ll need several packages. I used Python 3.7 with:
scipy
1.3.1numpy
1.17.3matplotlib
3.0.2moviepy
1.0.1
Introduction to cellular automata
Consider a square grid, where each cell can be in one of some predefined states. For this post we will restrict ourselves to automata with only two possible states, that is a cell can be either alive or dead. To each cell, one can assign its neighbourhood, i.e. set of cells that are directly in contact (via border or a corner) with it. The figure below shows a 4x5 grid with a red cell along with its neighbours.
Note the special case of the cells lying at the edge or in the corner of the grid. We have basically two choices:
 considering only the “real” neighbours. In this approach, the neighbourhood of such cells is smaller than that of ones lying in the interior.
 virtually connecting the opposite edges. This wraps our grid, making it a torus.
Below we can see how the neighbourhood of the corner cell looks like depending whether we wrap the grid or not.
We start with some initial state of each cell in the grid. We can transform the grid by applying a set of rules determining the next state of each cell as a function of its neighbours’ states.
For instance, the set of rules below defines a cellular automaton known as Conway’s Game of Life:
 If a dead cell has at least three neighbours, it becomes alive.
 If an alive cell has two or three neighbours it stays alive.
 All other alive cells become dead. That is, cells with less than two or more than four alive neighbours die.
One can see how the above rules describe a very simplified evolution of some primitive population. Each member of the population needs just the right number of neighbours to survive, otherwise, it dies as a result of under or overpopulation. New members of the population may also be born if the conditions are favourable.
If we restrict ourselves only to rules based on the number of neighbours, the behaviour of automaton can be completely specified by its initial state and two lists:
 list of numbers of neighbours defining when the dead cell becomes alive (in Game of Life this is [3]),
 list of numbers of neighbours defining when the living cell survives (in Game of Life this is [2, 3]).
Cells that are not born and do not survive can be considered dead since they either were dead before and haven’t become alive, or they died because they didn’t match survival criterion. One can write the rules concisely using a socalled rule string. There exist several formats of rule strings, but probably the most widely used takes the form Bx/Sy, where x denotes all numbers from the first list and y denotes all numbers from the second one. For instance, in this notation the rule string for Game of Life is B3/S23.
It would be now straightforward to implement a simple automaton using Python. However, before we do that, let’s discuss an alternative view on the neighbourcounting process that will allow us to vectorize some computations.
Counting neighbours as a 2D convolution
Let’s represent our grid as a matrix with 0s and 1s marking dead and alive cells respectively. Our goal is to produce a matrix with the same size as an input grid, but with all entries equal to the number of alive neighbours of the corresponding cell. For instance, assuming that we don’t wrap edges, this is an example of transformation we wish to perform:
And here is how the same transformation looks like if we glue the opposite edges:
Let’s forget about the edge cells for a moment and focus only on the ones in the interior. It turns out that for those cells the process boils down to the following steps:

Define matrix $K$:
 Create an empty matrix $M$ of the same size as an input grid.
 For each cell in the interior of the input grid:
 Position matrix $K$ on top of the input grid such that the current cell lies below $0$.
 Multiply elements of $K$ with corresponding elements of the grid.
 Sum the results of this multiplication and place it in the corresponding place in $M$.
See the equation below for an example of the above process.
As you may already see, what we perform here is a 2D convolution of our grid with $K$ as a kernel. Note that normally when performing convolution we reverse the order of rows and columns in a kernel, but $K$ is invariant under such reversal.
This is great because 2D convolution is already implemented in SciPy stack. Indeed, running the following code produces the result as in the example above.
The mode="same"
tells convolve2d
to return an array of the same shape as the grid, otherwise,
the full output larger than the grid is returned. Specifying boundary
determines how the boundary cells
are handled. There are several options, but the ones we are interested in are the following:
"fill"
: virtually extend the grid (by default with zeros), so the convolution can be performed as for interior elements."wrap"
: glue opposite edges of the grid together.
Now that we know all the necessary tools, we can move on to the actual implementation of our automata.
Implementing evolution of cellular automata using numpy
Let us first write down the requirements for our implementation. While it might be tempting to write
a class like CellularAutomaton
, I think a cleaner solution is to just write a generator the yields
consecutive states of the grid. What input data do we need?
 An initial state of the grid.
 The rules (specified by a rule string).
 A boolean indicating whether we should wrap the edges.
Therefore, our generator could have a signature like this:
Before we implement it, however, let’s write a simple helper that parses a rule string and produces two lists of integers  numbers of neighbours needed for birth and survival of a cell.
This is pretty simple:
 We first match the input string to the suitable regular expression. This expression matches
exactly those strings that have a letter B, followed by some number of digits, then by a
slash, then by a letter S and finally again by some number of digits. Note that any of
the groups of digits can be empty because we used
*
instead of+
.  If the expression matched, we have exactly two groups in the match  they correspond to the first and the second sequence of digits.
 We convert each digit in those two groups into an integer, and in the end, return two lists of numbers.
Now let’s see how to evolve the grid. If we had boolean arrays marking the newly born cells
and the surviving cells, then the new state is simply a logical elementwise or
between them.
Obtaining those arrays is straightforward using NumPy’s isin
function, which is
a vectorized version of the builtin in
operator. The following example shows how it works.
import numpy as np
arr = np.array([[1, 0, 2], [4, 2, 1]])
print(np.isin(arr, [1, 2]))
[[ True False True]
[False True True]]
isin
function in action. The output array has the same shape as arr
,
and each of its entry tells whether the corresponding element in arr
is equal to 1 or 2.
A boolean array marking newly cells can be constructed as follows:
 Compute an array with neighbour counts.
 Using
isin
, construct a boolean array marking cells that have sufficient neighbours to be born.  Use logical
and
to combine this boolean array with the negation of the current grid.
The last step is necessary because only the cells that were previously dead can be born. Marking surviving cells is similar, the only significant difference is that we combine it with the grid instead of its negation (because only already living cells can survive).
Here’s how the code for evolving our automaton might look like.
from functools import reduce
import numpy as np
from scipy.signal import convolve2d
def evolve(initial_state: np.ndarray, rules: str, wrap: bool = False) > np.ndarray:
birth_list, survival_list = parse_rules(rules)
kernel = np.array([[1, 1, 1], [1, 0, 1], [1, 1, 1]])
grid = np.array(initial_state)
boundary = "wrap" if wrap else "fill"
while True:
yield grid
neighbours = convolve2d(grid, kernel, mode="same", boundary=boundary)
birth_mask = np.isin(neighbours, birth_list)
survival_mask = np.isin(neighbours, survival_list)
grid = (birth_mask & ~grid)  (survival_mask & grid)
We can easily verify that the above code produces correct results, but it would require manual inspection of the produced arrays. This is quite boring, so let’s focus on visualizing our automata.
Drawing state of an automaton using matplotlib
Drawing state of our automaton is simple. Basically, we draw our image using imshow
and add grid positioned at
halfinteger marks. We also remove any other elements like spines and tick labels.
from matplotlib import pylab as plt
from matplotlib.cm import binary
def plot_state(state, axes, line_width=2):
axes.imshow(state, interpolation="none", cmap=binary, zorder=0)
axes.set_xticks(np.arange(0.5, state.shape[1]+0.5, 1), minor=True)
axes.set_yticks(np.arange(0.5, state.shape[0]+0.5, 1), minor=True)
axes.yaxis.grid(True, which='minor', linewidth=line_width)
axes.xaxis.grid(True, which='minor', linewidth=line_width)
for spine in ["left", "right", "top", "bottom"]:
axes.spines[spine].set_visible(False)
for axis in (axes.xaxis, axes.yaxis):
axis.set_ticklabels([])
axis.set_ticks_position('none')
In the below script we utilize plot_state
function to plot several states of some automaton.
from itertools import islice
initial_state = np.array(
[
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0],
[0, 0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 0, 0, 0],
[0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
]
)
states = list(islice(evolve(initial_state, "B3/S23"), 8, 11, None))
fig, axes = plt.subplots(1, 3, figsize=(16, 3))
for state, axis in zip(states, axes):
plot_state(state, axis)
plt.show()
This looks pretty nice, but it would be much nicer if we created an animated sequence of images or a movie clip of the automaton progressing through its states, which is what we’ll do next.
Creating animation of evolving automata using MoviePy
To produce an animation using MoviePy we first need to create a list of frames that will be combined
into a single clip. We could construct required frames ourselves using figure’s canvas tostring_rgb()
method but luckily MoviePy has us covered and provides mplfig_to_npimage
specifically for that purpose.
Once we have a list of frames, we use it to initialize an instance of ImageSequenceClip
. In turn,
this instance can be used to produce an animated gif or video file, or even displayed in Jupyter Notebook.
The process is shown in the snippet below.
from itertools import islice
from moviepy.editor import ImageSequenceClip
from moviepy.video.io.bindings import mplfig_to_npimage
def animate_automaton(
states: np.ndarray, line_width=2, fps=2, **fig_kwargs
) > ImageSequenceClip:
fig, axes = plt.subplots(1, 1, **fig_kwargs)
frames = []
for state in states:
plot_state(state, axes, line_width=line_width)
frames.append(mplfig_to_npimage(fig))
return ImageSequenceClip(frames, fps=fps)
initial_state = np.array([
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0],
[0, 0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 0, 0, 0],
[0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
])
states = list(islice(evolve(initial_state, "B3/S23"), 15))
clip = animate_automaton(states, figsize=(16, 9), dpi=30, fps=2, line_width=4)
clip.write_gif("automaton.gif")
islice
to select only
first fifteen first states of the evolution. The write_gif
method is used
to save our animation as a gif, alternatively write_videofile
could be used
to save animation as a video file.
Gallery of examples
Depite their simplicity, behaviour of cellular automata can be surprisingly rich and complex. So, to conclude, lets see some nice example animations.