The Game of Life, Inheritance, and Cellular Automata¶
Preamble: Run the cells below to import the necessary Python packages
This notebook created by William Gilpin. Consult the course website for all content and GitHub repository for raw files and runnable online code.
import numpy as np
# Wipe all outputs from this notebook
from IPython.display import Image, clear_output, display
clear_output(True)
# Import local plotting functions and in-notebook display functions
import matplotlib.pyplot as plt
%matplotlib inline
# Add the course directory to the Python path
import sys
sys.path.append('../')
import cphy.plotting as cplot
The Game of Life and Cellular Automata¶
- Cellular automata are systems where all physical laws and fields are discretized: space is a lattice, and time is discrete
- Physical laws take the form of IF-THEN statements
- Basic cellular automata are Markovian: given current field values at some set of lattice points, perform an instantaneous update based on a set of rules
The Game of Life¶
- Conway, 1970
- Complex dynamics emerge from simple rules
- Markovian: next state of universe depends only on current state of universe
- Local: next state at a lattice point depends only on current site and nearest neighbors
- Synchronous: we update all sites at once, rather than raster scanning
- Finite fields: A site can have values of 1 (alive) or 0 (dead)
Game life rules:¶
- Underpopulation: Any live cell with fewer than two live neighbours dies
- Survival: Any live cell with two or three live neighbours lives stays alive
- Overpopulation: Any live cell with more than three live neighbours dies
- Reproduction: Any dead cell with exactly three live neighbours becomes a live cell
These simple rules give rise to surprisingly complex behaviors; in fact, the Game of Life has been shown to support universal computation, given a suitable encoding scheme based on initial conditions. We've included some examples of fixed-point, limit-cycle, and soliton like dynamics below
.
Images from Scholarpedia
Implementing Cellular Automata using inheritance in Python¶
Below, we are going to define a "base" class for cellular automata, and then define subclasses for various specific automata. The base class will handle state counting, initialization, and will contain a simulate
method for the simulation loop (where we repeatedly call the update rule).
In the base class, we will define a next_state
method because the simulate
method depends on it. But because the state rule is different for each automaton, we will leave the implementation of next_state
to be overwritten by the subclasses (the base class raises an exception)
This is an example of the Template Method Pattern.
# A base allows us to impose structure on specific cases. We will use this as a parent
# for subsequent classes
class CellularAutomaton:
"""
A base class for cellular automata. Subclasses must implement the step method.
Parameters
n (int): The number of cells in the system
n_states (int): The number of states in the system
random_state (None or int): The seed for the random number generator. If None,
the random number generator is not seeded.
initial_state (None or array): The initial state of the system. If None, a
random initial state is used.
"""
def __init__(self, n, n_states, random_state=None, initial_state=None):
self.n_states = n_states
self.n = n
self.random_state = random_state
np.random.seed(random_state)
## The universe is a 2D array of integers
if initial_state is None:
self.initial_state = np.random.choice(self.n_states, size=(self.n, self.n))
else:
self.initial_state = initial_state
self.state = self.initial_state
self.history = [self.state]
def next_state(self):
"""
Output the next state of the entire board
"""
return NotImplementedError
def simulate(self, n_steps):
"""
Iterate the dynamics for n_steps, and return the results as an array
"""
for i in range(n_steps):
self.state = self.next_state()
self.history.append(self.state)
return self.state
Just to show how inheritance works, we will start by defining two boring automata---an extinction automaton that sends every single input condition to the same output state, and a random automaton that randomly assigns each input state to either 0 or 1.
class ExtinctionAutomaton(CellularAutomaton):
"""
A cellular automaton that simulates the extinction of a species.
"""
def __init__(self, n, **kwargs):
super().__init__(n, 2, **kwargs)
def next_state(self):
"""
Output the next state of the entire board
"""
next_state = np.zeros_like(self.state)
return next_state
class RandomAutomaton(CellularAutomaton):
"""
A cellular automaton with random updates
"""
def __init__(self, n, random_state=None, **kwargs):
super().__init__(n, 2, **kwargs)
self.random_state = random_state
def next_state(self):
"""
Output the next state of the entire board
"""
next_state = np.random.choice(self.n_states, size=(self.n, self.n))
return next_state
model = ExtinctionAutomaton(40, random_state=0)
model.simulate(200)
plt.figure()
plt.imshow(model.initial_state, cmap="gray")
plt.title("Initial state")
plt.figure()
plt.imshow(model.state, cmap="gray")
plt.title("Final state")
Text(0.5, 1.0, 'Final state')
model = RandomAutomaton(40, random_state=0)
model.simulate(200)
plt.figure()
plt.imshow(model.initial_state, cmap="gray")
plt.title("Initial state")
plt.figure()
plt.imshow(model.state, cmap="gray")
plt.title("Final state")
Text(0.5, 1.0, 'Final state')
We are now ready to define the Game of Life as a subclass of the base class.
# This child class inherits methods from the parent class
class GameOfLife(CellularAutomaton):
"""
An implementation of Conway's Game of Life in Python
Args:
n (int): The number of cells in the system
**kwargs: Additional keyword arguments passed to the base CellularAutomaton class
"""
def __init__(self, n, **kwargs):
# the super method calls the parent class's __init__ method
# and passes the arguments to it. It is a constructor
super().__init__(n, 2, **kwargs)
def next_state(self):
"""
Compute the next state of the lattice
"""
# Compute the next state
next_state = np.zeros_like(self.state) # preallocate the next state
for i in range(self.n):
for j in range(self.n):
# Compute the number of neighbors
n_neighbors = np.sum(
self.state[i - 1:i + 2, j - 1:j + 2]
) - self.state[i, j]
# Update the next state
if self.state[i, j] == 1:
if n_neighbors in [2, 3]:
next_state[i, j] = 1
else:
if n_neighbors == 3:
next_state[i, j] = 1
return next_state
model = GameOfLife(40, random_state=0)
model.simulate(500)
# print timing results
print("Average time per iteration: ", end="")
%timeit model.next_state()
plt.figure()
plt.imshow(model.initial_state, cmap="gray")
plt.title("Initial state")
plt.figure()
plt.imshow(model.state, cmap="gray")
plt.title("Final state")
Average time per iteration: 3.7 ms ± 58.8 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
Text(0.5, 1.0, 'Final state')
from matplotlib.animation import FuncAnimation
from IPython.display import HTML
# Assuming frames is a numpy array with shape (num_frames, height, width)
frames = np.array(model.history).copy()
fig = plt.figure(figsize=(6, 6))
img = plt.imshow(frames[0], vmin=0, vmax=1, cmap="gray");
plt.xticks([]); plt.yticks([])
# tight margins
plt.margins(0,0)
plt.gca().xaxis.set_major_locator(plt.NullLocator())
def update(frame):
img.set_array(frame)
ani = FuncAnimation(fig, update, frames=frames, interval=50)
HTML(ani.to_jshtml())
How long did that take?¶
Memory Scaling¶
- To avoid 2D vs 3D, let's just say there are $N$ total sites. So our 2D universe is $\sqrt{N} \times \sqrt{N}$ in span, etc
- For synchronous update rules, each site value $s_{ij}(t+1)$ at timestep $t + 1$ depends on the lattice $s_{ij}(t)$. That means that each timestep requires a full copy of the entire universe, into which we write the next states one-by-one
- The original copy requires space $\sim N$, on top of the $\sim N$ memory requirements to store the initial lattice. So if our original lattice was 1 MP, then we need to set aside 2 MP of memory
- So-called "Big-O" notation usually ignores prefactors, and so we would say that the asymptotic space requirements to simulate the Game of Life are $\mathcal{O}(N)$
Runtime scaling¶
Let's look at what our algorithm does carefully: at each lattice site, we poll the nearest neighbors (including diagonal neighbors), and compute a sum. We then execute a conditional on that sum to write to the future state.
Since every site gets updated, we know that we have to perform at least $N$ operations in order to update the entire universe
For each site, let's keep a running sum of all neighbors, and a check of the current value. that's $9$ operations per site. For a generalized neighborhood, that's $K$ operations per site
So the runtime looks something like $N \times K$. Since $K \ll N$ we can treat it as a prefactor, and say that the runtime is $\mathcal{O}(N)$
However, what if we had a global cellular automaton, like an N-body simulation, where we need to compute pairwise forces between all lattice sites? In this limit $K \sim N$ because we have to check every other site every time we update a single site. So a global CA has runtime $\mathcal{O}(N^2)$
Low-level optimizations (vectorizations, convolution, etc) decrease the runtime by reducing the prefactor, but we cannot change the Big-O scaling properties, which are intrinstic to algorithms
Can we improve our implementation?¶
A key idea of cellular automata is locality: the update of each site depends only on the state of its neighbors. The Game of Life uses a Moore neighborhood, in which the first nearest neighbors (including diagonals) are considered. The Game of Life update rule can therefore be seen as first implementing convolution of the grid state $f$ with a kernel $g$ of the form
$$ g \equiv \begin{bmatrix} 1 & 1 & 1 \\ 1 & 0 & 1 \\ 1 & 1 & 1 \end{bmatrix} $$
This discrete convolution essentially takes the sum of the neighbors, and uses this sum plus the state of the center cell in order to decide what to do. This makes the Game of Life an example of an outer totalistic cellular automaton. Recall that a two-dimensional discrete convolution has the form
$$ (f * g)(x, y) = \sum_{m=-\infty}^{\infty} \sum_{n=-\infty}^{\infty} f(x-m, y-n) g(m, n) $$
However, when the kernel has compact support (as in our nearest-neighborhood update rule), this double summation simplifies,
$$ (f * g)(x, y) = \sum_{m=-K}^{K} \sum_{n=-K}^{K} f(x-m, y-n) g(m, n) $$
For our Game of Life example, $K=1$. The periodic boundary conditions can be implemented by padding each dimension of $f$ with values from the other side of the array, and then cropping the convolution result to the original size of $f$. We can estimate the number of multiplications required to perform a discrete convolution by counting the number of elements in the kernel, $M$, and the number of elements in the grid, $N$, for an overall time complexity of $\mathcal{O}(M\,N)$.
from scipy.signal import convolve2d
class GameOfLife(CellularAutomaton):
def __init__(self, n, **kwargs):
# the super method calls the parent class's __init__ method
# and passes the arguments to it. It is a constructor
super().__init__(n, 2, **kwargs)
def next_state(self):
"""
Compute the next state of the lattice using discrete convolutions to implement
the rules of the Game of Life.
"""
# Compute the next state
next_state = np.zeros_like(self.state)
# Compute the number of neighbors for each cell
totalistic_kernel = np.ones((3, 3))
totalistic_kernel[1, 1] = 0 # don't count the center cell itself
# convolve with periodic boundary conditions
n_neighbors = convolve2d(self.state, totalistic_kernel, mode='same', boundary='wrap')
next_state = np.logical_or(
np.logical_and(self.state == 1, n_neighbors == 2),
n_neighbors == 3
)
return next_state
model = GameOfLife(40, random_state=0)
model.simulate(200)
print("Average time per iteration: ", end="")
%timeit model.next_state()
plt.figure()
plt.imshow(model.initial_state, cmap="gray")
plt.title("Initial state")
plt.figure()
plt.imshow(model.state, cmap="gray")
plt.title("Final state")
Average time per iteration: 44 µs ± 357 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)
Text(0.5, 1.0, 'Final state')
Can we go even faster?¶
Based on the discrete convolution equation, we estimate that the runtime of a single update step is proportional to both the Moore neighborhood size and the total size of the grid, $\mathcal{O}(M\,N)$. Because the update rule is local, however, $M/N \ll 1$ and so the runtime is dominated by the factor $\mathcal{O}(N)$.
It turns out that we can speed up the convolution operation slightly by performing it in frequency space. Recall the equivalence between convolution and multiplication in the Fourier domain:
$$ \mathcal{F}\{f * g\} = \mathcal{F}\{f\} \cdot \mathcal{F}\{g\} $$
where $\mathcal{F}$ is the Fourier transform operator. This means that we can compute the convolution of two functions by multiplying their Fourier transforms and then taking the inverse Fourier transform of the result. Since the product is taken elementwise over the transformed array and kernel, we expect the product to require $\mathcal{O}(N)$ multiplication operations, but we gain the cost of performing and inverting the Fourier transform. The Fourier Transform integral in two dimensions has the form
$$ \mathcal{F}\{f(x,y)\} = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f(x,y) e^{-i\,k_x\,x} e^{-i\,k_y\,y} dx\,dy $$
where $k_x$ and $k_y$ are the spatial frequencies in the $x$ and $y$ directions, respectively. This integral takes $\mathcal{O}(N^2)$ operations to compute, but it turns out that we can compute the Fourier transform in $\mathcal{O}(N\,\log N)$ operations using the Fast Fourier Transform algorithm.
While the asymptotic (and thus intrinsic) complexities are the same for convolutions in real and frequency space, the Fourier approach often has smaller prefactors and enables parallelization. In fact, many 2D convolution functions in modern languages use the FFT under the hood to speed up the computation.
Extending the Game of Life to a continuous domain¶
- We saw that the update rule for the game of life is a function of the number of neighbors of a cell.
- We can replace the discrete neighbor kernel with a continuous one.
- This idea is extended and refined in Lenia, a fully continuous-time and continuous-space cellular automaton that exhibits surprisingly complex, life-like behavior. I highly encourage checking out the author's interactive notebook tutorial. The following code is a highly-simplified version capturing some ideas from Lenia.
from scipy.stats import multivariate_normal
from scipy.signal import convolve2d
# We can use type hints to indicate the types of the arguments and return values
# of functions. This is not enforced by Python, but it can be useful for
# documentation and for static type checkers like mypy.
def gaussian_kernel(size: int, sigma: float):
"""
Make a gaussian kernel of given size and standard deviation.
Args:
size (int): the size of the kernel
sigma (float): the standard deviation of the gaussian
Returns:
ndarray: the gaussian kernel
"""
kernel = np.fromfunction(
lambda x, y: (1/ (2 * np.pi * sigma**2)) *
np.exp(- ((x - (size - 1) / 2) ** 2 + (y - (size - 1) / 2) ** 2) / (2 * sigma**2)),
(size, size)
)
return kernel / np.sum(kernel)
def gaussian1d(x, mu, sigma):
return np.exp(-(x - mu)**2/(2 * sigma**2))/(sigma * np.sqrt(2 * np.pi))
class SmoothGameOfLife(CellularAutomaton):
"""
A class for the Game of Life on a continuous-valued lattice
"""
def __init__(self, n, r=1, **kwargs):
# the super method calls the parent class's __init__ method
# and passes the arguments to it. It is a constructor
super().__init__(n, 2, **kwargs)
self.r = r
def next_state(self):
"""
Compute the next state of the lattice using discrete convolutions to implement
the rules of the Game of Life.
"""
# Compute the next state
next_state = np.zeros_like(self.state)
# Make a gaussian convolutional kernel of 10 x 10 size
# with standard deviation 1
neighbors_kernel = gaussian_kernel(20, self.r)
neighbors_kernel /= np.max(neighbors_kernel)
self_kernel = gaussian_kernel(20, 2)
self_kernel /= np.max(self_kernel)
totalistic_kernel = neighbors_kernel - self_kernel
totalistic_kernel /= np.sum(totalistic_kernel)
# convolve with periodic boundary conditions
n_neighbors = convolve2d(self.state, totalistic_kernel, mode='same', boundary='wrap')
n_neighbors3 = gaussian1d(n_neighbors, 3, 0.5)
n_neighbors3 /= np.max(n_neighbors3)
n_neighbors2 = gaussian1d(n_neighbors, 2, 0.5)
n_neighbors2 /= np.max(n_neighbors2)
#next_state = gaussian1d(self.state * n_neighbors2 + n_neighbors3, 1, 0.5)
next_state = gaussian1d(self.state * n_neighbors2 + n_neighbors3, 1, 0.8)
next_state /= np.max(next_state)
return next_state
model = SmoothGameOfLife(800, 5, random_state=14)
model.simulate(500)
print("Average time per iteration: ", end="")
%timeit model.next_state()
plt.figure()
plt.imshow(model.initial_state, cmap="gray")
plt.title("Initial state")
plt.figure()
plt.imshow(model.state, cmap="gray")
plt.title("Final state")
Average time per iteration: 3.16 ms ± 26.7 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
Text(0.5, 1.0, 'Final state')
## Make an interactive video
from ipywidgets import interact, interactive, fixed, interact_manual, Layout
import ipywidgets as widgets
def plotter(i):
# plt.close()
fig = plt.figure(figsize=(6, 6))
plt.imshow(model.history[i], cmap="YlGnBu")
plt.show()
interact(
plotter,
i=widgets.IntSlider(0, 0, len(model.history) - 1, 1, layout=Layout(width='800px'))
)
interactive(children=(IntSlider(value=0, description='i', layout=Layout(width='800px'), max=500), Output()), _…
<function __main__.plotter(i)>
<function cphy.plotting.vanish_axes(gca=None)>
# This code saves the sliding frames as a movie
frames = model.history[4:]
vmin, vmax = np.percentile(frames, 1), np.percentile(frames, 99)
for i in range(len(frames) - 1):
out_path = "../private_dump/life/frame" + str(i).zfill(4) + ".png"
plt.figure()
plt.imshow(frames[i], vmin=vmin, vmax=vmax, cmap="YlGnBu")
ax = plt.gca()
ax.set_axis_off()
ax.xaxis.set_major_locator(plt.NullLocator())
ax.yaxis.set_major_locator(plt.NullLocator())
ax.set_aspect(1, adjustable='box')
cplot.vanish_axes()
plt.savefig(out_path, bbox_inches='tight', pad_inches=0.0, dpi=300)
plt.close()
Cellular automata¶
Cellular automata represent an intriguing example of a physical system where surprising complexity emerges from relatively simple rules.
Cellular automata have a rich history in the physics literature, although their usage has limitations due to a lack of analytical tools for working with their governing equations.
One common criticism of Game-of-Life like CAs comes from their requirement of a synchronous update rule: all of the cells update to a new state at the same time: the value of a cell $x_{ij}(t+1)$ depends on the "frozen" values $x_{ij}(t)$ and its neighbors at the same instant $t$. In a real-world system, updates are asynchronous: different cells will update at different times. One way of thinking about asynchronous updates is that we only update a single cell in a given timestep, and we gradually raster-scan across the entire image to update all cells sequentially. Asynchronous update rules usually have qualitatively distinct dynamics than synchronous CA.
Complexity of CA dynamics vs the entropy of the rule table. Image from Langton, 1990
- In addition to having interesting theoretical properties, there are a few natural systems that exhibit discrete, local space and time dynamics that can be well described by CA equations.
Scale patterning in lizards. Image from Manukyan et al. (2017)
Appendix¶
## Make an interactive video
from ipywidgets import interact, interactive, fixed, interact_manual, Layout
import ipywidgets as widgets
def plotter(i):
# plt.close()
fig = plt.figure(figsize=(6, 6))
plt.imshow(model.history[i], vmin=0, vmax=1, cmap="gray")
plt.show()
interact(
plotter,
i=widgets.IntSlider(0, 0, len(model.history) - 1, 1, layout=Layout(width='800px'))
)