Chapter 20: The Performance Problem — and the Solution¶
Part III — Foundations of Spatial Simulation
Learning Objectives¶
By the end of this chapter you will be able to:
- Profile Python code with
cProfileinstead of guessing where time goes - Explain why a per-cell Python loop is fundamentally slow, not just unlucky
- Vectorize a cellular automaton rule with NumPy array operations
- Measure the resulting speedup, and validate that the fast version computes the same answer as the slow one
# Standard imports — add chapter-specific imports below
import random
import time
import cProfile
import pstats
import numpy as np
Chapter 18 ended with a number and a promise: a 100×100 Game of Life grid took over a hundred milliseconds for 20 generations, and this chapter would explain why, then fix it. Chapter 19 added scheduling flexibility but left that number completely untouched — salabim changes when code runs, not how fast the code itself is. This chapter is where the actual computation gets faster, and — just as importantly — where "faster" gets proven, not assumed.
Profiling Tools: cProfile and line_profiler¶
Guessing where time goes is exactly the trap Chapter 18's final exercise set up on purpose — most guesses about performance are wrong, including experienced programmers' guesses, which is the entire reason profiling tools exist. cProfile, in the standard library, needs no installation and answers "which function ate the time" precisely:
def count_live_neighbors(grid, r, c, rows, cols):
count = 0
for dr in (-1, 0, 1):
for dc in (-1, 0, 1):
if dr == 0 and dc == 0:
continue
nr, nc = r + dr, c + dc
if 0 <= nr < rows and 0 <= nc < cols:
count += grid[nr][nc]
return count
def step(grid, rows, cols):
next_grid = [[0] * cols for _ in range(rows)]
for r in range(rows):
for c in range(cols):
alive = grid[r][c] == 1
n = count_live_neighbors(grid, r, c, rows, cols)
if alive and n in (2, 3):
next_grid[r][c] = 1
elif not alive and n == 3:
next_grid[r][c] = 1
return next_grid
line_profiler (pip install line_profiler, not in the standard library) goes one level deeper than cProfile — timing individual lines inside a function, not just whole function calls. cProfile alone is already enough to answer this chapter's central question, so that's where the measurement below starts.
Measuring the Naive CA¶
Wrap Chapter 18's 20-generation run in cProfile and let the profiler answer, precisely, the question Chapter 18's last exercise only asked you to guess at:
random.seed(0)
rows = cols = 100
grid = [[1 if random.random() < 0.3 else 0 for _ in range(cols)] for _ in range(rows)]
def run_naive():
g = grid
for _ in range(20):
g = step(g, rows, cols)
profiler = cProfile.Profile()
profiler.enable()
run_naive()
profiler.disable()
stats = pstats.Stats(profiler).sort_stats("cumulative")
stats.print_stats(4)
200415 function calls (200413 primitive calls) in 0.465 seconds
Ordered by: cumulative time
List reduced from 110 to 4 due to restriction <4>
ncalls tottime percall cumtime percall filename:lineno(function)
20 0.069 0.003 0.425 0.021 /tmp/ipykernel_404817/1922488936.py:13(step)
200000 0.389 0.000 0.389 0.000 /tmp/ipykernel_404817/1922488936.py:1(count_live_neighbors)
1 0.000 0.000 0.315 0.315 /tmp/ipykernel_404817/934666796.py:5(run_naive)
1 0.000 0.000 0.135 0.135 /home/sergio/develop/github/lambdageo/ebooks/geospatial-modeling-python/.venv/lib/python3.12/site-packages/IPython/core/history.py:157(wrapper)
<pstats.Stats at 0x7f7118093140>
count_live_neighbors alone accounts for roughly 80% of the total time — 200,000 calls (20 generations × 10,000 cells), each one a small Python function call plus an 8-iteration inner loop, each iteration a bounds check and a list index. step's own outer double loop, by contrast, barely registers. If you guessed "the neighbor-counting function," Chapter 18's exercise was pointing you toward exactly the right place.
Why Python Loops Are Slow¶
Every one of those 200,000 count_live_neighbors calls pays the same fixed costs, over and over: a Python function call has real overhead (building a new frame, resolving grid, r, c as names each time); every grid[nr][nc] is a full Python object lookup — a list of lists stores boxed integer objects, not raw machine integers, so every access means dereferencing a pointer, checking a type, then finally reading a value; every +=, every comparison, is a full Python bytecode dispatch, not a single CPU instruction.
None of this is a bug, and none of it is unique to this particular function — it's the fixed cost of any Python-level loop over any large collection, the identical tax paid by Chapter 5's early loops, and precisely the tax Chapter 2's very first vectorization example (computing distances between cities with NumPy instead of a loop) was already quietly avoiding, six chapters before this one named the problem directly. NumPy's answer is to stop asking Python to loop at all — push the entire loop down into compiled C code, operating on one large, contiguous block of raw numbers instead of ten thousand individually-boxed Python objects.
Vectorization with NumPy¶
The trick that makes an 8-neighbor sum vectorizable: instead of asking, for each cell, "what are my neighbors' values," ask the opposite question for the whole grid at once — "if I shift the entire grid one step in each of the 8 directions, what value now sits where each cell used to be." Summing those 8 shifted copies gives every cell's neighbor count in one pass, no per-cell loop anywhere:
def step_numpy(grid):
padded = np.pad(grid, 1, mode="constant", constant_values=0)
neighbor_count = np.zeros_like(grid)
for dr in (-1, 0, 1):
for dc in (-1, 0, 1):
if dr == 0 and dc == 0:
continue
neighbor_count += padded[1 + dr : 1 + dr + grid.shape[0], 1 + dc : 1 + dc + grid.shape[1]]
alive = grid == 1
survive = alive & ((neighbor_count == 2) | (neighbor_count == 3))
born = (~alive) & (neighbor_count == 3)
return (survive | born).astype(grid.dtype)
The outer for dr / for dc loop still runs in Python — but only 8 times total, once per neighbor direction, regardless of whether the grid is 100×100 or 10,000×10,000. np.pad handles the edge case Chapter 18's manual bounds-checking (0 <= nr < rows) did by hand: padding the grid with a border of zeros means a shifted slice never actually runs off the array, so no per-cell edge check is needed at all. Everything inside the loop — the slicing, the addition — runs as a single compiled operation over the entire array at once, the same category of operation Chapter 8's raster algebra and Chapter 16's distance_transform_edt both leaned on.
Before trusting the speed, confirm the two versions agree — a vectorized rewrite that's fast but wrong is worse than the slow original, exactly the validation discipline Chapter 21's ModelExecutor argued for in a different context:
grid_np = np.array(grid, dtype=np.int8)
g_naive, g_fast = grid, grid_np.copy()
for gen in range(5):
g_naive = step(g_naive, rows, cols)
g_fast = step_numpy(g_fast)
match = np.array_equal(np.array(g_naive), g_fast)
print(f"generation {gen}: matches naive version = {match}")
generation 0: matches naive version = True generation 1: matches naive version = True generation 2: matches naive version = True generation 3: matches naive version = True generation 4: matches naive version = True
Benchmarking: Loop vs Array¶
With correctness confirmed, the timing comparison finally means something:
random.seed(0)
grid_list = [[1 if random.random() < 0.3 else 0 for _ in range(cols)] for _ in range(rows)]
grid_np = np.array(grid_list, dtype=np.int8)
start = time.perf_counter()
g = grid_list
for _ in range(20):
g = step(g, rows, cols)
naive_ms = (time.perf_counter() - start) * 1000
start = time.perf_counter()
g = grid_np
for _ in range(20):
g = step_numpy(g)
vectorized_ms = (time.perf_counter() - start) * 1000
print(f"Naive Python loop: {naive_ms:8.2f} ms")
print(f"Vectorized NumPy: {vectorized_ms:8.2f} ms")
print(f"Speedup: {naive_ms / vectorized_ms:.0f}x")
Naive Python loop: 273.85 ms Vectorized NumPy: 1.72 ms Speedup: 159x
On a 100×100 grid, the vectorized version comes in well over a hundred times faster — and the gap only widens as the grid grows, since the naive version's cost scales with the number of cells, while the vectorized version's fixed 8-direction Python loop stays exactly 8 iterations no matter how large the array gets underneath it. This is precisely the number Chapter 24 later measures again, on FireModel specifically, where the same vector-versus-raster gap shows up at a similar order of magnitude in DisSModel's own real code — this chapter's toy benchmark and that chapter's production one are the identical phenomenon.
What We Still Need¶
Vectorizing step() solved the arithmetic — counting neighbors and applying Conway's rule now runs at NumPy speed. It solved nothing else. The grid is still a bare NumPy array with no notion of geographic coordinates, no CRS, no way to represent an irregular boundary the way Chapter 7's GeoDataFrame can. Scheduling still means a manual Python for loop calling step() repeatedly — none of Chapter 19's salabim flexibility (multiple independent components, irregular timing) carries over automatically. And every one of Chapters 18 through 19's three separate concerns — the rule itself, the neighbor mechanics, and the scheduler driving it forward — still lives in three disconnected pieces of code with no shared contract between them.
Chapter 21 is where all three merge into one framework: the scheduler Chapter 19 explored, the vectorized substrate this chapter built, and a rule(idx) contract that lets the exact same Game of Life logic run on either a GeoDataFrame or a raw NumPy array, with no manual wiring between them at all.
Exercises¶
- Profile the vectorized version. Run
cProfileonstep_numpythe same way Measuring the Naive CA did forstep. Where does the (much smaller) remaining time go? - Scale the grid. Rerun Benchmarking: Loop vs Array at
rows = cols = 200instead of 100. Does the speedup ratio grow, shrink, or stay about the same? Explain what you observe in terms of how each version's cost scales with cell count. - A padding bug, on purpose. Change
mode="constant", constant_values=0tomode="wrap"instep_numpyand rerun the correctness check againststep(). Does it still match? What doesmode="wrap"represent physically, and which chapter's TerraME/DisSModel comparison already named this exact difference? - line_profiler, for real. Install
line_profiler(pip install line_profiler) and use its%lprunmagic (or the@profiledecorator with thekernprofcommand-line tool) oncount_live_neighbors. Which single line inside it costs the most?
# Your code here
Summary¶
Key concepts introduced¶
cProfile, measuring precisely rather than guessing where time goes — confirmingcount_live_neighborsas roughly 80% of the naive implementation's runtime- Why a Python-level loop is fundamentally slow: per-call overhead, boxed objects, bytecode dispatch, paid by every large Python loop, not a flaw specific to this function
- Vectorization: replacing "for each cell, check its neighbors" with "shift the whole grid 8 times and sum," pushing the loop into compiled NumPy code
- Validating a vectorized rewrite against the original before trusting its speed — a fast, silently-wrong function is worse than a slow, correct one
- A measured, order-of-magnitude speedup (roughly 150x on a 100×100 grid), foreshadowing the real
FireModelbenchmark Chapter 24 measures later - What vectorization alone doesn't solve: geographic coordinates, flexible scheduling, and a shared contract across rule/substrate/scheduler — the three things Chapter 21's framework provides together
This closes Part III. Every technique from here — the CA rule, the scheduler, the vectorized substrate — reappears inside DisSModel starting next chapter, this time built once, rather than reinvented per model.
Further Reading¶
- Python documentation, The Python Profilers: https://docs.python.org/3/library/profile.html
line_profileron GitHub: https://github.com/pyutils/line_profiler- NumPy documentation, NumPy fundamentals — the array model underlying every vectorization technique in this chapter: https://numpy.org/doc/stable/user/basics.html
- VanderPlas, J. (2016). Python Data Science Handbook, Chapter 2 ("Introduction to NumPy"). O'Reilly