race.py
12.5 kB · python · 356 lines
1import statistics2import time3from collections import deque4from itertools import product5from random import Random67import numpy as np89NUMPY_SEED = 2026072510PYTHON_SEED_BASE = 1000111213def corners(base, dimension, drop):14 return [c for c in product(range(base), repeat=dimension) if not drop(c)]151617DESIGNS = {18 "gasket": (2, 2, [(0, 0), (0, 1), (1, 0)]),19 "diagonal": (2, 2, [(0, 0), (1, 1)]),20 "antidiagonal": (2, 2, [(0, 1), (1, 0)]),21 "seven-of-eight": (2, 3, corners(2, 3, lambda c: c == (1, 1, 1))),22 "carpet": (3, 2, corners(3, 2, lambda c: c == (1, 1))),23 "sponge": (3, 3, corners(3, 3, lambda c: c.count(1) > 1)),24}2526SAME = {"antidiagonal": "diagonal"}2728PASS_A = [29 ("gasket", 5, 400), ("gasket", 6, 400),30 ("antidiagonal", 5, 400), ("antidiagonal", 6, 400),31 ("seven-of-eight", 3, 400), ("seven-of-eight", 4, 400),32 ("carpet", 3, 400), ("carpet", 4, 200),33 ("sponge", 3, 200), ("sponge", 4, 25),34]3536PASS_B = [37 ("gasket", 5, 400), ("gasket", 6, 400), ("gasket", 7, 200), ("gasket", 8, 100),38 ("diagonal", 5, 400), ("diagonal", 6, 400), ("diagonal", 7, 200),39 ("seven-of-eight", 3, 400), ("seven-of-eight", 4, 200),40 ("carpet", 3, 400), ("carpet", 4, 200),41 ("sponge", 3, 100), ("sponge", 4, 20),42]434445def kron_power(base, dimension, cells, level):46 tile = np.zeros((base,) * dimension, dtype=np.uint8)47 for cell in cells:48 tile[cell] = 149 out = tile50 for _ in range(level - 1):51 out = np.kron(out, tile)52 return out.astype(bool)535455def row_strides(shape):56 out = [1] * len(shape)57 for axis in range(len(shape) - 2, -1, -1):58 out[axis] = out[axis + 1] * shape[axis + 1]59 return out606162def edge_ranks(occupied, grid, cells):63 shape = occupied.shape64 stride = row_strides(shape)65 lows = []66 for axis in range(len(shape)):67 lo = [slice(None)] * len(shape)68 hi = [slice(None)] * len(shape)69 lo[axis] = slice(0, shape[axis] - 1)70 hi[axis] = slice(1, shape[axis])71 both = occupied[tuple(lo)] & occupied[tuple(hi)]72 left = grid[tuple(lo)][both]73 lows.append((left, left + stride[axis]))74 a = np.concatenate([pair[0] for pair in lows])75 b = np.concatenate([pair[1] for pair in lows])76 return np.searchsorted(cells, a), np.searchsorted(cells, b)777879def union_stats(n, ea, eb):80 parent = list(range(n))81 weight = [1] * n82 count = n83 largest = 184 for i, j in zip(ea, eb):85 while parent[i] != i:86 parent[i] = parent[parent[i]]87 i = parent[i]88 while parent[j] != j:89 parent[j] = parent[parent[j]]90 j = parent[j]91 if i == j:92 continue93 if weight[i] < weight[j]:94 i, j = j, i95 parent[j] = i96 weight[i] += weight[j]97 count -= 198 if weight[i] > largest:99 largest = weight[i]100 return count, largest101102103def kron_stats(occupied, grid, dimension):104 cells = np.flatnonzero(occupied.ravel())105 n = int(cells.size)106 ranks_a, ranks_b = edge_ranks(occupied, grid, cells)107 edges = int(ranks_a.size)108 count, largest = union_stats(n, ranks_a.tolist(), ranks_b.tolist())109 return n, count, largest / n, (2 * dimension * n - 2 * edges) / n110111112def substitution(base, dimension, cells, level):113 out = [(0,) * dimension]114 for _ in range(level):115 out = [tuple(c[a] * base + f[a] for a in range(dimension))116 for c in out for f in cells]117 return out118119120def digit_rule(base, dimension, cells, level):121 keep = set(cells)122 side = base ** level123 out = []124 for index in range(side ** dimension):125 coord = []126 left = index127 for _ in range(dimension):128 coord.append(left % side)129 left //= side130 ok = True131 for position in range(level):132 if tuple((c // base ** position) % base for c in coord) not in keep:133 ok = False134 break135 if ok:136 out.append(index)137 return out138139140def padded_steps(side, dimension):141 return [(side + 2) ** a for a in range(dimension)]142143144def place(coords, side, dimension, step):145 occupied = bytearray((side + 2) ** dimension)146 flat = [sum((c[a] + 1) * step[a] for a in range(dimension)) for c in coords]147 for index in flat:148 occupied[index] = 1149 return occupied, flat150151152def spread_out(picks, side, dimension, step):153 out = []154 for index in picks:155 left = index156 total = 0157 for a in range(dimension):158 total += (left % side + 1) * step[a]159 left //= side160 out.append(total)161 return out162163164def bfs_stats(occupied, flat, step, dimension):165 n = len(flat)166 edges = 0167 for index in flat:168 for s in step:169 if occupied[index + s]:170 edges += 1171 moves = step + [-s for s in step]172 seen = bytearray(len(occupied))173 components = 0174 largest = 0175 for start in flat:176 if seen[start]:177 continue178 components += 1179 seen[start] = 1180 queue = deque([start])181 size = 0182 while queue:183 index = queue.popleft()184 size += 1185 for s in moves:186 neighbour = index + s187 if occupied[neighbour] and not seen[neighbour]:188 seen[neighbour] = 1189 queue.append(neighbour)190 if size > largest:191 largest = size192 return n, components, largest / n, (2 * dimension * n - 2 * edges) / n193194195def run_a(name, level, seeds):196 base, dimension, cells = DESIGNS[name]197 design = kron_power(base, dimension, cells, level)198 grid = np.arange(design.size, dtype=np.int64).reshape(design.shape)199 n, dc, df, db = kron_stats(design, grid, dimension)200 rng = np.random.default_rng(NUMPY_SEED)201 comps, frac, bound = [], [], []202 for _ in range(seeds):203 flat = np.zeros(design.size, dtype=bool)204 flat[rng.choice(design.size, n, replace=False)] = True205 _, c, f, b = kron_stats(flat.reshape(design.shape), grid, dimension)206 comps.append(float(c))207 frac.append(f)208 bound.append(b)209 return pack(name, level, base ** level, dimension, n, design.size, seeds,210 dc, df, db, comps, frac, bound,211 lambda v: (float(np.mean(v)), float(np.std(v, ddof=1))))212213214def draw_b(name, level, seeds):215 base, dimension, cells = DESIGNS[name]216 side = base ** level217 total = side ** dimension218 step = padded_steps(side, dimension)219 design = substitution(base, dimension, cells, level)220 occupied, flat = place(design, side, dimension, step)221 n, dc, df, db = bfs_stats(occupied, flat, step, dimension)222 comps, frac, bound = [], [], []223 for seed in range(seeds):224 picks = Random(PYTHON_SEED_BASE + seed).sample(range(total), n)225 holes = spread_out(picks, side, dimension, step)226 blank = bytearray((side + 2) ** dimension)227 for index in holes:228 blank[index] = 1229 _, c, f, b = bfs_stats(blank, holes, step, dimension)230 comps.append(float(c))231 frac.append(f)232 bound.append(b)233 return pack(name, level, side, dimension, n, total, seeds,234 dc, df, db, comps, frac, bound,235 lambda v: (statistics.mean(v), statistics.stdev(v)))236237238def pack(name, level, side, dimension, n, total, seeds, dc, df, db,239 comps, frac, bound, moments):240 return {241 "name": name, "level": level, "side": side, "dimension": dimension,242 "cells": n, "density": n / total, "seeds": seeds,243 "design": (dc, df, db),244 "comps": moments(comps), "frac": moments(frac), "bound": moments(bound),245 "tie_comps": sum(1 for v in comps if v <= dc),246 "tie_bound": sum(1 for v in bound if v >= db),247 }248249250def show(result):251 grid = "x".join([str(result["side"])] * result["dimension"])252 print(" {:<15} L={} grid {:<12} cells {:>7} density {:.4f} seeds {}".format(253 result["name"], result["level"], grid, result["cells"],254 result["density"], result["seeds"]))255 dc, df, db = result["design"]256 print(" components design {:>10} random {:12.4f} +/- {:.4f} ties {}/{}".format(257 dc, result["comps"][0], result["comps"][1],258 result["tie_comps"], result["seeds"]))259 print(" largest design {:>10.4f} random {:12.4f} +/- {:.4f}".format(260 df, result["frac"][0], result["frac"][1]))261 print(" boundary design {:>10.4f} random {:12.4f} +/- {:.4f} ties {}/{}".format(262 db, result["bound"][0], result["bound"][1],263 result["tie_bound"], result["seeds"]))264265266def main():267 start = time.time()268 print("DOMAIN")269 print(" pass A: 2D grids to 81x81, 3D to 81^3, 400 seeds thinning to 25")270 print(" pass B: 2D grids to 256x256, 3D to 81^3, 400 seeds thinning to 20")271 print(" ties: random draws at or below the design count, at or above the design boundary")272 print()273274 print("SUBSTITUTION AGAINST DIGIT RULE")275 agree = True276 for name in ("gasket", "diagonal", "seven-of-eight", "carpet", "sponge"):277 base, dimension, cells = DESIGNS[name]278 for level in (1, 2, 3):279 side = base ** level280 a = sorted(digit_rule(base, dimension, cells, level))281 raw = substitution(base, dimension, cells, level)282 b = sorted({sum(c[x] * side ** x for x in range(dimension)) for c in raw})283 same = a == b and len(a) == len(cells) ** level284 agree = agree and same285 print(" {:<15} L={} cells {:>7} routes agree {}".format(286 name, level, len(a), same))287 print(" every route agrees and the fill count is multiplicative: {}".format(agree))288 print()289290 print("PASS A - KRONECKER POWERS, UNION-FIND, NUMPY PCG64")291 a_results = {}292 for name, level, seeds in PASS_A:293 result = run_a(name, level, seeds)294 show(result)295 a_results[(SAME.get(name, name), level)] = result296 print()297298 print("PASS B - SUBSTITUTION, BREADTH-FIRST SEARCH, PYTHON MERSENNE TWISTER")299 b_results = {}300 for name, level, seeds in PASS_B:301 result = draw_b(name, level, seeds)302 show(result)303 b_results[(SAME.get(name, name), level)] = result304 print()305306 print("DISPERSING EXTREME - COMPONENTS ARE NOT THE METRIC")307 for level in (5, 6, 7):308 result = b_results[("diagonal", level)]309 print(" diagonal L={} side {:>4} cells {:>5} seeds {:>3} design comps {:>5} random comps {:10.4f} +/- {:.4f}".format(310 level, result["side"], result["cells"], result["seeds"],311 result["design"][0], result["comps"][0], result["comps"][1]))312 print()313314 print("MAXIMUM BOUNDARY")315 for key in (("diagonal", 6), ("gasket", 6), ("carpet", 4),316 ("sponge", 3), ("seven-of-eight", 4)):317 result = b_results[key]318 print(" {:<15} L={} boundary per cell {:.4f} of the maximum {}".format(319 key[0], key[1], result["design"][2], 2 * result["dimension"]))320 print()321322 print("THE TWO PASSES AGAINST EACH OTHER")323 worst = 0.0324 for key in sorted(set(a_results) & set(b_results)):325 left, right = a_results[key], b_results[key]326 for field in ("comps", "frac", "bound"):327 gap = abs(left[field][0] - right[field][0])328 sd = max(left[field][1], right[field][1])329 worst = max(worst, gap / sd if sd else 0.0)330 print(" {:<15} L={} {:<6} pass A {:12.4f} pass B {:12.4f} gap {:.4f} sd".format(331 key[0], key[1], field, left[field][0], right[field][0],332 gap / sd if sd else 0.0))333 print(" every comparison agrees to within one standard deviation: {}".format(worst <= 1.0))334 print(" the widest gap is {:.4f} standard deviations".format(worst))335 print()336337 print("WITNESSES")338 gasket8 = b_results[("gasket", 8)]339 gasket5 = b_results[("gasket", 5)]340 sponge4 = b_results[("sponge", 4)]341 print(" gasket 256x256 design comps {} random {:.2f} +/- {:.2f}".format(342 gasket8["design"][0], gasket8["comps"][0], gasket8["comps"][1]))343 print(" gasket 256x256 design boundary {:.4f} random {:.4f} +/- {:.4f}".format(344 gasket8["design"][2], gasket8["bound"][0], gasket8["bound"][1]))345 print(" sponge 81^3 design comps {} random {:.2f} +/- {:.2f}".format(346 sponge4["design"][0], sponge4["comps"][0], sponge4["comps"][1]))347 print(" gasket 32x32 random draws reaching one component {}/{}".format(348 gasket5["tie_comps"], gasket5["seeds"]))349 print()350 print(" every check passed: {}".format(agree and worst <= 1.0))351 print(" wall {:.1f} s".format(time.time() - start))352 if not agree or worst > 1.0:353 raise SystemExit(1)354355356main()