L01 · Wildfire Percolation

"""L01 - Wildfire Percolation.

Part A: a fire lit on the left edge of a random forest (tree density p) spreads to the 4-neighbourhood.
        Monte Carlo over p and system size L shows the percolation phase transition at p_c = 0.5927
        and its finite-size scaling.
Part B: the Drossel-Schwabl forest-fire model (trees grow, lightning strikes) drives itself to a
        critical state: fire sizes follow a power law without anyone tuning p.

NumPy + Matplotlib only.  Run from the repo root:  uv run python L01_wildfire_percolation/sim.py
"""
import os
os.environ.setdefault("OMP_NUM_THREADS", "4")  # CPU etiquette (NumPy here is single-threaded anyway)

import time
from pathlib import Path

import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation, PillowWriter

OUT = Path(__file__).parent / "output"

# Literature values we verify against
P_C = 0.592746    # square-lattice site percolation threshold (Newman & Ziff 2000: 0.59274621)
NU = 4 / 3        # correlation-length exponent, exact in 2D (den Nijs 1979, Nienhuis 1982)
D_MIN = 1.1306    # shortest-path (chemical) fractal dimension (Zhou, Yang, Deng & Ziff 2012)
RHO_DS = 0.397636  # Drossel-Schwabl mean tree density at p/f = 1000 (Grassberger 2002, Table 1)
TAU_FIRE = (1.11, 1.19)  # DS fire-size exponent = cluster exponent (2.11-2.19, Clar et al. 1996; Grassberger 2002) - 1

# Dark theme (chrome from the dataviz reference palette); ordinal fire ramp: small L dim -> large L bright
BG, INK, INK2, MUTED, GRID = "#1a1a19", "#ffffff", "#c3c2b7", "#898781", "#2c2c2a"
L_COLORS = ["#b04a17", "#de6a26", "#f59a45", "#fcd08a"]
plt.rcParams.update({
    "figure.facecolor": BG, "axes.facecolor": BG, "savefig.facecolor": BG,
    "axes.edgecolor": "#383835", "axes.labelcolor": INK2, "text.color": INK, "axes.titlecolor": INK,
    "xtick.color": MUTED, "ytick.color": MUTED, "axes.grid": True, "grid.color": GRID, "grid.linewidth": 0.7,
    "axes.spines.top": False, "axes.spines.right": False, "legend.frameon": False,
    "font.size": 11, "axes.titlesize": 12.5, "axes.titleweight": "bold", "axes.titlelocation": "left",
    "lines.linewidth": 2, "lines.markersize": 5,
})
# Cell colours for the fire maps
EMPTY, TREE = np.array([0.05, 0.05, 0.06]), np.array([0.16, 0.36, 0.22])
FLAME = np.array([[1.0, 0.97, 0.80], [1.0, 0.80, 0.30], [1.0, 0.52, 0.12], [0.86, 0.26, 0.09],
                  [0.58, 0.16, 0.08], [0.40, 0.15, 0.10], [0.28, 0.14, 0.11]])  # age 0 -> cold ash
AGE_COLOURS = np.array([[0.60, 0.79, 0.38], [0.30, 0.56, 0.27], [0.12, 0.34, 0.19]])  # young / mid / old trees


# ----------------------------------------------------------------------------------------------------
# The fire itself: one cellular-automaton rule, simulated front by front
# ----------------------------------------------------------------------------------------------------
def spread(trees, front, neighbours):
    """Yield the burning front at t = 0, 1, 2, ...  CA rule: a tree ignites if a 4-neighbour is
    burning; a burning cell becomes ash next step.  Only front cells can change, so we update those
    instead of all cells (same automaton, O(front) instead of O(L^2) per step).
    `trees` is a flat bool array that the fire consumes in place (burned tree -> not a tree)."""
    front = front[trees[front]]
    trees[front] = False
    while front.size:
        yield front
        nb = neighbours(front)
        front = np.unique(nb[trees[nb]])
        trees[front] = False


def arrival_times(trees):
    """Ignite column 0 of each forest (bool array (H, W) or (S, H, W)); return the step at which each
    cell caught fire, -1 = never.  Open boundaries (a border of empty cells stops the fire)."""
    shape = trees.shape
    pad = np.pad(trees.reshape(-1, *shape[-2:]), ((0, 0), (1, 1), (1, 1)))
    S, Hp, Wp = pad.shape
    seeds = (np.arange(S)[:, None] * Hp * Wp + np.arange(1, Hp - 1) * Wp + 1).ravel()
    t_arr = np.full(pad.size, -1, np.int32)
    for t, front in enumerate(spread(pad.ravel(), seeds,
                                     lambda f: np.concatenate((f - 1, f + 1, f - Wp, f + Wp)))):
        t_arr[front] = t
    return t_arr.reshape(S, Hp, Wp)[:, 1:-1, 1:-1].reshape(shape)


def arrival_times_dense(trees):
    """Textbook version of the same CA: update every cell every step.  Slow but obviously correct;
    verify() checks that arrival_times() reproduces it exactly."""
    burning = trees.copy()
    burning[:, 1:] = False
    unburnt = trees & ~burning
    t_arr = np.where(burning, 0, -1)
    t = 0
    while burning.any():
        t += 1
        nb = np.zeros_like(burning)
        nb[1:] |= burning[:-1]; nb[:-1] |= burning[1:]; nb[:, 1:] |= burning[:, :-1]; nb[:, :-1] |= burning[:, 1:]
        burning = nb & unburnt
        unburnt &= ~burning
        t_arr[burning] = t
    return t_arr


# ----------------------------------------------------------------------------------------------------
# Part A: Monte Carlo sweep over density p and system size L
# ----------------------------------------------------------------------------------------------------
def sweep(Ls=(32, 64, 128, 256), samples=(2500, 1200, 500, 200), seed=1):
    """For each L draw S forests ONCE as uniform numbers U (tree where U < p) and reuse them at every p
    (common random numbers: curves are smooth and monotone, differences between p are pure signal)."""
    ps = np.unique(np.round(np.r_[np.arange(0.40, 0.801, 0.025), np.arange(0.55, 0.6401, 0.003)], 4))
    rng = np.random.default_rng(seed)
    out = {}
    for L, S in zip(Ls, samples):
        t0 = time.perf_counter()
        U = rng.random((S, L, L), dtype=np.float32)
        cross, burned, dur = [], [], []
        for p in ps:
            trees = U < p
            t_arr = arrival_times(trees)
            cross.append((t_arr[:, :, -1] >= 0).any(1))
            burned.append((t_arr >= 0).sum((1, 2)) / np.maximum(trees.sum((1, 2)), 1))
            dur.append(t_arr.max((1, 2)))
        out[L] = dict(cross=np.array(cross), burned=np.array(burned), dur=np.array(dur), S=S)
        print(f"  L={L:4d}  S={S:5d} forests x {len(ps)} densities  {time.perf_counter() - t0:6.1f} s")
    return ps, out


def logistic_fit(ps, P, n):
    """Weighted least squares on logit(P) = (p - p_half) / w  ->  (p_half, w).  No SciPy needed."""
    m = (P > 0.02) & (P < 0.98)
    a, b = np.polyfit(ps[m], np.log(P[m] / (1 - P[m])), 1, w=np.sqrt(n * P[m] * (1 - P[m])))
    return -b / a, 1 / a


def estimate(ps, res, pick):
    """Finite-size-scaling analysis using forests pick[L] of each size."""
    Ls = np.array(sorted(res))
    p_half, width = np.array([logistic_fit(ps, res[L]["cross"][:, pick[L]].mean(1), len(pick[L])) for L in Ls]).T
    i_c = np.argmin(abs(ps - P_C))
    dur_c = np.array([res[L]["dur"][i_c, pick[L]].mean() for L in Ls])
    lin = np.polyfit(Ls ** (-1 / NU), p_half, 1)                     # p_half = p_c + a L^(-1/nu)
    return dict(p_c=lin[1], lin=lin, p_half=p_half, width=width, dur_c=dur_c,
                nu=-1 / np.polyfit(np.log(Ls), np.log(width), 1)[0],  # width ~ L^(-1/nu)
                d_min=np.polyfit(np.log(Ls), np.log(dur_c), 1)[0])    # duration at p_c ~ L^d_min


def analyse(ps, res, n_boot=300, seed=3):
    """Point estimates + bootstrap standard errors: resample the forests with replacement and redo the
    whole analysis.  (Resampling whole forests keeps the correlation that common random numbers
    create between densities, which a naive per-point error bar would ignore.)"""
    rng = np.random.default_rng(seed)
    a = estimate(ps, res, {L: np.arange(r["S"]) for L, r in res.items()})
    boot = [estimate(ps, res, {L: rng.integers(r["S"], size=r["S"]) for L, r in res.items()})
            for _ in range(n_boot)]
    a["err"] = {k: np.std([b[k] for b in boot], axis=0) for k in ("p_c", "nu", "d_min", "p_half", "width", "dur_c")}
    a["Ls"] = np.array(sorted(res))
    a["P_at_pc"] = {L: np.interp(P_C, ps, res[L]["cross"].mean(1)) for L in a["Ls"]}
    return a


# ----------------------------------------------------------------------------------------------------
# Part B: Drossel-Schwabl forest-fire model (self-organised criticality)
# ----------------------------------------------------------------------------------------------------
def drossel_schwabl(L, theta, strikes, warmup, frame_every=0, seed=7):
    """theta = p/f -> large limit (Drossel & Schwabl 1992; algorithm as in Grassberger 1993):
    between two lightning strikes theta random sites try to grow a tree (empty -> tree), then lightning
    hits one random site and, if it holds a tree, that whole connected cluster burns instantly.
    Periodic boundaries.  Returns fire sizes after warm-up, density per strike, optional RGB frames."""
    rng = np.random.default_rng(seed)
    N = L * L
    idx = np.arange(N).reshape(L, L)
    nbr = np.stack([np.roll(idx, s, a) for s in (1, -1) for a in (0, 1)], -1).reshape(N, 4)
    neighbours = lambda f: nbr[f].ravel()
    forest = np.zeros(N, bool)
    planted = np.zeros(N, np.int32)          # strike index when each tree grew (for its age colour)
    burnt = np.full(N, -10**8, np.int32)     # strike index when each site last burned
    sizes, rho, frames = np.zeros(strikes, np.int64), np.zeros(strikes), []
    for k in range(strikes):
        sites = rng.integers(N, size=theta)
        new = sites[~forest[sites]]
        forest[new] = True
        planted[new] = k
        rho[k] = np.count_nonzero(forest) / N
        for front in spread(forest, rng.integers(N, size=1), neighbours):
            burnt[front] = k
            sizes[k] += front.size
        if frame_every and k >= warmup and (k - warmup) % frame_every == 0:
            frames.append(ds_rgb(forest, planted, burnt, k, L, theta, frame_every))
    frames.append(ds_rgb(forest, planted, burnt, k, L, theta, frame_every or N // (40 * theta)))  # final state
    return sizes[warmup:], rho, frames


def ds_rgb(forest, planted, burnt, k, L, theta, every):
    """Trees in 3 age classes (few colours keep the GIF small); fresh burns glow and cool to soil."""
    regrowths = (k - planted) * theta / (L * L)          # age in units of the forest regrowth time N/theta
    trees = AGE_COLOURS[np.searchsorted([0.3, 0.8], regrowths, side="right")]
    fire = FLAME[np.clip((k - burnt) // every, 0, len(FLAME) - 1)]
    fire[k - burnt >= every * (len(FLAME) - 1)] = EMPTY + 0.03
    return (np.where(forest[:, None], trees, fire).reshape(L, L, 3) * 255).astype(np.uint8)


def power_law(sizes, s_lo=10, s_hi=2000):
    """Log-binned fire-size distribution n(s) (a raw histogram is useless on log-log: most big-s bins
    hold 0 or 1 fire) and a least-squares slope over the scaling range s_lo..s_hi."""
    s = sizes[sizes > 0]
    edges = np.unique(np.logspace(0, np.log10(s.max() + 1), 45).astype(int))
    counts = np.histogram(s, edges)[0]
    keep = counts > 0
    centres = np.sqrt(edges[:-1] * edges[1:])[keep]
    n = (counts / np.diff(edges) / s.size)[keep]
    fit = (centres >= s_lo) & (centres <= s_hi)
    slope, icpt = np.polyfit(np.log(centres[fit]), np.log(n[fit]), 1)
    return dict(s=centres, n=n, fit=fit, tau=-slope, icpt=icpt, fires=s.size)


# ----------------------------------------------------------------------------------------------------
# Figures and animations
# ----------------------------------------------------------------------------------------------------
def fire_rgb(trees, t_arr, t, step=1):
    """RGB frame of a percolation fire at time t: empty / unburned tree / flame cooling into ash."""
    age = t - t_arr
    rgb = np.where(trees[..., None], TREE, EMPTY)
    lit = (t_arr >= 0) & (age >= 0)
    rgb[lit] = FLAME[np.minimum(age[lit] // step, len(FLAME) - 1)]
    return rgb


def hero(path, H=450, W=800):
    """Arrival-time map of a fire that just barely crosses a forest at p = p_c."""
    for seed in range(100):
        trees = np.random.default_rng(seed).random((H, W)) < P_C
        t_arr = arrival_times(trees)
        if (t_arr[:, -1] >= 0).any():
            break
    rgb = np.where(trees[..., None], TREE * 0.6, EMPTY * 0.7)
    lit = t_arr >= 0
    rgb[lit] = plt.cm.inferno(0.2 + 0.8 * t_arr[lit] / t_arr.max())[:, :3]
    fig = plt.figure(figsize=(16, 9), dpi=100)
    fig.add_axes([0, 0, 1, 1]).imshow(rgb, interpolation="nearest", aspect="auto")
    fig.axes[0].set_axis_off()
    fig.text(0.025, 0.93, "Wildfire at the percolation threshold", fontsize=30, weight="bold",
             bbox=dict(fc=BG, ec="none", alpha=0.8, pad=8))
    fig.text(0.025, 0.875, f"{W} x {H} random forest, tree density p = p_c = 0.5927.  The fire is lit on the left "
             f"edge; colour = time it arrived (purple early, yellow after {t_arr.max()} steps).",
             fontsize=13, color=INK2, bbox=dict(fc=BG, ec="none", alpha=0.8, pad=6))
    fig.savefig(path)
    plt.close(fig)
    return seed, int(t_arr.max()), lit.sum() / trees.sum()


def gif_fire(path, L=160, ps=(0.55, P_C, 0.65), frames=150):
    """Three forests drawn from the SAME random numbers, below / at / above p_c, burning in lockstep."""
    for seed in range(100):  # pick a forest whose critical fire makes it across (P ~ 1/2 at p_c)
        U = np.random.default_rng(seed).random((L, L))
        if (arrival_times(U < P_C)[:, -1] >= 0).any():
            break
    forests = [U < p for p in ps]
    arrivals = [arrival_times(tr) for tr in forests]
    T = max(a.max() for a in arrivals) + 1
    step = max(1, T // frames)
    times = list(range(0, T + 12 * step, step))
    px, gap, top, bot = 2 * L, 24, 70, 40                     # 2 screen px per cell, exact
    Wfig, Hfig = 3 * px + 4 * gap, px + top + bot
    fig = plt.figure(figsize=(Wfig / 100, Hfig / 100), dpi=100)
    labels = ["sub-critical", "critical", "super-critical"]
    ims, txts = [], []
    for i, (tr, a, p, lab) in enumerate(zip(forests, arrivals, ps, labels)):
        ax = fig.add_axes([(gap + i * (px + gap)) / Wfig, bot / Hfig, px / Wfig, px / Hfig])
        ims.append(ax.imshow(fire_rgb(tr, a, 0), interpolation="nearest"))
        ax.set_axis_off()
        ax.set_title(f"p = {p:.4f}  {lab}", fontsize=12, pad=6, loc="center")
        txts.append(fig.text((gap + i * (px + gap) + px / 2) / Wfig, 12 / Hfig, "", ha="center",
                             fontsize=11, color=INK2))
    fig.text(0.5, 1 - 22 / Hfig, "Fire lit on the left edge of three random forests (same seed, different density)",
             ha="center", fontsize=13, weight="bold")

    def draw(t):
        for im, tx, tr, a in zip(ims, txts, forests, arrivals):
            im.set_data(fire_rgb(tr, a, t, step))
            tx.set_text(f"t = {min(t, a.max()):4d}   burned {((a >= 0) & (a <= t)).sum() / tr.sum():5.1%}")
        return ims + txts

    FuncAnimation(fig, draw, frames=times, blit=False).save(path, writer=PillowWriter(fps=20))
    plt.close(fig)


def gif_ds(path, frames, rho, warmup, every, L, theta):
    """Forest map (top) and the tree density trace (bottom) of a Drossel-Schwabl run."""
    px = 2 * L
    Wfig, Hfig = px + 40, px + 250
    fig = plt.figure(figsize=(Wfig / 100, Hfig / 100), dpi=100)
    ax = fig.add_axes([20 / Wfig, 210 / Hfig, px / Wfig, px / Hfig])
    im = ax.imshow(frames[0], interpolation="nearest")
    ax.set_axis_off()
    ax.set_title(f"Drossel-Schwabl forest, L = {L}, p/f = {theta}", fontsize=11, pad=6, loc="center")
    ax2 = fig.add_axes([70 / Wfig, 55 / Hfig, (px - 70) / Wfig, 120 / Hfig])
    k = warmup + every * np.arange(len(frames))
    sweeps = (k - warmup) * theta / L**2
    ax2.axhline(P_C, color=MUTED, ls="--", lw=1)
    ax2.text(0.01 * sweeps[-1], P_C - 0.012, "percolation threshold $p_c$ (never reached)", va="top",
             fontsize=9, color=MUTED)
    (line,) = ax2.plot([], [], color=AGE_COLOURS[0], lw=1.5)
    ax2.set(xlim=(0, sweeps[-1]), ylim=(0.25, 0.62), xlabel="time (growth attempts per site)", ylabel="tree density")

    def draw(i):
        im.set_data(frames[i])
        line.set_data(sweeps[:i + 1], rho[k[:i + 1]])
        return im, line

    FuncAnimation(fig, draw, frames=len(frames), blit=False).save(path, writer=PillowWriter(fps=15))
    plt.close(fig)


def log_ticks(ax, xt, yt):
    """Log-log axes labelled with plain numbers at the given ticks (instead of crowded 10^x minors)."""
    ax.set(xscale="log", yscale="log")
    ax.minorticks_off()
    ax.set_xticks(xt, [f"{v:g}" for v in xt])
    ax.set_yticks(yt, [f"{v:g}" for v in yt])


def plot_sweep(path, ps, res):
    fig, axes = plt.subplots(1, 3, figsize=(16, 4.8), constrained_layout=True)
    for L, c in zip(sorted(res), L_COLORS):
        for ax, key in zip(axes, ("cross", "burned", "dur")):
            ax.plot(ps, res[L][key].mean(1), "-o", color=c, ms=3, label=f"L = {L}")
    for ax in axes:
        ax.axvline(P_C, color=MUTED, ls="--", lw=1)
        ax.set_xlim(0.40, 0.80)
        ax.set_xlabel("tree density $p$")
    axes[0].text(P_C + 0.005, 0.03, "$p_c = 0.5927$", color=MUTED, fontsize=10)
    axes[0].set(ylabel="P(fire crosses the forest)", title="A  Crossing probability sharpens with L")
    axes[1].set(ylabel="fraction of trees burned", title="B  Fraction of the forest burned")
    axes[2].set(ylabel="fire duration (CA steps)", yscale="log", title="C  Critical slowing down")
    for ax, loc in zip(axes, ("upper left", "upper left", "lower right")):
        ax.legend(loc=loc)
    fig.savefig(path, dpi=150)
    plt.close(fig)


def plot_fss(path, ps, res, a):
    fig, axes = plt.subplots(1, 4, figsize=(20, 4.9), constrained_layout=True)
    Ls, e = a["Ls"], a["err"]
    for L, c in zip(Ls, L_COLORS):
        axes[0].plot((ps - P_C) * L ** (1 / NU), res[L]["cross"].mean(1), "-o", color=c, ms=3, label=f"L = {L}")
    axes[0].set(xlim=(-4, 4), xlabel=r"rescaled density $(p - p_c)\,L^{1/\nu}$,  $\nu = 4/3$",
                ylabel="P(fire crosses)", title="A  Data collapse: one universal curve")
    axes[0].legend(loc="upper left")

    x, c = np.array([24, 340]), L_COLORS[2]
    axes[1].errorbar(Ls, a["width"], e["width"], fmt="o", color=c, ms=8, label="measured transition width")
    axes[1].plot(x, a["width"][0] * (x / Ls[0]) ** (-1 / a["nu"]), color=c, lw=1.5,
                 label=rf"fit $\nu$ = {a['nu']:.3f} $\pm$ {e['nu']:.3f}")
    axes[1].plot(x, a["width"][0] * (x / Ls[0]) ** (-1 / NU), "--", color=MUTED, lw=1, label=r"exact $\nu = 4/3$")
    log_ticks(axes[1], Ls, [0.005, 0.01, 0.02])
    axes[1].set(xlabel="system size L", ylabel="width $w$ of the P(p) transition",
                title=r"B  Width shrinks as $L^{-1/\nu}$")
    axes[1].legend(loc="lower left")

    xs = Ls ** (-1 / NU)
    xx = np.linspace(0, xs.max() * 1.08, 50)
    axes[2].axhline(P_C, color=MUTED, ls="--", lw=1, label="$p_c$ = 0.592746 (Newman & Ziff 2000)")
    axes[2].plot(xx, np.polyval(a["lin"], xx), color=c, lw=1.5)
    axes[2].errorbar(xs, a["p_half"], e["p_half"], fmt="o", color=c, ms=8, label="midpoint of P(p) for each L")
    axes[2].errorbar([0], [a["p_c"]], [e["p_c"]], fmt="D", color=INK, ms=8, capsize=4,
                     label=f"extrapolated: {a['p_c']:.4f} $\\pm$ {e['p_c']:.4f}")
    for L, xv, yv in zip(Ls, xs, a["p_half"]):
        axes[2].annotate(f"L={L}", (xv, yv), xytext=(6, -14), textcoords="offset points", fontsize=9, color=INK2)
    axes[2].set(xlabel=r"$L^{-1/\nu}$   ($L \to \infty$ at 0)", ylabel="p where P(cross) = 1/2",
                xlim=(-0.004, xs.max() * 1.2), ylim=(P_C - 0.0035, P_C + 0.0035),
                title=r"C  Extrapolate to $L \to \infty$")
    axes[2].legend(loc="upper right", fontsize=9)

    axes[3].errorbar(Ls, a["dur_c"], e["dur_c"], fmt="o", color=c, ms=8, label="mean fire duration at $p_c$")
    axes[3].plot(x, a["dur_c"][0] * (x / Ls[0]) ** a["d_min"], color=c, lw=1.5,
                 label=rf"fit $d_{{min}}$ = {a['d_min']:.3f} $\pm$ {e['d_min']:.3f}")
    axes[3].plot(x, a["dur_c"][0] * (x / Ls[0]) ** D_MIN, "--", color=MUTED, lw=1, label="literature 1.1306")
    log_ticks(axes[3], Ls, [30, 100, 300, 1000])
    axes[3].set(xlabel="system size L", ylabel="CA steps", title=r"D  Fire duration $\sim L^{d_{min}}$")
    axes[3].legend(loc="upper left")
    fig.savefig(path, dpi=150)
    plt.close(fig)


def plot_ds(path, window, rho, warmup, theta, L, pl, rho_ss, rho_se):
    fig, axes = plt.subplots(1, 3, figsize=(17, 5.2), constrained_layout=True,
                             gridspec_kw=dict(width_ratios=[1, 1.35, 1.35]))
    axes[0].imshow(window, interpolation="nearest")
    axes[0].set_axis_off()
    axes[0].set_title(f"A  {window.shape[0]} x {window.shape[1]} window: fire scars")

    t = np.arange(rho.size) * theta / L**2
    tw = warmup * theta / L**2
    axes[1].plot(t, rho, color=AGE_COLOURS[0], lw=0.6)
    axes[1].axvspan(0, tw, color="#383835", alpha=0.6, lw=0)
    axes[1].text(tw / 2, 0.03, "warm-up\n(discarded)", ha="center", fontsize=9, color=MUTED)
    axes[1].axhline(P_C, color=MUTED, ls="--", lw=1)
    axes[1].text(t[-1], P_C + 0.01, "percolation threshold $p_c$ = 0.593", ha="right", fontsize=9, color=MUTED)
    axes[1].text(t[-1], 0.30, rf"steady state $\langle\rho\rangle$ = {rho_ss:.4f} $\pm$ {rho_se:.4f}" "\n"
                 f"Grassberger (2002), p/f = 1000: {RHO_DS}", ha="right", va="top", fontsize=10, color=INK2)
    axes[1].set(xlim=(0, t[-1]), ylim=(0, 0.66), xlabel="time (growth attempts per site)",
                ylabel=r"tree density $\rho$", title="B  The forest settles far below $p_c$")

    s, n, fit = pl["s"], pl["n"], pl["fit"]
    axes[2].loglog(s, n, "o", color=L_COLORS[2], ms=6, label=f"{pl['fires']:,} fires (log-binned)")
    axes[2].loglog(s[fit], np.exp(pl["icpt"]) * s[fit] ** -pl["tau"], color=INK, lw=1.5, ls="--",
                   label=rf"fit $s^{{-\tau}}$ on 10-2000: $\tau$ = {pl['tau']:.3f}" "\n(literature 1.11-1.19)")
    axes[2].set(xlabel="fire size s (trees burned)", ylabel="n(s), probability density",
                title="C  Fire sizes follow a power law")
    axes[2].legend(loc="lower left")
    fig.suptitle(f"Drossel-Schwabl forest-fire model, L = {L}, p/f = {theta}: self-organised criticality",
                 fontsize=14, weight="bold")
    fig.savefig(path, dpi=150)
    plt.close(fig)


# ----------------------------------------------------------------------------------------------------
def verify(a, rho_ss, rho_se, tau):
    """Every number the README claims is checked here against a known result."""
    rng = np.random.default_rng(0)
    for p in (0.5, P_C, 0.7):  # fast front-based CA == textbook dense CA, cell for cell
        tr = rng.random((64, 64)) < p
        assert (arrival_times(tr) == arrival_times_dense(tr)).all()
    assert arrival_times(np.ones((1, 8, 8), bool)).max() == 7  # full forest: fire crosses 8 cells in 7 steps
    print("  [PASS] front-based CA reproduces the dense CA exactly")
    e = a["err"]
    checks = [  # name, measured, +/- (bootstrap or batch means), expected, tolerance
        ("p_c   percolation threshold", a["p_c"], e["p_c"], P_C, 0.02),
        ("nu    correlation-length exponent", a["nu"], e["nu"], NU, 0.10),
        ("d_min shortest-path dimension", a["d_min"], e["d_min"], D_MIN, 0.05),
        ("rho   DS mean density, p/f=1000", rho_ss, rho_se, RHO_DS, 0.005),
    ]
    for name, got, se, want, tol in checks:
        ok = abs(got - want) < tol
        print(f"  [{'PASS' if ok else 'FAIL'}] {name}: {got:.4f} +/- {se:.4f}  expected {want:.4f}  (tol {tol})")
        assert ok, name
    lo, hi = TAU_FIRE
    print(f"  [{'PASS' if lo <= tau <= hi else 'NOTE'}] tau   DS fire-size exponent: {tau:.3f}  "
          f"literature {lo}-{hi} (an effective exponent, see README)")
    print("  P(cross) at p_c: " + ", ".join(f"L={L}: {v:.3f}" for L, v in a["P_at_pc"].items())
          + "   (Cardy 1992: -> 0.5 for a square)")


def main():
    OUT.mkdir(exist_ok=True)
    t0 = time.perf_counter()
    print("Part A: percolation sweep")
    ps, res = sweep()
    a = analyse(ps, res)
    plot_sweep(OUT / "percolation_sweep.png", ps, res)
    plot_fss(OUT / "finite_size_scaling.png", ps, res, a)
    print(f"  p_half(L) = {np.round(a['p_half'], 4)}, width(L) = {np.round(a['width'], 4)}, "
          f"duration at p_c = {np.round(a['dur_c'], 1)}")
    seed, T, frac = hero(OUT / "hero.png")
    print(f"  hero: seed {seed}, fire lasted {T} steps, burned {frac:.1%} of trees")
    gif_fire(OUT / "fire_spread.gif")
    print(f"  Part A done at {time.perf_counter() - t0:.0f} s")

    print("Part B: Drossel-Schwabl forest fire")
    L, theta, warm = 1024, 1000, 12_000
    t1 = time.perf_counter()
    sizes, rho, final = drossel_schwabl(L, theta, strikes=warm + 80_000, warmup=warm)
    pl = power_law(sizes)
    rho_ss = rho[warm:].mean()
    rho_se = rho[warm:].reshape(20, -1).mean(1).std(ddof=1) / np.sqrt(20)  # batch means: rho is autocorrelated
    print(f"  L={L}, p/f={theta}: {pl['fires']} fires, largest {sizes.max()}, <rho> = {rho_ss:.4f}, "
          f"tau = {pl['tau']:.3f}  ({time.perf_counter() - t1:.0f} s)")
    plot_ds(OUT / "forest_fire_soc.png", final[-1][:256, :256], rho, warm, theta, L, pl, rho_ss, rho_se)
    Lm, thm, every = 200, 300, 3
    wm = 20 * Lm * Lm // thm
    _, rho_m, frames = drossel_schwabl(Lm, thm, strikes=wm + 150 * every, warmup=wm, frame_every=every, seed=11)
    gif_ds(OUT / "forest_fire_soc.gif", frames[:-1], rho_m, wm, every, Lm, thm)
    print(f"  Part B done at {time.perf_counter() - t0:.0f} s")

    print("Verification")
    verify(a, rho_ss, rho_se, pl["tau"])
    print(f"Total runtime {time.perf_counter() - t0:.0f} s wall, {time.process_time():.0f} s CPU")


if __name__ == "__main__":
    main()