Python Cellular Automata — Deep Dive

High-Performance Game of Life with NumPy

Convolution-Based Neighbor Counting

Instead of rolling arrays 8 times, use scipy.signal.convolve2d with a kernel that counts neighbors in one operation:

import numpy as np
from scipy.signal import convolve2d

KERNEL = np.array([[1, 1, 1],
                   [1, 0, 1],
                   [1, 1, 1]])

def step_convolve(grid):
    neighbors = convolve2d(grid, KERNEL, mode='same', boundary='wrap')
    birth = (grid == 0) & (neighbors == 3)
    survive = (grid == 1) & ((neighbors == 2) | (neighbors == 3))
    return (birth | survive).astype(np.uint8)

This is 3-5× faster than the np.roll approach because convolve2d runs in optimized C code.

Performance Comparison (1000×1000 grid, 100 steps)

Method	Time
Pure Python nested loops	~45 seconds
NumPy roll (8 rolls)	~2.5 seconds
scipy convolve2d	~0.8 seconds
Numba JIT	~0.3 seconds
CuPy (GPU)	~0.05 seconds

Wolfram Elementary CA Implementation

def elementary_ca(rule_number, width=201, steps=100):
    """Generate a 1D elementary cellular automaton."""
    # Decode rule number into lookup table
    rule_bin = format(rule_number, '08b')
    rules = {(int(rule_bin[7-i])): None for i in range(8)}
    for i in range(8):
        neighborhood = ((i >> 2) & 1, (i >> 1) & 1, i & 1)
        rules[neighborhood] = int(rule_bin[7 - i])

    # Initialize grid
    grid = np.zeros((steps, width), dtype=np.uint8)
    grid[0, width // 2] = 1  # single cell in center

    for t in range(1, steps):
        for x in range(width):
            left = grid[t-1, (x-1) % width]
            center = grid[t-1, x]
            right = grid[t-1, (x+1) % width]
            grid[t, x] = rules[(left, center, right)]

    return grid

# Vectorized version
def elementary_ca_fast(rule_number, width=1001, steps=500):
    rule_bin = np.array([int(b) for b in format(rule_number, '08b')][::-1], dtype=np.uint8)

    grid = np.zeros((steps, width), dtype=np.uint8)
    grid[0, width // 2] = 1

    for t in range(1, steps):
        left = np.roll(grid[t-1], 1)
        center = grid[t-1]
        right = np.roll(grid[t-1], -1)
        index = (left << 2) | (center << 1) | right
        grid[t] = rule_bin[index]

    return grid

Hashlife: Exponential Speedup

For very large grids evolved over millions of generations, the Hashlife algorithm (Bill Gosper, 1984) exploits the self-similar structure of Life patterns:

1. Quadtree representation — the grid is recursively divided into quadrants.
2. Memoization — identical subpatterns share the same node. A 256×256 region that appears in multiple places is stored once.
3. Temporal acceleration — a memoized node of size 2^n can compute 2^(n-2) steps in one operation.

from functools import lru_cache

class QuadNode:
    """Immutable quadtree node for Hashlife."""
    __slots__ = ('level', 'nw', 'ne', 'sw', 'se', 'population', '_hash')

    def __init__(self, level, nw=None, ne=None, sw=None, se=None, alive=None):
        self.level = level
        if level == 0:
            self.population = 1 if alive else 0
            self.nw = self.ne = self.sw = self.se = None
        else:
            self.nw, self.ne, self.sw, self.se = nw, ne, sw, se
            self.population = nw.population + ne.population + \
                              sw.population + se.population
        self._hash = hash((level, id(nw), id(ne), id(sw), id(se),
                          alive if level == 0 else None))

    def __hash__(self):
        return self._hash

    def __eq__(self, other):
        if self.level != other.level:
            return False
        if self.level == 0:
            return self.population == other.population
        return (self.nw is other.nw and self.ne is other.ne and
                self.sw is other.sw and self.se is other.se)

DEAD = QuadNode(0, alive=False)
ALIVE = QuadNode(0, alive=True)

@lru_cache(maxsize=1_000_000)
def make_node(level, nw, ne, sw, se):
    return QuadNode(level, nw, ne, sw, se)

@lru_cache(maxsize=1_000_000)
def advance(node):
    """Advance a node by 2^(level-2) generations."""
    if node.level < 2:
        raise ValueError("Cannot advance level < 2")

    if node.population == 0:
        return make_node(node.level - 1, DEAD, DEAD, DEAD, DEAD) if node.level > 2 \
               else QuadNode(0, alive=False)

    # ... recursive advance using 9 sub-quadrants
    # This is the core of Hashlife — memoized recursive computation
    pass

Hashlife can simulate the Gosper Glider Gun for 2^64 generations in seconds, because repeated patterns are cached and temporal jumps skip redundant computation.

GPU-Accelerated CA with CuPy

import cupy as cp

KERNEL_GPU = cp.array([[1, 1, 1],
                       [1, 0, 1],
                       [1, 1, 1]], dtype=cp.float32)

def step_gpu(grid_gpu):
    from cupyx.scipy.signal import convolve2d
    neighbors = convolve2d(grid_gpu, KERNEL_GPU, mode='same', boundary='wrap')
    birth = (grid_gpu == 0) & (neighbors == 3)
    survive = (grid_gpu == 1) & ((neighbors == 2) | (neighbors == 3))
    return (birth | survive).astype(cp.uint8)

# Initialize on GPU
size = 4096
grid_gpu = cp.random.randint(0, 2, (size, size), dtype=cp.uint8)

# Run 1000 steps
for _ in range(1000):
    grid_gpu = step_gpu(grid_gpu)

# Transfer back to CPU for visualization
result = cp.asnumpy(grid_gpu)

For a 4096×4096 grid, GPU execution is 50-100× faster than CPU NumPy.

Custom Rule Systems

Larger-than-Life

Extends Life’s rules to larger neighborhoods. Instead of 3×3 (range 1), use range R with customizable birth and survival thresholds:

def larger_than_life_step(grid, R=5, birth_min=34, birth_max=45,
                          survive_min=34, survive_max=58):
    # Build circular kernel of radius R
    y, x = np.ogrid[-R:R+1, -R:R+1]
    kernel = (x*x + y*y <= R*R).astype(int)
    kernel[R, R] = 0  # exclude center

    neighbors = convolve2d(grid, kernel, mode='same', boundary='wrap')
    birth = (grid == 0) & (neighbors >= birth_min) & (neighbors <= birth_max)
    survive = (grid == 1) & (neighbors >= survive_min) & (neighbors <= survive_max)
    return (birth | survive).astype(np.uint8)

Larger-than-Life produces smooth, blob-like patterns useful for procedural terrain and organic textures.

Continuous CAs (Lenia)

Lenia (Bert Chan, 2018) generalizes Life to continuous space and states:

def lenia_step(grid, kernel, growth_func, dt=0.1):
    """One step of a Lenia simulation."""
    # Convolution gives neighborhood potential
    from scipy.signal import fftconvolve
    potential = fftconvolve(grid, kernel, mode='same')

    # Growth function maps potential to growth rate
    growth = growth_func(potential)

    # Continuous state update
    return np.clip(grid + dt * growth, 0, 1)

def gaussian_growth(x, mu=0.15, sigma=0.015):
    return 2 * np.exp(-((x - mu) ** 2) / (2 * sigma ** 2)) - 1

# Ring-shaped kernel
def ring_kernel(R=13, r_inner=0.5):
    y, x = np.ogrid[-R:R+1, -R:R+1]
    dist = np.sqrt(x*x + y*y) / R
    kernel = np.exp(-((dist - r_inner) ** 2) / 0.1)
    kernel[(dist > 1)] = 0
    return kernel / kernel.sum()

Lenia produces lifelike creatures that move, divide, and interact — blurring the line between simulation and artificial life.

Practical Applications

Procedural Cave Generation

def generate_cave(width, height, fill_prob=0.45, steps=5):
    # Random initialization
    grid = (np.random.random((height, width)) < fill_prob).astype(np.uint8)

    # Border walls
    grid[0, :] = grid[-1, :] = grid[:, 0] = grid[:, -1] = 1

    for _ in range(steps):
        neighbors = convolve2d(grid, np.ones((3,3)), mode='same', boundary='fill')
        # 4-5 rule: wall if 5+ neighbors are walls, or fewer than 3 are walls
        grid = ((neighbors >= 5) | (neighbors <= 2)).astype(np.uint8)
        grid[0, :] = grid[-1, :] = grid[:, 0] = grid[:, -1] = 1

    return grid

After 4-5 iterations, the random noise resolves into organic cave systems with smooth walls and connected chambers.

Forest Fire Model

EMPTY, TREE, FIRE = 0, 1, 2

def forest_fire_step(grid, p_grow=0.01, p_ignite=0.0001):
    new_grid = grid.copy()
    rows, cols = grid.shape

    # Trees adjacent to fire catch fire
    fire_neighbors = convolve2d(
        (grid == FIRE).astype(int),
        np.ones((3,3)), mode='same', boundary='fill'
    )
    new_grid[(grid == TREE) & (fire_neighbors > 0)] = FIRE

    # Random lightning strikes
    lightning = np.random.random((rows, cols)) < p_ignite
    new_grid[(grid == TREE) & lightning] = FIRE

    # Fire burns out
    new_grid[grid == FIRE] = EMPTY

    # Empty cells grow trees
    growth = np.random.random((rows, cols)) < p_grow
    new_grid[(grid == EMPTY) & growth] = TREE

    return new_grid

This model exhibits self-organized criticality — the system naturally reaches a state where fires of all sizes occur with a power-law distribution, matching real wildfire statistics.

Visualization

Matplotlib Animation

import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

fig, ax = plt.subplots()
grid = np.random.randint(0, 2, (200, 200), dtype=np.uint8)
img = ax.imshow(grid, cmap='binary', interpolation='nearest')

def animate(frame):
    global grid
    grid = step_convolve(grid)
    img.set_data(grid)
    return [img]

anim = FuncAnimation(fig, animate, frames=500, interval=50, blit=True)
plt.show()

Pygame Real-Time Display

For interactive exploration (zoom, pan, place cells), Pygame provides a direct pixel buffer that updates at 60+ FPS for grids up to ~2000×2000.

The one thing to remember: Cellular automata leverage simple local rules to generate complex emergent behavior — and Python’s NumPy, SciPy, and GPU libraries make it possible to simulate millions of cells in real time, powering everything from cave generators to artificial life research.