#!/usr/bin/env python3 """ percolation-threshold.py Science Journaling Club, Volume 1 Issue 2 (Winter 2025), theme "How Things Spread". THE QUESTION ------------ Fill the cells of a square grid at random, each cell independently occupied with probability p. At small p you get scattered islands. At large p you get one blob with holes in it. Somewhere in between there is a value of p at which a connected path of occupied cells first appears that runs all the way from the top edge to the bottom edge, and in the limit of an infinite grid that value is sharp. It is called the site percolation threshold, p_c. For the square lattice it has no known closed form, so the only way to know it is to measure it. The accepted value from large specialised simulations is 0.59274605. We ask two things: (1) Can a student laptop, running ordinary CPython, measure that constant well enough for the measurement to be worth anything, and how big is the honest error bar? (2) At the threshold, how does the size of the largest connected cluster grow with the side length L of the grid? The prediction is a power law with exponent 91/48 = 1.895833..., the fractal dimension of the incipient infinite cluster. THE MODEL (this is a computation, not an observation) ----------------------------------------------------- There is no experiment here in the physical sense. Nothing was poured through anything. The lattice is an array in memory and "percolation" means that a graph search finds a path. Everything below is a statement about a random graph. The reason it is interesting is that random graphs of this kind are believed, and in some cases proved, to share their critical exponents with real disordered systems: conduction through a random resistor network, gelation, spread on a spatial contact network. The threshold value 0.59274605 is specific to this lattice and universal to nothing. The exponents are the universal part. Lattice L x L square grid, L in {16, 32, 64, 128, 256, 512, 1024}. Neighbours 4-connected (up, down, left, right). Diagonals do NOT connect. Boundaries free (open) in both directions. No wrap-around, no periodicity. Spanning at least one connected cluster of occupied sites contains a site in row 0 and a site in row L-1. Left-right crossing is not counted. Occupation each site independently occupied with probability p. ALGORITHM --------- Weighted union-find (disjoint-set forest) with union by size and path compression, driven by the Newman-Ziff sweep: * draw a uniformly random permutation of the L*L sites; * add sites one at a time, unioning each new site with whichever of its four neighbours are already present; * each root carries two bits, "this cluster touches the top row" and "this cluster touches the bottom row"; the bits are OR-ed on every union, so the first moment a root carries both bits is exactly the first moment the lattice spans; * record n_span, the number of occupied sites at that moment; * separately record S_max, the size of the largest cluster at the fixed occupancy n = round(p_c * L*L), which is the microcanonical critical point. Spanning is monotone: once a spanning cluster exists, adding more sites cannot destroy it. So a single sweep yields the whole microcanonical crossing curve R_L(n) = P(spanning | exactly n sites occupied) = P(n_span <= n), and the sweep can stop as soon as both n_span and the S_max checkpoint have been passed. The quantity actually wanted is the fixed-p crossing probability, obtained exactly by binomial convolution (the "canonical" ensemble): R_L(p) = sum_n C(N, n) p^n (1-p)^(N-n) R_L(n), N = L*L. That is an identity, not an approximation. Every sweep therefore contributes to every point of every crossing curve, which is the whole reason for using this sweep. The threshold estimate for a finite lattice is p*(L), defined by R_L(p*) = 1/2, and the extrapolation to infinite size uses the standard finite-size scaling form p*(L) = p_c + a L^(-1/nu), nu = 4/3 exactly in two dimensions. ASSUMPTIONS ----------- 1. Sites are independent. No correlation, no clustering, no lattice defects. 2. The correlation-length exponent nu = 4/3 is taken from the exactly solved 2D percolation universality class rather than fitted, for the headline number. We also fit it freely as a check and report what comes out. 3. The finite-size correction is taken to be a single power law. Corrections to scaling (a second, faster-decaying term) are ignored in the headline fit; we test the damage by refitting with the smallest lattices dropped. 4. numpy's PCG64 generator is treated as a source of independent uniform variates. Every stream is seeded from a stated master seed. 5. S_max is measured at fixed occupancy n = round(p_c N), not at fixed p. The two differ by finite-size corrections that vanish as L grows. LIMITATIONS, STATED PLAINLY --------------------------- * Sweep counts are not equal across L. The machine used here drops its clock by about a factor of three after a few seconds of sustained load, and CPython is orders of magnitude slower than the C codes that produced the reference value. We ran 40000 sweeps at L=16 and 320 at L=1024. The L=1024 point is the weakest number in this study and we say so rather than hide it. * Free boundaries were chosen because they are the easiest thing to define. Periodic boundaries in the transverse direction converge faster and are what most published codes use. That choice, not the arithmetic, is the largest single source of residual bias here. * Only top-to-bottom crossing is counted. Counting "either direction" changes the crossing probability at every finite L, though not the extrapolated threshold. * The published value 0.59274605 comes from simulations of the same kind, run far harder. This is a replication, not an independent measurement of a physical constant. VALIDATION PERFORMED -------------------- A. The union-find is checked against an independent breadth-first flood fill on thousands of random configurations at six lattice sizes. Spanning verdict and largest cluster size must agree exactly, every time. B. The incremental sweep is checked by rebuilding the configuration at n_span and at n_span - 1 and asking the flood fill whether each spans. The first must span, the second must not. C. For L=3 and L=4 every one of the 2^9 and 2^16 configurations is enumerated by brute force, giving the crossing probability as an exact polynomial in p. The Monte Carlo pipeline is compared against it end to end. D. The crossing probability evaluated at the accepted threshold is printed for every L to test whether it converges to a size-independent number (Cardy's formula gives exactly 1/2 for a square in the continuum limit). E. The extrapolated threshold is printed beside the accepted value with the difference in units of the club's own bootstrap error bar. F. The fractal dimension of the largest cluster at threshold is fitted and compared to 91/48. Seeds: master seed 20251215. Every worker stream is np.random.SeedSequence with entropy [20251215, L, chunk_index], so the output is reproducible regardless of how the work is scheduled across processes. Usage: python analysis/percolation-threshold.py > analysis/percolation-threshold-output.txt """ import os for _v in ("OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "MKL_NUM_THREADS", "NUMEXPR_NUM_THREADS", "VECLIB_MAXIMUM_THREADS"): os.environ.setdefault(_v, "1") import concurrent.futures as cf import math import sys import time import numpy as np # ---------------------------------------------------------------------------------- # constants # ---------------------------------------------------------------------------------- P_C_ACCEPTED = 0.59274605 # square-lattice site percolation, literature value NU = 4.0 / 3.0 # 2D percolation correlation-length exponent (exact) INV_NU = 1.0 / NU # 0.75 D_FRACTAL = 91.0 / 48.0 # 1.8958333..., exact in 2D MASTER_SEED = 20251215 EMPTY = -9 # sentinel in the union-find parent array BLK = 1 << 14 # site-order block size (keeps the int churn local) MAIN_SIZES = [(16, 40000), (32, 40000), (64, 20000), (128, 8000), (256, 5000), (512, 1200), (1024, 320)] EXACT_SIZES = [3, 4] EXACT_SWEEPS = 200000 N_BOOT = 400 N_BOOT_CONV = 100 WORKERS = 8 # ---------------------------------------------------------------------------------- # the union-find sweep # ---------------------------------------------------------------------------------- def sweep(L, order_np, n_target): """One Newman-Ziff sweep on an L x L lattice with free boundaries. order_np is a permutation of the padded interior cell indices. The lattice is stored with a one-cell border that is never occupied, so the four neighbours of cell s are s-1, s+1, s-W, s+W with no bounds tests. parent array `ptr`: EMPTY (-9) site not yet occupied >= 0 index of parent -1, -2, -3, -4 this site is a root; -(value)-1 is the two-bit flag (bit 0 = cluster touches top row, bit 1 = touches bottom row) Returns (n_span, smax) where n_span is the number of occupied sites when a spanning cluster first exists (-1 if the sweep ended first), and smax is the largest cluster size at exactly n_target occupied sites (-1 if that point was never reached). """ W = L + 2 NW = W * W W2 = 2 * W # first index below row 1 LW = L * W # first index of row L (the last real row) ptr = [EMPTY] * NW sz = [0] * NW nspan = -1 smax = -1 mx = 0 k = 0 ntot = order_np.shape[0] for b0 in range(0, ntot, BLK): for s in order_np[b0:b0 + BLK].tolist(): k += 1 if s < W2: ptr[s] = -2 # flags = 1, touches top row elif s >= LW: ptr[s] = -3 # flags = 2, touches bottom row else: ptr[s] = -1 # flags = 0 sz[s] = 1 if not mx: mx = 1 # a lone site is already a cluster of size one rs = s for j in (s - 1, s + 1, s - W, s + W): if ptr[j] != EMPTY: # find the root of the growing cluster, compressing as we go r = rs while ptr[r] >= 0: r = ptr[r] while ptr[rs] >= 0: nx = ptr[rs] ptr[rs] = r rs = nx rs = r # find the root of the neighbour, compressing as we go b = j while ptr[b] >= 0: b = ptr[b] r2 = b b = j while ptr[b] >= 0: nx = ptr[b] ptr[b] = r2 b = nx b = r2 if b != rs: if sz[rs] < sz[b]: # union by size rs, b = b, rs f = (-ptr[rs] - 1) | (-ptr[b] - 1) ptr[b] = rs ptr[rs] = -1 - f t = sz[rs] + sz[b] sz[rs] = t if t > mx: mx = t if f == 3 and nspan < 0: nspan = k # the checkpoint and the early exit are both taken only after every # neighbour of site k has been merged, so mx is complete if k == n_target: smax = mx if nspan > 0 and k >= n_target: return nspan, smax return nspan, smax def cells_of(L): """Padded indices of the L*L interior cells.""" W = L + 2 idx = np.arange(1, L + 1, dtype=np.int64) return ((idx[:, None]) * W + idx[None, :]).ravel() def chunk_worker(args): L, nsweeps, seed_key, n_target = args rng = np.random.default_rng(np.random.SeedSequence(list(seed_key))) cells = cells_of(L) ns = np.empty(nsweeps, dtype=np.int64) sm = np.empty(nsweeps, dtype=np.int64) for i in range(nsweeps): a, b = sweep(L, rng.permutation(cells), n_target) ns[i] = a sm[i] = b return ns, sm # ---------------------------------------------------------------------------------- # independent reference implementations, used only for validation # ---------------------------------------------------------------------------------- def uf_static(occ): """Spanning verdict and largest cluster for a fixed configuration, via the same union-find code the sweep uses (sites fed in index order).""" L = occ.shape[0] W = L + 2 rr, cc = np.nonzero(occ) if rr.size == 0: return False, 0 order = ((rr + 1) * W + (cc + 1)).astype(np.int64) nspan, smax = sweep(L, order, order.shape[0]) return (nspan > 0), smax def bfs_static(occ): """Spanning verdict and largest cluster by an independent flood fill.""" L = occ.shape[0] seen = np.zeros((L, L), dtype=bool) smax = 0 spans = False for r0 in range(L): for c0 in range(L): if occ[r0, c0] and not seen[r0, c0]: seen[r0, c0] = True stack = [(r0, c0)] size = 0 top = False bot = False while stack: y, x = stack.pop() size += 1 if y == 0: top = True if y == L - 1: bot = True if y > 0 and occ[y - 1, x] and not seen[y - 1, x]: seen[y - 1, x] = True stack.append((y - 1, x)) if y < L - 1 and occ[y + 1, x] and not seen[y + 1, x]: seen[y + 1, x] = True stack.append((y + 1, x)) if x > 0 and occ[y, x - 1] and not seen[y, x - 1]: seen[y, x - 1] = True stack.append((y, x - 1)) if x < L - 1 and occ[y, x + 1] and not seen[y, x + 1]: seen[y, x + 1] = True stack.append((y, x + 1)) if size > smax: smax = size if top and bot: spans = True return spans, smax def spans_mask(mask, L): """Does the configuration encoded in this bitmask span top to bottom?""" seen = 0 stack = [] for c in range(L): if (mask >> c) & 1: seen |= 1 << c stack.append(c) while stack: i = stack.pop() r, c = divmod(i, L) if r == L - 1: return True if r > 0: j = i - L if ((mask >> j) & 1) and not ((seen >> j) & 1): seen |= 1 << j stack.append(j) if r < L - 1: j = i + L if ((mask >> j) & 1) and not ((seen >> j) & 1): seen |= 1 << j stack.append(j) if c > 0: j = i - 1 if ((mask >> j) & 1) and not ((seen >> j) & 1): seen |= 1 << j stack.append(j) if c < L - 1: j = i + 1 if ((mask >> j) & 1) and not ((seen >> j) & 1): seen |= 1 << j stack.append(j) return False def exact_crossing_counts(L): """A[k] = number of configurations with exactly k occupied sites that span.""" N = L * L A = np.zeros(N + 1, dtype=np.float64) for mask in range(1 << N): if spans_mask(mask, L): A[bin(mask).count("1")] += 1.0 return A def exact_R(A, N, p): """A[k] already counts whole configurations, so each one carries probability p^k (1-p)^(N-k) with no binomial coefficient in front of it.""" ks = np.arange(N + 1, dtype=np.float64) lw = ks * math.log(p) + (N - ks) * math.log1p(-p) return float(np.sum(A * np.exp(lw))) # ---------------------------------------------------------------------------------- # microcanonical -> canonical conversion # ---------------------------------------------------------------------------------- class Convolver: """Exact binomial convolution of a microcanonical crossing curve, on a window of n wide enough that the binomial mass outside it is numerically negligible.""" def __init__(self, N, nspan, pmin, pmax): sd = math.sqrt(0.25 * N) n1 = int(math.floor(N * pmin - 14.0 * sd)) n2 = int(math.ceil(N * pmax + 14.0 * sd)) n1 = max(0, min(n1, int(np.min(nspan)))) n2 = min(N, max(n2, int(np.max(nspan)))) self.N = N self.n1 = n1 self.n2 = n2 n = np.arange(n1, n2 + 1, dtype=np.float64) self.n = n logC0 = (math.lgamma(N + 1) - math.lgamma(n1 + 1) - math.lgamma(N - n1 + 1)) if n.size > 1: inc = np.log(N - n[:-1]) - np.log(n[:-1] + 1.0) self.logC = logC0 + np.concatenate(([0.0], np.cumsum(inc))) else: self.logC = np.array([logC0]) self.width = n.size def cdf(self, nspan): """Empirical P(n_span <= n) on the window.""" counts = np.bincount(np.asarray(nspan) - self.n1, minlength=self.width) return np.cumsum(counts).astype(np.float64) / float(len(nspan)) def R(self, p, F): lw = self.logC + self.n * math.log(p) + (self.N - self.n) * math.log1p(-p) w = np.exp(lw) s = w.sum() return float(w.dot(F) / s), float(s) def solve_half(self, F, lo, hi, tol=1e-11): a, b = lo, hi if self.R(a, F)[0] - 0.5 > 0 or self.R(b, F)[0] - 0.5 < 0: return float("nan") while b - a > tol: m = 0.5 * (a + b) if self.R(m, F)[0] - 0.5 < 0.0: a = m else: b = m return 0.5 * (a + b) def wls(x, y, sig): """Weighted straight-line fit y = c + m x. Returns c, m, chi2, dof.""" w = 1.0 / np.asarray(sig) ** 2 S = w.sum() Sx = (w * x).sum() Sy = (w * y).sum() Sxx = (w * x * x).sum() Sxy = (w * x * y).sum() det = S * Sxx - Sx * Sx m = (S * Sxy - Sx * Sy) / det c = (Sxx * Sy - Sx * Sxy) / det resid = y - (c + m * x) chi2 = float((w * resid ** 2).sum()) return float(c), float(m), chi2, len(x) - 2 def hr(title): print() print("=" * 78) print(title) print("=" * 78) # ---------------------------------------------------------------------------------- def main(): t_start = time.time() print("SITE PERCOLATION ON THE SQUARE LATTICE: MEASURING p_c ON A LAPTOP") print("Science Journaling Club, Volume 1 Issue 2, Winter 2025") print("python %s | numpy %s" % (sys.version.split()[0], np.__version__)) print("master seed = %d | worker processes = %d" % (MASTER_SEED, WORKERS)) print("accepted value for square-lattice site percolation: p_c = %.8f" % P_C_ACCEPTED) # ------------------------------------------------------------------ validation A hr("VALIDATION A. UNION-FIND AGAINST AN INDEPENDENT FLOOD FILL") print("Random configurations, occupation probability drawn uniformly from") print("[0.30, 0.80]. Both the spanning verdict and the largest cluster size must") print("agree exactly. A single disagreement invalidates everything downstream.") print() print(" L configs span mismatches S_max mismatches") vrng = np.random.default_rng(np.random.SeedSequence([MASTER_SEED, 1])) tot_cfg = 0 tot_bad = 0 for L, reps in ((4, 4000), (6, 3000), (8, 2000), (12, 1200), (16, 800), (24, 400)): bad_s = 0 bad_m = 0 for _ in range(reps): p = float(vrng.uniform(0.30, 0.80)) occ = vrng.random((L, L)) < p s1, m1 = uf_static(occ) s2, m2 = bfs_static(occ) if s1 != s2: bad_s += 1 if m1 != m2: bad_m += 1 tot_cfg += reps tot_bad += bad_s + bad_m print("%4d %9d %17d %18d" % (L, reps, bad_s, bad_m)) print() print("total mismatches over %d configurations: %d" % (tot_cfg, tot_bad)) # ------------------------------------------------------------------ validation B hr("VALIDATION B. THE INCREMENTAL SWEEP LANDS ON THE EXACT CROSSING SITE") print("For each sweep we rebuild the configuration from the first n_span sites and") print("ask the flood fill whether it spans (it must), then remove the last site and") print("ask again (it must not). This tests the sweep, the flag bookkeeping and the") print("path compression together.") print() print(" L sweeps failures") for L, reps in ((8, 400), (16, 300), (24, 150)): W = L + 2 cells = cells_of(L) bad = 0 for _ in range(reps): order = vrng.permutation(cells) nspan, _ = sweep(L, order, L * L) if nspan < 0: bad += 1 continue idx = order[:nspan] rr = idx // W - 1 cc = idx % W - 1 occ = np.zeros((L, L), dtype=bool) occ[rr, cc] = True s_at, _ = bfs_static(occ) occ[rr[-1], cc[-1]] = False s_before, _ = bfs_static(occ) if not (s_at and not s_before): bad += 1 print("%4d %10d %10d" % (L, reps, bad)) tot_bad += bad print() print("Validations A and B together: %d failures." % tot_bad) if tot_bad != 0: print("STOPPING: the connectivity code is wrong, nothing below can be trusted.") return # ------------------------------------------------------------- the main campaign hr("MONTE CARLO CAMPAIGN") print("Newman-Ziff sweeps. One sweep = one uniformly random permutation of all") print("L*L sites, added one at a time, recording the occupancy at which the lattice") print("first spans and the largest cluster at occupancy round(p_c * L*L).") print("Every sweep contributes to every point of that lattice's crossing curve, so") print("the sweep count is the sample size at every p.") print() print(" L N sweeps n_target wall(s) sweeps/s") data = {} small = {} with cf.ProcessPoolExecutor(max_workers=WORKERS) as ex: for L, nsweeps in MAIN_SIZES: N = L * L n_target = int(round(P_C_ACCEPTED * N)) per = max(1, nsweeps // (WORKERS * 3)) tasks = [] left = nsweeps ci = 0 while left > 0: m = min(per, left) tasks.append((L, m, (MASTER_SEED, L, ci), n_target)) left -= m ci += 1 t0 = time.time() res = list(ex.map(chunk_worker, tasks)) dt = time.time() - t0 ns = np.concatenate([r[0] for r in res]) sm = np.concatenate([r[1] for r in res]) assert ns.min() > 0, "a sweep finished without ever spanning" assert sm.min() > 0, "S_max checkpoint never reached" data[L] = (ns, sm) print("%4d %9d %10d %10d %11.1f %10.1f" % (L, N, nsweeps, n_target, dt, nsweeps / dt)) # tiny lattices, for the exact-enumeration check for L in EXACT_SIZES: N = L * L per = max(1, EXACT_SWEEPS // (WORKERS * 3)) tasks = [] left = EXACT_SWEEPS ci = 0 while left > 0: m = min(per, left) tasks.append((L, m, (MASTER_SEED, 1000 + L, ci), N)) left -= m ci += 1 res = list(ex.map(chunk_worker, tasks)) small[L] = np.concatenate([r[0] for r in res]) # -------------------------------------------------------- microcanonical summary hr("MICROCANONICAL RESULT: OCCUPANCY AT FIRST SPANNING") print("n_span / N is the fraction of sites present when the lattice first crosses.") print("Its mean is itself an estimator of p_c, with the same finite-size correction") print("as the crossing-probability estimator but a different amplitude.") print() print(" L sweeps mean n_span/N s.d. std.err min/N max/N") micro = {} for L, _ in MAIN_SIZES: ns, _ = data[L] N = float(L * L) f = ns / N mu = float(f.mean()) sd = float(f.std(ddof=1)) se = sd / math.sqrt(len(f)) micro[L] = (mu, se, sd) print("%4d %10d %16.6f %10.6f %11.6f %8.4f %8.4f" % (L, len(f), mu, sd, se, f.min(), f.max())) # -------------------------------------------------------------- crossing curves hr("CANONICAL CROSSING CURVES R_L(p)") print("Exact binomial convolution of the microcanonical curve. The window of n used") print("for the convolution must carry essentially all of the binomial mass at every") print("p evaluated; the 'mass' column is the worst case seen for that L.") print() conv = {} grids = {} curves = {} worst_mass = {} for L, _ in MAIN_SIZES: ns, _ = data[L] N = L * L hw = min(0.34, max(0.006, 3.2 * L ** (-INV_NU))) pmin, pmax = P_C_ACCEPTED - hw, P_C_ACCEPTED + hw cv = Convolver(N, ns, pmin, pmax) F = cv.cdf(ns) grid = np.linspace(pmin, pmax, 41) vals = [] wm = 1.0 for p in grid: r, s = cv.R(float(p), F) vals.append(r) wm = min(wm, s) conv[L] = (cv, F) grids[L] = grid curves[L] = np.array(vals) worst_mass[L] = wm print(" L window width worst binomial mass in window R(pmin) R(pmax)") for L, _ in MAIN_SIZES: cv, F = conv[L] print("%4d %14d %31.12f %9.5f %9.5f" % (L, cv.width, worst_mass[L], curves[L][0], curves[L][-1])) # ------------------------------------------------------------------ validation C hr("VALIDATION C. BRUTE-FORCE ENUMERATION ON TINY LATTICES") print("Every configuration of a 3x3 and a 4x4 lattice is enumerated (512 and 65536") print("of them). That gives the crossing probability as an exact polynomial in p.") print("The Monte Carlo pipeline, run end to end on the same lattice sizes with") print("%d sweeps each, must reproduce it." % EXACT_SWEEPS) print() print(" L p exact R(p) Monte Carlo R(p) diff diff/s.e.") for L in EXACT_SIZES: N = L * L A = exact_crossing_counts(L) ns = small[L] cv = Convolver(N, ns, 0.02, 0.98) F = cv.cdf(ns) for p in (0.40, 0.50, P_C_ACCEPTED, 0.70, 0.80): ex_r = exact_R(A, N, p) mc_r, _ = cv.R(p, F) se = math.sqrt(max(ex_r * (1 - ex_r), 1e-12) / len(ns)) print("%4d %8.5f %14.9f %18.9f %11.2e %11.2f" % (L, p, ex_r, mc_r, mc_r - ex_r, (mc_r - ex_r) / se)) # ------------------------------------------------------------------ validation D hr("VALIDATION D. CROSSING PROBABILITY AT THE ACCEPTED THRESHOLD") print("If the crossing probability at p_c converges to a size-independent number,") print("that is the universal crossing probability. For a square (aspect ratio 1)") print("Cardy's formula gives exactly 1/2 in the continuum limit.") print() print(" L R_L(p_c) std.err R_L(p_c) - 0.5 in s.e.") r_at_pc = {} for L, _ in MAIN_SIZES: cv, F = conv[L] r, _ = cv.R(P_C_ACCEPTED, F) M = len(data[L][0]) se = math.sqrt(max(r * (1 - r), 1e-12) / M) r_at_pc[L] = (r, se) print("%4d %12.6f %11.6f %17.6f %9.2f" % (L, r, se, r - 0.5, (r - 0.5) / se)) # --------------------------------------------------- p*(L) and the bootstrap hr("FINITE-SIZE THRESHOLD ESTIMATES p*(L)") print("p*(L) is defined by R_L(p*) = 1/2, found by bisection on the convolved curve.") print("The error bar is a bootstrap over sweeps: %d resamples of the sweep set at" % N_BOOT) print("each L, with the whole convolution and bisection redone each time.") print() boot_rng = np.random.default_rng(np.random.SeedSequence([MASTER_SEED, 7])) pstar = {} pstar_se = {} boot_pstar = {} for L, _ in MAIN_SIZES: cv, F = conv[L] ns = data[L][0] lo, hi = float(grids[L][0]), float(grids[L][-1]) pstar[L] = cv.solve_half(F, lo, hi) M = len(ns) reps = np.empty(N_BOOT) for b in range(N_BOOT): idx = boot_rng.integers(0, M, M) reps[b] = cv.solve_half(cv.cdf(ns[idx]), lo, hi) boot_pstar[L] = reps pstar_se[L] = float(np.std(reps, ddof=1)) print(" L sweeps p*(L) boot s.e. L^(-1/nu) p*(L) - p_c") for L, _ in MAIN_SIZES: print("%4d %10d %13.7f %12.7f %11.6f %14.7f" % (L, len(data[L][0]), pstar[L], pstar_se[L], L ** (-INV_NU), pstar[L] - P_C_ACCEPTED)) # --------------------------------------------------------- finite-size scaling hr("FINITE-SIZE SCALING EXTRAPOLATION TO INFINITE LATTICE") print("Model: p*(L) = p_c + a * L^(-1/nu), with 1/nu = 3/4 held fixed at the exact") print("two-dimensional value. Weighted least squares, weights 1/s.e.^2.") print() Ls_all = [L for L, _ in MAIN_SIZES] def fit_set(Ls, pst, se, invnu=INV_NU): x = np.array([L ** (-invnu) for L in Ls]) y = np.array([pst[L] for L in Ls]) s = np.array([se[L] for L in Ls]) return wls(x, y, s) print(" lattices used p_c estimate slope a chi2 dof") fits = {} for drop in (0, 1, 2): Ls = Ls_all[drop:] c, m, chi2, dof = fit_set(Ls, pstar, pstar_se) fits[drop] = (Ls, c, m, chi2, dof) print(" L >= %-5d (%d points) %14.7f %14.6f %8.2f %5d" % (Ls[0], len(Ls), c, m, chi2, dof)) print() print("Bootstrapping the extrapolation: each of the %d replicates resamples the" % N_BOOT) print("sweeps at every L independently, re-solves p*(L) everywhere, and refits.") boot_pc = {} for drop in (0, 1, 2): Ls = fits[drop][0] x = np.array([L ** (-INV_NU) for L in Ls]) s = np.array([pstar_se[L] for L in Ls]) vals = np.empty(N_BOOT) for b in range(N_BOOT): y = np.array([boot_pstar[L][b] for L in Ls]) vals[b] = wls(x, y, s)[0] boot_pc[drop] = vals print() print(" lattices used p_c estimate boot s.e. 95% interval") for drop in (0, 1, 2): Ls = fits[drop][0] c = fits[drop][1] v = boot_pc[drop] se = float(np.std(v, ddof=1)) lo, hi = np.percentile(v, [2.5, 97.5]) print(" L >= %-5d %14.7f %13.7f [%.7f, %.7f]" % (Ls[0], c, se, lo, hi)) HEAD = 1 Ls_h = fits[HEAD][0] pc_hat = fits[HEAD][1] pc_se = float(np.std(boot_pc[HEAD], ddof=1)) hr("THE HEADLINE NUMBER") print("Club estimate (lattices L = %s, 1/nu = 3/4 fixed):" % ", ".join(str(L) for L in Ls_h)) print() print(" p_c (club) = %.7f +/- %.7f (1 s.e., bootstrap)" % (pc_hat, pc_se)) print(" p_c (accepted) = %.8f" % P_C_ACCEPTED) print(" difference = %+.7f" % (pc_hat - P_C_ACCEPTED)) print(" difference = %+.2f of the club's own error bars" % ((pc_hat - P_C_ACCEPTED) / pc_se)) print() print(" relative precision of the club estimate: %.2e" % (pc_se / pc_hat)) print(" decimal places agreed with the accepted value: %d" % max(0, int(-math.log10(abs(pc_hat - P_C_ACCEPTED))))) print() print("For comparison, the same arithmetic applied to the microcanonical estimator") print("(mean occupancy at first spanning) rather than the R = 1/2 criterion:") mic_val = {L: micro[L][0] for L in Ls_all} mic_se = {L: micro[L][1] for L in Ls_all} cmi, mmi, chimi, dofmi = fit_set(Ls_h, mic_val, mic_se) print(" p_c (club, mean n_span/N route) = %.7f difference %+.7f" % (cmi, cmi - P_C_ACCEPTED)) print(" slope a = %+.6f, chi2/dof = %.1f/%d" % (mmi, chimi, dofmi)) # ---------------------------------------------- the correlation length exponent hr("LETTING THE EXPONENT FLOAT") print("Same data, same fit, but now 1/nu is scanned instead of fixed. The") print("chi-square is minimised on a grid; the reported nu is 1/(best exponent).") print() ys = np.arange(0.30, 1.6001, 0.002) best = None for yy in ys: c, m, chi2, dof = fit_set(Ls_h, pstar, pstar_se, invnu=float(yy)) if best is None or chi2 < best[3]: best = (float(yy), c, m, chi2) print(" best 1/nu = %.3f (exact value 0.750)" % best[0]) print(" implied nu = %.3f (exact value %.4f)" % (1.0 / best[0], NU)) print(" p_c at that exponent = %.7f" % best[1]) print(" chi2 at minimum = %.3f" % best[3]) print() print(" The chi-square surface is shallow in this direction. A laptop-sized data") print(" set constrains p_c far better than it constrains nu, because p_c is the") print(" intercept and nu only reshapes the approach to it.") print() print(" 1/nu p_c from that exponent chi2") for yy in (0.60, 0.65, 0.70, 0.75, 0.80, 0.85, 0.90, 1.00): c, m, chi2, dof = fit_set(Ls_h, pstar, pstar_se, invnu=yy) print(" %.2f %22.7f %11.3f" % (yy, c, chi2)) # ---------------------------------------------------------- fractal dimension hr("VALIDATION F. LARGEST CLUSTER AT THRESHOLD") print("S_max is measured at fixed occupancy n = round(p_c * N). Theory says") print("S_max ~ L^D with D = 91/48 = %.7f exactly." % D_FRACTAL) print() print(" L mean S_max std.err S_max/N S_max / L^(91/48)") xs = [] ysv = [] sgs = [] for L, _ in MAIN_SIZES: sm = data[L][1].astype(np.float64) mu = float(sm.mean()) se = float(sm.std(ddof=1) / math.sqrt(len(sm))) xs.append(math.log(L)) ysv.append(math.log(mu)) sgs.append(se / mu) print("%4d %15.2f %11.2f %11.6f %21.6f" % (L, mu, se, mu / (L * L), mu / L ** D_FRACTAL)) c, m, chi2, dof = wls(np.array(xs), np.array(ysv), np.array(sgs)) print() print(" log-log weighted fit over all seven sizes:") print(" D (club) = %.5f" % m) print(" D (exact) = %.7f (91/48)" % D_FRACTAL) print(" difference = %+.5f" % (m - D_FRACTAL)) print(" chi2 / dof = %.1f / %d" % (chi2, dof)) c2, m2, chi22, dof2 = wls(np.array(xs[2:]), np.array(ysv[2:]), np.array(sgs[2:])) print(" dropping L = 16 and 32, where corrections to scaling bite hardest:") print(" D (club) = %.5f difference %+.5f chi2/dof %.1f/%d" % (m2, m2 - D_FRACTAL, chi22, dof2)) # ------------------------------------------------------------------ convergence hr("CONVERGENCE: WHAT THE ANSWER LOOKED LIKE AFTER m SWEEPS") print("Using only the first m sweeps at every lattice size (capped at what that") print("size actually has), the whole pipeline is re-run: convolution, bisection,") print("refit. The error bar is a %d-replicate bootstrap of that reduced data set." % N_BOOT_CONV) print() ladder = [50, 100, 200, 320, 500, 1000, 2000, 5000, 20000, 40000] print(" m p*(256) at m p_c(m) [L>=32] boot s.e. p_c(m) - accepted") conv_rows = [] crng = np.random.default_rng(np.random.SeedSequence([MASTER_SEED, 9])) for m in ladder: pst_m = {} se_m = {} bt_m = {} for L in Ls_h + [256]: if L in pst_m: continue ns = data[L][0] mm = min(m, len(ns)) sub = ns[:mm] cv = conv[L][0] lo, hi = float(grids[L][0]), float(grids[L][-1]) pst_m[L] = cv.solve_half(cv.cdf(sub), lo, hi) reps = np.empty(N_BOOT_CONV) for b in range(N_BOOT_CONV): idx = crng.integers(0, mm, mm) reps[b] = cv.solve_half(cv.cdf(sub[idx]), lo, hi) bt_m[L] = reps se_m[L] = float(np.std(reps, ddof=1)) x = np.array([L ** (-INV_NU) for L in Ls_h]) s = np.array([se_m[L] for L in Ls_h]) y = np.array([pst_m[L] for L in Ls_h]) cm = wls(x, y, s)[0] bvals = np.array([wls(x, np.array([bt_m[L][b] for L in Ls_h]), s)[0] for b in range(N_BOOT_CONV)]) cse = float(np.std(bvals, ddof=1)) conv_rows.append((m, pst_m[256], cm, cse, cm - P_C_ACCEPTED)) print("%7d %14.7f %16.7f %13.7f %20.7f" % (m, pst_m[256], cm, cse, cm - P_C_ACCEPTED)) # -------------------------------------------------------- machine-readable block hr("FIGURE DATA (machine readable)") print("# CURVE L then 41 pairs p,R") for L, _ in MAIN_SIZES: g = grids[L] c_ = curves[L] print("CURVE %d %s" % (L, " ".join("%.6f,%.5f" % (a, b) for a, b in zip(g, c_)))) print("# SIZE L sweeps pstar se meanNspanFrac se R_at_pc se meanSmax se") for L, _ in MAIN_SIZES: sm = data[L][1].astype(np.float64) print("SIZE %d %d %.7f %.7f %.7f %.7f %.6f %.6f %.3f %.3f" % (L, len(data[L][0]), pstar[L], pstar_se[L], micro[L][0], micro[L][1], r_at_pc[L][0], r_at_pc[L][1], sm.mean(), sm.std(ddof=1) / math.sqrt(len(sm)))) print("# CONV m pstar256 pc se pc-accepted") for row in conv_rows: print("CONV %d %.7f %.7f %.7f %.7f" % row) print("# FIT pc se slope chi2 dof Lmin") print("FIT %.7f %.7f %.6f %.3f %d %d" % (pc_hat, pc_se, fits[HEAD][2], fits[HEAD][3], fits[HEAD][4], Ls_h[0])) print("# BOOTHIST 30 bins lo,hi,count over the bootstrap p_c replicates") hist, edges = np.histogram(boot_pc[HEAD], bins=30) print("BOOTHIST " + " ".join("%.7f,%.7f,%d" % (edges[i], edges[i + 1], hist[i]) for i in range(30))) print("# NSPAN L then 40 bins lo,count of n_span/N over [0.40,0.80]") for L, _ in MAIN_SIZES: f = data[L][0] / float(L * L) h, e = np.histogram(f, bins=40, range=(0.40, 0.80)) print("NSPAN %d %s" % (L, " ".join("%.4f,%d" % (e[i], h[i]) for i in range(40)))) print() print("total runtime: %.1f s" % (time.time() - t_start)) if __name__ == "__main__": main()