#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ Why small fires are common and big fires are not rare enough ============================================================ Science Journaling Club, Volume 1, Issue 2 (Winter 2025), theme "How Things Spread". THE QUESTION ------------ The Drossel-Schwabl forest fire model is two rules on a grid: trees appear at random, lightning falls at random, and whatever connected patch of trees the lightning lands in burns to the ground. It produces fires of every size from one cell to tens of thousands, and it is usually quoted as a textbook case of self-organised criticality, meaning that the distribution of fire sizes is said to be a power law with no characteristic scale. We ask three things. 1. Is the fire size distribution really a power law, tested rather than asserted, using maximum likelihood and a goodness-of-fit test instead of a straight line drawn through a log-log plot? 2. What is the exponent, and how does it depend on the ratio theta = p/f of the tree growth rate to the lightning rate? 3. Are lattices of the size a school club can run (512 x 512 up to 1024 x 1024) large enough to say anything about the long-running argument over whether this model is critical at all? THE MODEL, EXACTLY AS IMPLEMENTED --------------------------------- Square lattice of L x L cells, each empty or holding a tree, with helical boundary conditions: cells carry one index i, the neighbours of i are i-1, i+1, i-L and i+L, and i + L*L is identified with i. That is the same wiring Grassberger (2002) used, so our densities are directly comparable with his. Time is handled in the standard separation-of-time-scales form of the model. Each elementary event is either a growth attempt (probability theta/(1+theta)) or a lightning strike (probability 1/(1+theta)), so the number of growth attempts between successive strikes is geometric with mean theta = p/f. growth attempt : pick a cell uniformly at random. If it is empty, a tree appears. If it already holds a tree, nothing happens. lightning : pick a cell uniformly at random. If it is empty, nothing happens and no fire is recorded. If it holds a tree, the entire 4-connected cluster of trees containing that cell is removed in one instant, and its size s is recorded as one fire event. The fire is instantaneous relative to growth: no tree grows while a fire burns. That is the standard double separation of time scales, and it is what makes the model tractable and also what makes it an idealisation. ASSUMPTIONS THIS MODEL MAKES, WHICH REALITY DOES NOT ---------------------------------------------------- * Fire spreads to every touching tree with probability 1 and to nothing else. There is no wind, no slope, no spotting of embers across gaps, no humidity, no species that resist burning. * Trees appear independently at uniformly random cells. Real forests have seed dispersal, so new trees cluster near old ones. * A burnt cell becomes instantly available for regrowth. There is no soil damage, no seed bank, no succession. * Lightning is uniform in space and time. Real lightning is clustered in storms and correlated with drought. * The lattice is flat and homogeneous, with no roads, rivers or firebreaks. * Everything is dimensionless. A "cell" is not an acre and a "fire" is not a hectare. Nothing here is a prediction about any real forest anywhere. WHAT THE COMPUTATION IS ----------------------- This is a simulation, run on a laptop. No forest was observed, no fire was measured, and no field data of any kind enters this study. The computation IS the experiment: the object under study is the model itself, and every number below is output printed by this script. A NOTE ON CONVENTIONS --------------------- We write the fire size distribution as P(s) proportional to s^(-tau). Much of the physics literature, Grassberger (2002) included, writes it as s^(1-tau), so their tau is our tau plus one. Their "tau = 2.15" is our "tau = 1.15". Every published value quoted in this file has been converted to our convention. LIMITATIONS OF THIS PARTICULAR STUDY ------------------------------------ * Our largest lattice is 1024 x 1024 and our largest theta is 5000. The published work that reopened the argument about this model used lattices up to 65536 x 65536 and theta up to 256000 (Grassberger 2002). We are three to four orders of magnitude short in area. We can measure the drift in the exponent; we cannot claim to have resolved the asymptotic behaviour. * At our largest theta the biggest fires eat a large fraction of the lattice, so those runs are limited by lattice size and not by theta. We measure how much that matters by repeating one of them on a lattice four times larger. * Fires in a single run are correlated in time, so an ordinary bootstrap misstates the uncertainty on the exponent. We use a moving-block bootstrap as well and report both. * Our goodness-of-fit test uses the Clauset-Shalizi-Newman semi-parametric procedure with the fitting window held fixed across synthetic replicates, rather than reselected each time. That makes the test slightly generous. * The runs are not all the same length. The main run records 120,000 fires; theta = 2000 records 20,000 and theta = 5000 records 10,000, because a fire at large theta costs roughly its own size in work. The exponent at the top of the theta sweep therefore carries about twice the error of the main run, and we print that error rather than hiding the imbalance. * The fitting window is fixed at s in [4, theta] and chosen before any fit is made. Section 3 prints what happens across thirteen other windows, and the spread there is far larger than any statistical error in this file. That spread, not the bootstrap, is the real uncertainty on the exponent. HOW TO RUN ---------- python forest-fire-criticality.py > forest-fire-criticality-output.txt Needs numpy. Runtime 7 to 12 minutes on a laptop, depending on load. Every random number comes from numpy's PCG64 generator seeded with SEED below; rerunning reproduces every digit printed here. """ import os # All the linear algebra in this file is least squares on two columns, so BLAS # threading buys nothing and its per-thread buffers are the largest allocation # numpy makes at import time. Pinning it to one thread lets the script start on # a laptop that is short of memory. It does not change any number printed here. os.environ.setdefault("OPENBLAS_NUM_THREADS", "1") os.environ.setdefault("OMP_NUM_THREADS", "1") import gc import math import platform import sys import time import numpy as np SEED = 20251124 # Target number of random cell indices held in memory per batch of strikes. Kept # small on purpose: this laptop is often short of memory and a big temporary # array is the first thing to fail. PLANT_CHUNK = 50000 T_START = time.perf_counter() def banner(text): print() print("=" * 78) print(text) print("=" * 78) def sub(text): print() print("-- " + text + " " + "-" * max(0, 74 - len(text))) def elapsed(): return time.perf_counter() - T_START # --------------------------------------------------------------------------- # 1. The forest fire model # --------------------------------------------------------------------------- def draw_cells(rng, N, size): """rng.integers with a retry, because the machine is sometimes out of memory. numpy allocates the output array before it draws anything, so a call that fails with MemoryError leaves the generator state untouched. Repeating the identical call therefore produces exactly the same numbers it would have produced first time, and the run stays reproducible. """ last = None for attempt in range(6): try: return rng.integers(0, N, size=size, dtype=np.int32) except MemoryError as exc: last = exc gc.collect() time.sleep(2.0 * (attempt + 1)) raise last def warmup_for(L, theta, turnovers=12): """Fires to discard: enough to burn and regrow the whole lattice many times. One lattice-worth of trees is planted every L*L / (theta*(1-rho)) fires, and rho sits near 0.4, so we use 0.6 for (1-rho) and multiply by the number of turnovers we want to throw away. """ return max(3000, int(turnovers * (L * L) / (theta * 0.6))) def simulate(L, theta, n_fires, warmup_fires, rng, trace_every=0): """Run the Drossel-Schwabl model and return the sizes of n_fires fires. The lattice is a bytearray, which Python indexes quickly, and a numpy view over the same memory does the planting in one vectorised assignment per lightning interval. Helical wrapping costs one comparison per neighbour. """ N = L * L lat = bytearray(N) arr = np.frombuffer(lat, dtype=np.uint8) # shares memory with lat sizes = np.empty(n_fires, dtype=np.int64) trace = [] total_target = warmup_fires + n_fires seen = 0 # fires so far, warm-up included got = 0 # fires stored strikes = 0 # lightning strikes after warm-up warm_strikes = 0 p_geo = 1.0 / (theta + 1.0) # Lightning strikes are handled in batches so the planting between them can # be drawn with one call. A batch of B strikes needs about theta*B random # cell indices at once, so B has to shrink as theta grows or the temporary # array becomes hundreds of megabytes. We hold it near a quarter of a # million indices, and draw them as int32, which is ample for L <= 46340. BATCH = int(min(256, max(8, PLANT_CHUNK // max(theta, 1)))) while seen < total_target: gaps = rng.geometric(p_geo, size=BATCH) - 1 n_plant = int(gaps.sum()) plants = draw_cells(rng, N, n_plant) bolts = draw_cells(rng, N, BATCH) gap_list = gaps.tolist() bolt_list = bolts.tolist() off = 0 for b in range(BATCH): k = gap_list[b] if k: arr[plants[off:off + k]] = 1 off += k if seen < warmup_fires: warm_strikes += 1 else: strikes += 1 seed_cell = bolt_list[b] if lat[seed_cell]: lat[seed_cell] = 0 stack = [seed_cell] push = stack.append pop = stack.pop n = 1 while stack: c = pop() d = c - 1 if d < 0: d += N if lat[d]: lat[d] = 0 push(d) n += 1 d = c + 1 if d >= N: d -= N if lat[d]: lat[d] = 0 push(d) n += 1 d = c - L if d < 0: d += N if lat[d]: lat[d] = 0 push(d) n += 1 d = c + L if d >= N: d -= N if lat[d]: lat[d] = 0 push(d) n += 1 seen += 1 if seen > warmup_fires and got < n_fires: sizes[got] = n got += 1 if trace_every and seen % trace_every == 0: trace.append((seen, float(arr.sum()) / N)) if seen >= total_target: break return { "L": L, "theta": theta, "sizes": sizes[:got], "strikes": strikes, "warm_strikes": warm_strikes, "fires": got, "rho": got / strikes if strikes else float("nan"), "trace": trace, } # --------------------------------------------------------------------------- # 2. Fitting machinery: discrete power law on a finite support # --------------------------------------------------------------------------- # # The exponents at stake here are near 1.15, well below 2, so the unbounded # discrete power law (the zeta distribution) is no use: sum x^-tau diverges for # tau <= 1, and in any case the data have a real upper cutoff. We therefore fit # the power law truncated to the integer support [xmin, xmax]: # # P(x) = x^-tau / Z(tau), Z(tau) = sum_{k=xmin}^{xmax} k^-tau # # The log-likelihood is l(tau) = -tau * sum(log x) - n * log Z(tau), and the # maximum satisfies mean(log x)_data = mean(log x)_model(tau), with the right # hand side monotonically decreasing in tau. One bisection finds it. def support(xmin, xmax): xs = np.arange(xmin, xmax + 1, dtype=np.float64) return xs, np.log(xs) def model_moments(tau, logxs): w = np.exp(-tau * logxs) Z = w.sum() m1 = float((w * logxs).sum() / Z) m2 = float((w * logxs * logxs).sum() / Z) return m1, m2, Z def mle_tau(mean_log_data, logxs, lo=0.001, hi=8.0, iters=70): """Maximum likelihood exponent for the truncated discrete power law.""" for _ in range(iters): mid = 0.5 * (lo + hi) m1, _, _ = model_moments(mid, logxs) if m1 > mean_log_data: lo = mid else: hi = mid return 0.5 * (lo + hi) def ks_distance(counts, tau, logxs): w = np.exp(-tau * logxs) cdf_m = np.cumsum(w) cdf_m /= cdf_m[-1] cdf_e = np.cumsum(counts) / counts.sum() return float(np.abs(cdf_e - cdf_m).max()) def fit_window(data, xmin, xmax): """Fit sizes in [xmin, xmax]. Returns a dict of fit results.""" sel = data[(data >= xmin) & (data <= xmax)] n = sel.size if n < 20: return None xs, logxs = support(xmin, xmax) mean_log = float(np.log(sel).mean()) tau = mle_tau(mean_log, logxs) m1, m2, Z = model_moments(tau, logxs) var_log = max(m2 - m1 * m1, 1e-12) se = 1.0 / math.sqrt(n * var_log) # asymptotic, from the Fisher info counts = np.bincount(sel - xmin, minlength=xs.size).astype(np.float64) ks = ks_distance(counts, tau, logxs) return {"tau": tau, "se": se, "n": n, "xmin": xmin, "xmax": xmax, "ks": ks, "mean_log": mean_log, "var_log": var_log} def sample_truncated(tau, xs, logxs, n, rng): w = np.exp(-tau * logxs) cdf = np.cumsum(w) cdf /= cdf[-1] return xs[np.searchsorted(cdf, rng.random(n))].astype(np.int64) def gof_pvalue(data, xmin, xmax, n_synth, rng): """Clauset-Shalizi-Newman goodness of fit, window held fixed. Simulate n_synth data sets of the same size from the fitted law, refit each and compare its KS distance with the observed one. p is the fraction of synthetic sets that fit their own model no better than the data fits its model. A small p means the power law is rejected. """ fit = fit_window(data, xmin, xmax) if fit is None: return None xs, logxs = support(xmin, xmax) worse = 0 ks_synth = np.empty(n_synth) for i in range(n_synth): syn = sample_truncated(fit["tau"], xs, logxs, fit["n"], rng) t_s = mle_tau(float(np.log(syn).mean()), logxs) cnt = np.bincount(syn - xmin, minlength=xs.size).astype(np.float64) d = ks_distance(cnt, t_s, logxs) ks_synth[i] = d if d >= fit["ks"]: worse += 1 return {"p": worse / n_synth, "ks": fit["ks"], "tau": fit["tau"], "n": fit["n"], "ks_synth_mean": float(ks_synth.mean()), "ks_synth_95": float(np.quantile(ks_synth, 0.95))} def choose_block(n, target_blocks=60, cap=1000): """Block length for the moving-block bootstrap. We want enough blocks for the resampling distribution to mean something (at least a few dozen) and blocks long enough to carry the correlation between consecutive fires. This keeps roughly `target_blocks` of them. """ return int(max(50, min(cap, n // target_blocks))) def block_bootstrap_se(data, xmin, xmax, block, n_boot, rng): """Standard error on tau from a moving-block bootstrap of the fire sequence. Consecutive fires are not independent: a big fire clears the lattice and the next fires are small. Resampling contiguous blocks keeps that structure. Blocks are summarised once as (count, sum of logs) inside the window, so each bootstrap fit is a single bisection. """ inwin = (data >= xmin) & (data <= xmax) logs = np.where(inwin, np.log(np.maximum(data, 1)), 0.0) nb = data.size // block cnt_blocks = np.array([inwin[i * block:(i + 1) * block].sum() for i in range(nb)], dtype=np.float64) sum_blocks = np.array([logs[i * block:(i + 1) * block].sum() for i in range(nb)], dtype=np.float64) xs, logxs = support(xmin, xmax) taus = np.empty(n_boot) for b in range(n_boot): pick = rng.integers(0, nb, size=nb) n = cnt_blocks[pick].sum() if n < 20: taus[b] = np.nan continue taus[b] = mle_tau(sum_blocks[pick].sum() / n, logxs) taus = taus[np.isfinite(taus)] return (float(taus.std(ddof=1)), float(np.quantile(taus, 0.025)), float(np.quantile(taus, 0.975))) def loglog_regression(data, xmin, xmax, n_bins=0): """The naive method: histogram, take logs, fit a straight line by least squares. With n_bins = 0 the histogram is the raw integer histogram, which is what most people plot. Empty bins are dropped, which is exactly the step that breaks the estimator. With n_bins > 0 the counts go into logarithmically spaced bins and are divided by bin width first, which helps but does not cure it. """ sel = data[(data >= xmin) & (data <= xmax)] if sel.size < 20: return None if n_bins == 0: counts = np.bincount(sel, minlength=xmax + 1)[xmin:xmax + 1].astype(float) xs = np.arange(xmin, xmax + 1, dtype=float) keep = counts > 0 x = np.log10(xs[keep]) y = np.log10(counts[keep] / counts.sum()) else: edges = np.unique(np.round(np.logspace(np.log10(xmin), np.log10(xmax + 1), n_bins + 1)).astype(np.int64)) counts, _ = np.histogram(sel, bins=edges) width = np.diff(edges).astype(float) centre = np.sqrt(edges[:-1].astype(float) * edges[1:].astype(float)) keep = counts > 0 x = np.log10(centre[keep]) y = np.log10(counts[keep] / (width[keep] * sel.size)) A = np.vstack([x, np.ones_like(x)]).T coef, *_ = np.linalg.lstsq(A, y, rcond=None) resid = y - A @ coef s2 = float((resid ** 2).sum() / max(x.size - 2, 1)) cov = s2 * np.linalg.inv(A.T @ A) return {"slope": float(coef[0]), "tau": float(-coef[0]), "se": float(math.sqrt(cov[0, 0])), "points": int(x.size)} def log_binned_pdf(data, n_per_decade=8): """Logarithmically binned estimate of P(s), for plotting.""" smax = int(data.max()) edges = np.unique(np.round(np.logspace( 0, math.log10(smax + 1), int(math.ceil(math.log10(smax + 1) * n_per_decade)) + 1)).astype(np.int64)) counts, _ = np.histogram(data, bins=edges) width = np.diff(edges).astype(float) centre = np.sqrt(edges[:-1].astype(float) * edges[1:].astype(float)) pdf = counts / (width * data.size) keep = counts > 0 return centre[keep], pdf[keep], counts[keep] def powerfit(x, y): """Least squares slope of log10 y on log10 x, with its standard error.""" A = np.vstack([np.log10(x), np.ones_like(x, dtype=float)]).T coef, *_ = np.linalg.lstsq(A, np.log10(y), rcond=None) resid = np.log10(y) - A @ coef s2 = float((resid ** 2).sum() / max(len(x) - 2, 1)) cov = s2 * np.linalg.inv(A.T @ A) return float(coef[0]), float(math.sqrt(cov[0, 0])) # --------------------------------------------------------------------------- # 3. Run # --------------------------------------------------------------------------- def main(): master = np.random.default_rng(SEED) rng_sim, rng_fit, rng_val = master.spawn(3) banner("FOREST FIRE CRITICALITY - Science Journaling Club") print("Drossel-Schwabl forest fire model on a square lattice, helical boundary.") print("Every number below is printed by this script. No field data is used.") print("Convention: P(s) ~ s^(-tau). Physics papers usually write s^(1-tau),") print("so their exponent is ours plus one. Published values quoted here are") print("already converted to our convention.") print() print(" master seed : %d" % SEED) print(" generator : numpy PCG64 via default_rng, 3 spawned streams") print(" python : %s" % sys.version.split()[0]) print(" numpy : %s" % np.__version__) print(" platform : %s" % platform.platform()) # ---------------------------------------------------------------- part 1 banner("PART 1. VALIDATING THE FITTING MACHINERY ON DATA WE MADE OURSELVES") print("Before fitting anything from the forest, we check that the estimator can") print("recover an exponent we already know. Synthetic draws from a truncated") print("discrete power law, fitted by the same code that will fit the fires.") print() print("%-10s %-9s %-14s %10s %10s %12s %8s" % ("tau_true", "n", "support", "tau_hat", "se", "difference", "z")) val_rows = [] for tau_true, n_syn, xmin_s, xmax_s in [ (1.16, 200000, 1, 10000), (1.16, 20000, 1, 10000), (1.50, 200000, 1, 10000), (2.50, 200000, 1, 10000), (1.16, 200000, 10, 3000)]: xs, logxs = support(xmin_s, xmax_s) syn = sample_truncated(tau_true, xs, logxs, n_syn, rng_val) fit = fit_window(syn, xmin_s, xmax_s) diff = fit["tau"] - tau_true z = diff / fit["se"] print("%-10.4f %-9d [%5d,%6d] %10.5f %10.5f %+12.5f %8.2f" % (tau_true, n_syn, xmin_s, xmax_s, fit["tau"], fit["se"], diff, z)) val_rows.append((tau_true, n_syn, xmin_s, xmax_s, fit["tau"], fit["se"], diff, z)) print() print("A |z| below about 2 means the estimator recovered the truth inside its") print("own stated error. That is the whole test.") sub("the same synthetic data, fitted the popular way instead") print("A straight line through a log-log histogram is what almost every") print("undergraduate lab does. On data with a KNOWN answer, here is what it gives.") print() print("%-10s %-12s %10s %10s %12s %10s %10s" % ("tau_true", "method", "tau_hat", "se", "error", "err/true", "points")) bias_rows = [] for tau_true in (1.16, 1.50, 2.50): xs, logxs = support(1, 10000) syn = sample_truncated(tau_true, xs, logxs, 200000, rng_val) f_mle = fit_window(syn, 1, 10000) r_raw = loglog_regression(syn, 1, 10000, n_bins=0) r_bin = loglog_regression(syn, 1, 10000, n_bins=20) for name, r in (("MLE", {"tau": f_mle["tau"], "se": f_mle["se"], "points": f_mle["n"]}), ("OLS raw", r_raw), ("OLS logbin", r_bin)): err = r["tau"] - tau_true print("%-10.2f %-12s %10.5f %10.5f %+12.5f %9.1f%% %10d" % (tau_true, name, r["tau"], r["se"], err, 100 * err / tau_true, r["points"])) bias_rows.append((tau_true, name, r["tau"], r["se"], err, 100 * err / tau_true)) print() # ---------------------------------------------------------------- part 2 banner("PART 2. THE MAIN RUN") L_MAIN = 512 THETA_MAIN = 500 N_MAIN = 120000 warm = warmup_for(L_MAIN, THETA_MAIN) print("lattice : %d x %d = %d cells, helical boundary" % (L_MAIN, L_MAIN, L_MAIN * L_MAIN)) print("theta = p/f : %d growth attempts per lightning strike, on average" % THETA_MAIN) print("warm-up fires : %d (discarded)" % warm) print("recorded fires : %d" % N_MAIN) t0 = time.perf_counter() main_run = simulate(L_MAIN, THETA_MAIN, N_MAIN, warm, rng_sim, trace_every=500) print("wall time : %.1f s" % (time.perf_counter() - t0)) sizes = main_run["sizes"] rho_hit = main_run["rho"] sub("steady state and the mass balance identity") trace = main_run["trace"] tr_idx = np.array([t[0] for t in trace], dtype=float) tr_rho = np.array([t[1] for t in trace], dtype=float) warm_mask = tr_idx <= warm rec_mask = ~warm_mask print("density of trees, sampled every 500 fires (just after each fire):") print(" first 10 samples of the warm-up %s" % " ".join("%.3f" % v for v in tr_rho[:10])) print(" during warm-up (%d samples) mean %.5f sd %.5f" % (warm_mask.sum(), tr_rho[warm_mask].mean(), tr_rho[warm_mask].std())) print(" first half of recording mean %.5f" % tr_rho[rec_mask][:rec_mask.sum() // 2].mean()) print(" second half of recording mean %.5f" % tr_rho[rec_mask][rec_mask.sum() // 2:].mean()) rho_bar = float(tr_rho[rec_mask].mean()) print(" recording, all of it mean %.5f sd %.5f" % (rho_bar, tr_rho[rec_mask].std())) print(" fraction of strikes that hit a tree %.5f" % rho_hit) print(" (that fraction IS the density a lightning bolt actually sees, and it") print(" is the number we compare with the published tables)") print() print("In a steady state every tree that grows must eventually burn. Per") print("lightning strike, theta*(1-rho) trees appear and rho* trees burn, so") print("the mean fire size is fixed by the density alone:") print() s_pred = THETA_MAIN * (1.0 - rho_hit) / rho_hit s_meas = float(sizes.mean()) s_se = float(sizes.std(ddof=1) / math.sqrt(sizes.size)) print(" analytic = theta (1 - rho) / rho = %10.3f" % s_pred) print(" club value measured = %10.3f +/- %.3f" % (s_meas, s_se)) print(" difference = %+10.3f (%.3f%% of the analytic value)" % (s_meas - s_pred, 100 * (s_meas - s_pred) / s_pred)) print(" difference in standard errors = %+10.2f" % ((s_meas - s_pred) / s_se)) print() print("This is an exact identity for the model, so it is a check on the simulator") print("rather than a discovery. If the lattice were not in a steady state, or the") print("burning code missed cells, the two numbers would part company.") sub("what the fires look like") print(" number of fires recorded : %d" % sizes.size) print(" lightning strikes : %d" % main_run["strikes"]) print(" total cells burnt : %d" % int(sizes.sum())) print(" smallest fire : %d" % sizes.min()) print(" largest fire : %d cells (%.3f%% of the lattice)" % (sizes.max(), 100.0 * sizes.max() / (L_MAIN * L_MAIN))) print(" median fire : %d" % int(np.median(sizes))) print(" mean fire : %.2f" % s_meas) print(" fires of size 1 : %d (%.2f%% of all fires)" % ((sizes == 1).sum(), 100.0 * (sizes == 1).sum() / sizes.size)) print(" fires of size >= 1000 : %d (%.3f%% of fires, %.1f%% of burnt area)" % ((sizes >= 1000).sum(), 100.0 * (sizes >= 1000).sum() / sizes.size, 100.0 * sizes[sizes >= 1000].sum() / sizes.sum())) m1 = float(sizes.mean()) m2 = float((sizes.astype(float) ** 2).mean()) print(" / (cutoff proxy) : %.1f" % (m2 / m1)) for q in (0.5, 0.9, 0.99, 0.999, 0.9999): print(" quantile %-8.4f : %d" % (q, int(np.quantile(sizes, q)))) # ---------------------------------------------------------------- part 3 banner("PART 3. IS IT A POWER LAW?") print("Three fits of the same 120,000 fires. The first covers everything from") print("single-cell fires to the largest fire seen. The second covers the window") print("the theory actually predicts a power law for, s from 4 up to theta, which") print("stops below the cutoff. The third lets an automatic rule pick the lower") print("edge, which turns out to be a bad idea here and we say why.") full = fit_window(sizes, 1, int(sizes.max())) sub("fit A: the whole range, s = 1 to %d" % int(sizes.max())) print(" tau (MLE) : %.4f +/- %.4f (asymptotic)" % (full["tau"], full["se"])) print(" events used : %d" % full["n"]) print(" KS distance : %.5f" % full["ks"]) t0 = time.perf_counter() g_full = gof_pvalue(sizes, 1, int(sizes.max()), 100, rng_fit) print(" goodness of fit p : %.3f (100 synthetic sets, %.1f s)" % (g_full["p"], time.perf_counter() - t0)) print(" synthetic KS, mean / 95%% : %.5f / %.5f" % (g_full["ks_synth_mean"], g_full["ks_synth_95"])) print(" The observed KS is %.0f times the 95th percentile of what the fitted" % (g_full["ks"] / g_full["ks_synth_95"])) print(" power law itself produces. Rejected, and not narrowly.") XMIN_MAIN, XMAX_MAIN = 4, THETA_MAIN best = fit_window(sizes, XMIN_MAIN, XMAX_MAIN) sub("fit B (our primary estimate): s = %d to %d" % (XMIN_MAIN, XMAX_MAIN)) print("The upper edge is tied to theta because the cutoff scales with theta, so") print("this window keeps the same position relative to the cutoff at every theta.") print() print(" tau (MLE) : %.4f" % best["tau"]) print(" asymptotic standard error : %.4f" % best["se"]) t0 = time.perf_counter() BLK_MAIN = choose_block(sizes.size) bs_se, bs_lo, bs_hi = block_bootstrap_se(sizes, XMIN_MAIN, XMAX_MAIN, block=BLK_MAIN, n_boot=400, rng=rng_fit) print(" block bootstrap SE : %.4f (400 resamples, blocks of %d fires, %.1f s)" % (bs_se, BLK_MAIN, time.perf_counter() - t0)) print(" block bootstrap 95%% CI : [%.4f, %.4f]" % (bs_lo, bs_hi)) print(" ratio to asymptotic SE : %.2f" % (bs_se / best["se"])) print(" events used : %d" % best["n"]) print(" KS distance : %.5f" % best["ks"]) t0 = time.perf_counter() g_win = gof_pvalue(sizes, XMIN_MAIN, XMAX_MAIN, 200, rng_fit) print(" goodness of fit p : %.3f (200 synthetic sets, %.1f s)" % (g_win["p"], time.perf_counter() - t0)) print(" synthetic KS, mean / 95%% : %.5f / %.5f" % (g_win["ks_synth_mean"], g_win["ks_synth_95"])) TAU_MAIN = best["tau"] SE_MAIN = bs_se sub("fit C: letting an automatic rule choose the lower edge") print("The Clauset-Shalizi-Newman recipe picks the s_min that minimises the KS") print("distance. It was designed for distributions with no upper cutoff. With") print("the upper edge pinned at theta, shrinking the window always lowers KS,") print("so the rule walks off to the right and lands on a window too narrow to") print("mean anything. Here is the whole scan, so you can see it happen.") print() print("%8s %10s %10s %10s %10s %10s" % ("s_min", "n in win", "decades", "tau", "se", "KS")) best_ks = None xmin_rows = [] for xmin in (1, 2, 3, 4, 6, 8, 12, 16, 24, 32, 48, 64, 96, 128): f = fit_window(sizes, xmin, XMAX_MAIN) if f is None: continue dec = math.log10(XMAX_MAIN / xmin) print("%8d %10d %10.2f %10.4f %10.4f %10.5f" % (xmin, f["n"], dec, f["tau"], f["se"], f["ks"])) xmin_rows.append((xmin, f["n"], dec, f["tau"], f["se"], f["ks"])) if best_ks is None or f["ks"] < best_ks["ks"]: best_ks = f print() print(" KS is smallest at s_min = %d, which leaves only %.2f decades of window" % (best_ks["xmin"], math.log10(XMAX_MAIN / best_ks["xmin"]))) print(" and gives tau = %.4f, which is %+.4f from the fit B value of %.4f and" % (best_ks["tau"], best_ks["tau"] - TAU_MAIN, TAU_MAIN)) print(" %.1f times the fit B standard error away from it." % (abs(best_ks["tau"] - TAU_MAIN) / best["se"])) print(" We report fit B and treat this as a warning about automatic rules: the") print(" window the rule lands on moves from run to run, and the exponent moves") print(" with it.") sub("the naive log-log regression, on our own fire data") naive_rows = [] for name, kw in (("raw histogram", dict(n_bins=0)), ("log-binned, 20 bins", dict(n_bins=20))): r = loglog_regression(sizes, XMIN_MAIN, XMAX_MAIN, **kw) print(" %-22s tau = %.4f +/- %.4f (%d points) difference from MLE %+0.4f" % (name, r["tau"], r["se"], r["points"], r["tau"] - TAU_MAIN)) naive_rows.append((name, r["tau"], r["se"], r["points"])) sub("how much the answer depends on where we put the window") print("%10s %10s %10s %10s %10s %10s" % ("s_min", "s_max", "n", "tau", "se", "KS")) sens_rows = [] for xmin, xmax in [(1, 100), (1, 500), (1, 2000), (1, int(sizes.max())), (4, 250), (4, 500), (4, 1000), (4, 2000), (4, 5000), (16, 500), (16, 2000), (32, 1000), (64, 4000)]: f = fit_window(sizes, xmin, xmax) if f is None: continue print("%10d %10d %10d %10.4f %10.4f %10.5f" % (xmin, xmax, f["n"], f["tau"], f["se"], f["ks"])) sens_rows.append((xmin, xmax, f["n"], f["tau"], f["se"], f["ks"])) taus_sens = np.array([r[3] for r in sens_rows]) print() print(" spread of tau across these %d windows: %.4f to %.4f, range %.4f" % (len(sens_rows), taus_sens.min(), taus_sens.max(), np.ptp(taus_sens))) print(" The choice of window moves the answer by far more than the statistical") print(" error on any one of them. That is the real uncertainty in this number.") sub("the goodness of fit test has enormous power with this many fires") print("With enough data, a test rejects any model that is even slightly wrong.") print("Same window, same fitted exponent, smaller random subsets of the fires:") print() print("%12s %12s %10s %10s" % ("subsample", "tau", "KS", "p")) power_rows = [] win_data = sizes[(sizes >= XMIN_MAIN) & (sizes <= XMAX_MAIN)] for nsub in (300, 1000, 3000, 10000, 30000, win_data.size): if nsub > win_data.size: continue sub_data = rng_fit.choice(win_data, size=nsub, replace=False) g = gof_pvalue(sub_data, XMIN_MAIN, XMAX_MAIN, 100, rng_fit) print("%12d %12.4f %10.5f %10.3f" % (sub_data.size, g["tau"], g["ks"], g["p"])) power_rows.append((int(sub_data.size), g["tau"], g["ks"], g["p"])) # ---------------------------------------------------------------- part 4 banner("PART 4. CONVERGENCE OF THE ESTIMATE") print("Running maximum likelihood exponent over the fit B window, as fires") print("accumulate. The error column is one asymptotic standard error.") print() print("%12s %12s %10s %12s %12s" % ("fires used", "n in window", "tau", "se", "tau-tau_final")) conv_rows = [] for frac in (0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.3, 0.5, 0.7, 0.85, 1.0): k = max(200, int(frac * sizes.size)) f = fit_window(sizes[:k], XMIN_MAIN, XMAX_MAIN) print("%12d %12d %10.4f %12.4f %+12.4f" % (k, f["n"], f["tau"], f["se"], f["tau"] - TAU_MAIN)) conv_rows.append((k, f["n"], f["tau"], f["se"])) # ---------------------------------------------------------------- part 5 banner("PART 5. DEPENDENCE ON THETA = p/f") THETAS = [50, 125, 500, 2000, 5000] NFIRES = {50: 100000, 125: 100000, 500: N_MAIN, 2000: 20000, 5000: 10000} print("Five values spanning two decades, each with its own run at L = %d." % L_MAIN) print("theta = 125 is there because Grassberger (2002) tabulates the tree") print("density at exactly that value, which gives us a direct comparison.") print("The fitting window is s in [4, theta] throughout.") print() sweep = {} for th in THETAS: if th == THETA_MAIN: sweep[th] = main_run continue w = warmup_for(L_MAIN, th) t0 = time.perf_counter() sweep[th] = simulate(L_MAIN, th, NFIRES[th], w, rng_sim) print(" theta = %-6d %7d fires warm-up %6d %6.1f s" % (th, NFIRES[th], w, time.perf_counter() - t0)) print() print("%7s %8s %9s %10s %10s %8s %8s %9s %8s %8s %7s" % ("theta", "fires", "rho", " meas", " exact", "diff%", "max s", "/", "tau", "se(bb)", "gof p")) sweep_rows = [] for th in THETAS: run = sweep[th] s = run["sizes"] rho = run["rho"] pred = th * (1 - rho) / rho mm1 = float(s.mean()) mm2 = float((s.astype(float) ** 2).mean()) f = fit_window(s, 4, th) bse, blo, bhi = block_bootstrap_se(s, 4, th, block=choose_block(s.size), n_boot=200, rng=rng_fit) g = gof_pvalue(s, 4, th, 100, rng_fit) print("%7d %8d %9.5f %10.2f %10.2f %+8.3f %8d %9.1f %8.4f %8.4f %7.3f" % (th, s.size, rho, mm1, pred, 100 * (mm1 - pred) / pred, s.max(), mm2 / mm1, f["tau"], bse, g["p"])) sweep_rows.append({"theta": th, "n": int(s.size), "rho": rho, "mean": mm1, "pred": pred, "maxs": int(s.max()), "cut": mm2 / mm1, "tau": f["tau"], "se": bse, "ci": (blo, bhi), "p": g["p"], "q999": float(np.quantile(s, 0.999)), "ks": f["ks"], "nwin": f["n"], "frac": s.max() / float(L_MAIN * L_MAIN)}) print() print("The last column is the fraction of the lattice the biggest single fire") print("consumed, and it is the reason to distrust the largest theta:") for r in sweep_rows: print(" theta %5d : largest fire covered %5.1f%% of the lattice" % (r["theta"], 100 * r["frac"])) sub("how the cutoff moves with theta") th_arr = np.array([r["theta"] for r in sweep_rows], dtype=float) cut_arr = np.array([r["cut"] for r in sweep_rows], dtype=float) q999 = np.array([r["q999"] for r in sweep_rows], dtype=float) mx = np.array([r["maxs"] for r in sweep_rows], dtype=float) cut_slope, cut_slope_se = powerfit(th_arr, cut_arr) for label, arr in (("/", cut_arr), ("99.9th percentile", q999), ("largest fire seen", mx)): sl, sle = powerfit(th_arr, arr) print(" %-20s scales as theta^%.3f +/- %.3f" % (label, sl, sle)) print() print(" Grassberger (2002) reports the cutoff growing roughly as theta^1.08 in") print(" the large-theta limit, on lattices far bigger than ours.") sub("does the exponent drift with theta?") tau_arr = np.array([r["tau"] for r in sweep_rows]) se_arr = np.array([r["se"] for r in sweep_rows]) print(" " + ", ".join("theta %d: %.4f+/-%.4f" % (r["theta"], r["tau"], r["se"]) for r in sweep_rows)) wts = 1.0 / se_arr ** 2 tau_w = float((tau_arr * wts).sum() / wts.sum()) chi2 = float((wts * (tau_arr - tau_w) ** 2).sum()) dof = len(tau_arr) - 1 print(" weighted mean tau : %.4f" % tau_w) print(" chi-square about that mean : %.1f on %d degrees of freedom" % (chi2, dof)) print(" a constant exponent is %s by these five runs." % ("consistent with" if chi2 < 2 * dof else "NOT supported")) A = np.vstack([np.log10(th_arr), np.ones_like(th_arr)]).T drift_coef, *_ = np.linalg.lstsq(A, tau_arr, rcond=None) drift = float(drift_coef[0]) print(" drift: d(tau)/d(log10 theta) : %+.4f per decade" % drift) print(" extrapolating that drift, tau would reach 1.19 at theta = 10^%.1f," % ((1.19 - drift_coef[1]) / drift)) print(" which is a straight-line guess well outside the range we measured and") print(" should not be believed. We print it to show how far away the answer is.") # ---------------------------------------------------------------- part 6 banner("PART 6. DOES THE LATTICE SIZE MATTER?") print("Same theta = %d, three lattice sizes. If the model were simply critical" % THETA_MAIN) print("and our lattice big enough, the exponent would not care.") print() print("%8s %10s %9s %9s %10s %9s %9s %9s" % ("L", "cells", "fires", "rho", "", "max s", "tau", "se(bb)")) size_rows = [] for L in (256, 512, 1024): if L == L_MAIN: run = main_run else: run = simulate(L, THETA_MAIN, 30000, warmup_for(L, THETA_MAIN), rng_sim) s = run["sizes"] f = fit_window(s, 4, THETA_MAIN) bse, _, _ = block_bootstrap_se(s, 4, THETA_MAIN, block=choose_block(s.size), n_boot=200, rng=rng_fit) print("%8d %10d %9d %9.5f %10.2f %9d %9.4f %9.4f" % (L, L * L, s.size, run["rho"], s.mean(), s.max(), f["tau"], bse)) size_rows.append({"L": L, "n": int(s.size), "rho": run["rho"], "mean": float(s.mean()), "tau": f["tau"], "se": bse, "maxs": int(s.max())}) sub("the same check at theta = 5000, where the lattice is clearly too small") t0 = time.perf_counter() big = simulate(1024, 5000, 6000, warmup_for(1024, 5000), rng_sim) sb = big["sizes"] fb = fit_window(sb, 4, 5000) bseb, _, _ = block_bootstrap_se(sb, 4, 5000, block=choose_block(sb.size), n_boot=200, rng=rng_fit) small = sweep[5000] print(" %-34s %14s %14s" % ("", "L = 512", "L = 1024")) print(" %-34s %14d %14d" % ("fires", small["sizes"].size, sb.size)) print(" %-34s %14.5f %14.5f" % ("density rho", small["rho"], big["rho"])) print(" %-34s %14.2f %14.2f" % ("mean fire size", small["sizes"].mean(), sb.mean())) print(" %-34s %14d %14d" % ("largest fire", small["sizes"].max(), sb.max())) print(" %-34s %13.1f%% %13.1f%%" % ("largest fire, % of lattice", 100 * small["sizes"].max() / 512 ** 2, 100 * sb.max() / 1024 ** 2)) print(" %-34s %14.1f %14.1f" % ("/", float((small["sizes"].astype(float) ** 2).mean() / small["sizes"].mean()), float((sb.astype(float) ** 2).mean() / sb.mean()))) print(" %-34s %14.4f %14.4f" % ("tau, window [4, 5000]", sweep_rows[-1]["tau"], fb["tau"])) print(" %-34s %14.4f %14.4f" % ("block bootstrap SE", sweep_rows[-1]["se"], bseb)) print(" wall time for the L = 1024 run: %.1f s" % (time.perf_counter() - t0)) print() print(" difference in tau between the two lattice sizes: %+.4f, which is %.1f" % (fb["tau"] - sweep_rows[-1]["tau"], abs(fb["tau"] - sweep_rows[-1]["tau"]) / math.hypot(bseb, sweep_rows[-1]["se"]))) print(" combined standard errors.") # ---------------------------------------------------------------- part 7 banner("PART 7. THE CLUB'S NUMBER BESIDE THE PUBLISHED ONES") print("%-46s %11s %11s %11s" % ("quantity", "club", "published", "difference")) print("-" * 82) def cmp_row(name, ours, theirs, note, fmt="%11.4f"): print(("%-46s " + fmt + " " + fmt + " " + fmt) % (name, ours, theirs, ours - theirs)) if note: print(" %s" % note) tau5000 = sweep_rows[-1]["tau"] cmp_row("tau at theta = 500 (main run)", TAU_MAIN, 1.15, "value from the pre-2002 simulations, quoted by Grassberger 2002") cmp_row("tau at theta = 5000 (widest window)", tau5000, 1.15, "same comparison at our largest theta") cmp_row("tau at theta = 5000 (widest window)", tau5000, 1.19, "Grassberger 2002 own estimate, tau_G = 2.19 +/- 0.01") cmp_row("tau at theta = 5000 (widest window)", tau5000, 1.111, "Grassberger 2002 envelope exponent, tau_G' = 2.111 +/- 0.006") rho125 = [r for r in sweep_rows if r["theta"] == 125][0]["rho"] cmp_row("tree density rho at theta = 125", rho125, 0.379837, "Grassberger 2002 Table 1, rho = 0.379837 +/- 0.000006", "%11.5f") cmp_row("cutoff scaling exponent", cut_slope, 1.08, "Grassberger 2002, s_max ~ theta^1.08", "%11.3f") cmp_row("mean fire size vs exact identity (%)", 100 * (s_meas - s_pred) / s_pred, 0.0, "exact steady-state mass balance, no fitting involved", "%11.4f") print() n_rho = [r for r in sweep_rows if r["theta"] == 125][0]["n"] rho_se = math.sqrt(rho125 * (1 - rho125) / (n_rho / rho125)) print("The density comparison is the sharpest one we have. Our theta = 125 run") print("gives rho = %.5f from %d lightning strikes, a binomial standard error of" % (rho125, int(n_rho / rho125))) print("%.5f, so we sit %.1f standard errors from Grassberger's tabulated %.6f." % (rho_se, abs(rho125 - 0.379837) / rho_se, 0.379837)) print() print("Our exponent at theta = 500 sits %+.4f from 1.15, with a block-bootstrap" % (TAU_MAIN - 1.15)) print("standard error of %.4f. At theta = 5000 it has climbed to %.4f. The" % (SE_MAIN, tau5000)) print("exponent is not a constant of the model at the sizes we can reach: it") print("drifts upward with theta by %+.4f per decade, and the published argument" % drift) print("is about where, or whether, that drift stops.") # ---------------------------------------------------------------- figures banner("PART 8. NUMBERS FOR THE FIGURES") sub("fig 1: fire size distribution, main run, log-binned, 8 bins per decade") c, p, k = log_binned_pdf(sizes, 8) print("%14s %16s %12s" % ("s (centre)", "P(s)", "count")) for i in range(len(c)): print("%14.3f %16.8e %12d" % (c[i], p[i], k[i])) sub("fig 2: compensated distributions s^tau P(s), tau = %.4f, all five theta" % TAU_MAIN) for th in THETAS: cc, pp, kk = log_binned_pdf(sweep[th]["sizes"], 6) print() print("theta = %d" % th) print("%14s %16s %16s %10s" % ("s", "P(s)", "s^tau P(s)", "count")) for i in range(len(cc)): print("%14.3f %16.8e %16.8e %10d" % (cc[i], pp[i], (cc[i] ** TAU_MAIN) * pp[i], kk[i])) sub("fig 3: convergence of tau with fires accumulated") print("%12s %12s %12s" % ("fires", "tau", "se")) for k_, nw, t_, se_ in conv_rows: print("%12d %12.5f %12.5f" % (k_, t_, se_)) sub("fig 3b: tree density trace, main run, every 500 fires, thinned") step = max(1, len(tr_idx) // 60) print("%12s %12s" % ("fire index", "density")) for i in range(0, len(tr_idx), step): print("%12d %12.5f" % (int(tr_idx[i]), tr_rho[i])) sub("fig 4: estimator bias, MLE against log-log regression") print("%10s %-12s %12s %12s" % ("tau_true", "method", "tau_hat", "error%")) for row in bias_rows: print("%10.2f %-12s %12.5f %+11.1f%%" % (row[0], row[1], row[2], row[5])) print() print("and on the club's own fire data, where nobody knows the true answer:") print(" %-22s tau = %.4f" % ("MLE, window [4, 500]", TAU_MAIN)) for name, tval, tse, npts in naive_rows: print(" %-22s tau = %.4f" % (name, tval)) sub("fig 5: exponent and cutoff against theta") print("%10s %10s %10s %12s %12s %12s %10s" % ("theta", "tau", "se", "/", "q99.9", "max s", "gof p")) for r in sweep_rows: print("%10d %10.4f %10.4f %12.1f %12.1f %12d %10.3f" % (r["theta"], r["tau"], r["se"], r["cut"], r["q999"], r["maxs"], r["p"])) sub("fig 5b: exponent against lattice size at theta = %d" % THETA_MAIN) print("%10s %10s %10s %12s" % ("L", "tau", "se", "")) for r in size_rows: print("%10d %10.4f %10.4f %12.2f" % (r["L"], r["tau"], r["se"], r["mean"])) sub("fig 5c: the window sensitivity table again, compactly") print("%10s %10s %10s %10s" % ("s_min", "s_max", "tau", "se")) for xmin, xmax, n, t_, se_, ks_ in sens_rows: print("%10d %10d %10.4f %10.4f" % (xmin, xmax, t_, se_)) banner("SUMMARY") total_fires = (sum(sweep[t]["sizes"].size for t in THETAS) + sum(r["n"] for r in size_rows if r["L"] != L_MAIN) + sb.size) print(" seed : %d" % SEED) print(" main run : L = %d, theta = %d, %d fires" % (L_MAIN, THETA_MAIN, sizes.size)) print(" fires recorded across every run : %d" % total_fires) print(" tau, primary window s in [%d, %d] : %.4f +/- %.4f (block bootstrap)" % (XMIN_MAIN, XMAX_MAIN, TAU_MAIN, SE_MAIN)) print(" tau, whole range : %.4f (rejected, p = %.3f)" % (full["tau"], g_full["p"])) print(" goodness of fit p, primary window : %.3f" % g_win["p"]) print(" tau at theta = 5000 : %.4f +/- %.4f" % (tau5000, sweep_rows[-1]["se"])) print(" naive log-log regression, same window : %.4f" % naive_rows[0][1]) print(" exponent drift per decade of theta : %+.4f" % drift) print(" cutoff scaling : theta^%.3f +/- %.3f" % (cut_slope, cut_slope_se)) print(" density at theta = 125 vs Grassberger : %.5f vs 0.379837" % rho125) print(" total wall time : %.1f s" % elapsed()) print() print("Nothing in this file was measured in a forest. It is the output of a") print("program, and the program is the experiment.") if __name__ == "__main__": main()