waves.py
17.3 kB · python · 381 lines
1import itertools2import sys3from collections import Counter4from fractions import Fraction5from math import ceil, floor67import numpy as np89N = 25610SEED = 2026011SEEDS = (20260, 4093)12STEPS = 6413DENSITY = 0.514KMAX = N // 215DEAD = 0.0216FULL = 0.9817GAIN = 3.018TOL = 1e-91920BUGS_B = (Fraction(34, 120), Fraction(45, 120))21BUGS_S = (Fraction(34, 120), Fraction(58, 120))22RULES = [23 ("bugs", BUGS_B, BUGS_S),24 ("wide-survive", BUGS_B, (Fraction(28, 100), Fraction(60, 100))),25 ("narrow-birth", (Fraction(30, 100), Fraction(34, 100)), BUGS_S),26]27GRID_B_LO = [Fraction(20, 100), Fraction(30, 100), Fraction(35, 100)]28GRID_B_HI = [Fraction(35, 100), Fraction(40, 100), Fraction(45, 100), Fraction(50, 100), Fraction(55, 100)]29GRID_S_LO = [Fraction(0), Fraction(30, 100)]30GRID_S_HI = [Fraction(45, 100), Fraction(50, 100), Fraction(55, 100), Fraction(60, 100), Fraction(70, 100), Fraction(80, 100)]313233def residue_corners(d, base=2):34 return [[(i // base ** (d - 1 - j)) % base for j in range(d)] for i in range(base ** d)]353637def tile(code, side=3):38 corners = residue_corners(2)39 filled = [(code >> i) & 1 for i in range(len(corners))]40 out = np.zeros((side, side), dtype=np.uint8)41 for x in range(side):42 for y in range(side):43 out[x, y] = filled[(x % 2) * 2 + (y % 2)]44 return out454647def kron_power(t, level):48 out = t49 for _ in range(1, level):50 out = np.kron(out, t)51 return out525354def design_mask(code, side, level):55 mask = kron_power(tile(code, side), level).copy()56 c = (side ** level - 1) // 257 mask[c, c] = 058 return mask596061def box_mask(r):62 mask = np.ones((2 * r + 1, 2 * r + 1), dtype=np.uint8)63 mask[r, r] = 064 return mask656667def torus_kernel(mask):68 side = mask.shape[0]69 c = (side - 1) // 270 kernel = np.zeros((N, N))71 xs, ys = np.nonzero(mask)72 kernel[(xs - c) % N, (ys - c) % N] = 173 return kernel747576FX = np.fft.fftfreq(N) * N77KX, KY = np.meshgrid(FX, FX, indexing="ij")78RING = np.rint(np.sqrt(KX ** 2 + KY ** 2)).astype(int)79COUNT = np.bincount(RING.ravel())808182def ring_sum(values):83 return np.bincount(RING.ravel(), weights=values.ravel(), minlength=COUNT.size)848586def ring_mean(values):87 return ring_sum(values) / COUNT888990def first_min(a, start):91 return next((k for k in range(start, KMAX) if a[k + 1] >= a[k]), None)929394def first_max(a, start):95 return next((k for k in range(start, KMAX) if a[k + 1] < a[k]), None)969798def mask_profile(full):99 a = ring_mean(np.abs(full))[: KMAX + 1]100 k_min = first_min(a, 1)101 k_2 = first_max(a, k_min + 1)102 k_min2 = first_min(a, k_2 + 1)103 half = a[k_2] / 2104 lobe = [k for k in range(k_min, k_min2 + 1) if a[k] >= half]105 g = ring_mean(np.real(full))[: KMAX + 1]106 k_neg = 1 + int(np.argmin(g[1:]))107 lo = k_neg108 while lo > 1 and g[lo - 1] < 0:109 lo -= 1110 hi = k_neg111 while hi < KMAX and g[hi + 1] < 0:112 hi += 1113 return a, k_min, k_2, k_min2, (lobe[0], lobe[-1]), g, k_neg, (lo, hi)114115116def riesz_product(code, side, level):117 t = tile(code, side)118 c = (side - 1) // 2119 xs, ys = np.nonzero(t)120 ex, ey = xs - c, ys - c121 out = np.ones((N, N), dtype=complex)122 for j in range(level):123 tx = side ** j * KX / N124 ty = side ** j * KY / N125 p = np.zeros((N, N), dtype=complex)126 for a, b in zip(ex, ey):127 p += np.exp(-2j * np.pi * (tx * a + ty * b))128 out *= p129 return out - int(t[c, c])130131132def thresholds(m, lo, hi):133 return ceil(lo * m), floor(hi * m)134135136def soup(seed):137 rng = np.random.default_rng(seed)138 return (rng.random((N, N)) < DENSITY).astype(np.uint8)139140141def evolve(half, m, birth, survive, seed=SEED):142 blo, bhi = thresholds(m, *birth)143 slo, shi = thresholds(m, *survive)144 grid = soup(seed)145 churn = None146 for _ in range(STEPS):147 count = np.rint(np.fft.irfft2(np.fft.rfft2(grid) * half, s=(N, N))).astype(int)148 born = (grid == 0) & (count >= blo) & (count <= bhi)149 kept = (grid == 1) & (count >= slo) & (count <= shi)150 new = (born | kept).astype(np.uint8)151 churn = int(np.count_nonzero(new != grid))152 grid = new153 return grid, churn154155156def spectrum(grid):157 f = np.fft.rfft2(grid - grid.mean())158 power = np.zeros((N, N))159 power[:, : N // 2 + 1] = np.abs(f) ** 2160 power[:, N // 2 + 1 :] = np.abs(f[np.r_[0, N - 1 : 0 : -1], N // 2 - 1 : 0 : -1]) ** 2161 s = ring_sum(power)[: KMAX + 1]162 d = s / COUNT[: KMAX + 1]163 if s[1:].sum() == 0:164 return 1, 0.0, 0.0165 k = 1 + int(np.argmax(d[1:]))166 share = s[k] / s[1:].sum()167 gain = d[k] / d[1:].mean()168 return k, share, gain169170171def classify(density, churn, k, gain, side):172 if density < DEAD:173 return "dead"174 if density > FULL:175 return "full"176 if churn > 0:177 return "active"178 if gain < GAIN:179 return "flat"180 if N / k > 2 * side:181 return "coarse"182 return "ring"183184185def fmt(fr):186 return f"{float(fr):.2f}"187188189def rule_label(birth, survive):190 return f"B[{fmt(birth[0])},{fmt(birth[1])}] S[{fmt(survive[0])},{fmt(survive[1])}]"191192193def main():194 t7 = tile(7)195 assert t7.sum() == 8 and t7[1, 1] == 0 and t7.shape == (3, 3)196 m7 = design_mask(7, 3, 2)197 assert m7.shape == (9, 9) and m7.sum() == 64198 assert design_mask(6, 3, 1).sum() == 4 and design_mask(6, 3, 1)[0, 1] == 1199 assert design_mask(9, 3, 1).sum() == 4 and design_mask(9, 3, 1)[0, 0] == 1200 print(f"domain: torus {N}x{N}, soup density {DENSITY}, seed {SEED}, {STEPS} steps, rings k = 1..{KMAX}")201 print(f"classes: dead is density < {DEAD}, full is density > {FULL}, active is churn > 0 at the last step, flat is gain < {GAIN}, coarse is N/k* > 2 side, ring is the rest")202 for seed in SEEDS:203 k0, share0, gain0 = spectrum(soup(seed))204 print(f"soup seed {seed} at step 0: k* {k0}, share {share0:.4f}, gain {gain0:.2f}")205 print()206 print("table 1: the Riesz product, max |DFT(mask) - (prod_j P(side^j xi / N) - tile centre)| over all frequencies")207 print(f"{'code':>4} {'level':>5} {'side':>4} {'ones':>5} {'centre':>6} {'max err':>10} {'tag':>8}")208 riesz_ok = True209 for code in (7, 6, 9):210 for level in (2, 3, 4):211 mask = design_mask(code, 3, level)212 err = np.max(np.abs(np.fft.fft2(torus_kernel(mask)) - riesz_product(code, 3, level)))213 ok = err < TOL214 riesz_ok &= ok215 print(f"{code:>4} {level:>5} {3 ** level:>4} {mask.sum():>5} {int(tile(code)[1, 1]):>6} {err:>10.2e} {'Verified' if ok else 'Refuted':>8}")216 print(f"Riesz product at codes 7, 6, 9, levels 2..4, tolerance {TOL:.0e}: {'Verified' if riesz_ok else 'Refuted'}")217 print()218 masks = [219 ("box r=5", box_mask(5), 11),220 ("box r=13", box_mask(13), 27),221 ("code 7 L2", design_mask(7, 3, 2), 9),222 ("code 7 L3", design_mask(7, 3, 3), 27),223 ("code 7 L4", design_mask(7, 3, 4), 81),224 ("code 6 L3", design_mask(6, 3, 3), 27),225 ("code 9 L3", design_mask(9, 3, 3), 27),226 ]227 print("mask table: ones m, first minimum k_min, first secondary maximum k_2, second minimum, half-height lobe of the ring-mean |DFT|, then the most negative ring k_neg and the negative band of the ring-mean signed DFT")228 print(f"{'mask':>10} {'side':>4} {'m':>5} {'k_min':>5} {'k_2':>5} {'k_min2':>6} {'lobe':>9} {'N/k_min':>8} {'N/k_2':>7} {'N/k_2/side':>10} {'A(k_min)/m':>10} {'A(k_2)/m':>9} {'k_neg':>5} {'neg band':>9} {'g(k_neg)/m':>10}")229 profiles = {}230 for name, mask, side in masks:231 kernel = torus_kernel(mask)232 m = int(mask.sum())233 a, k_min, k_2, k_min2, lobe, g, k_neg, neg = mask_profile(np.fft.fft2(kernel))234 profiles[name] = (np.fft.rfft2(kernel), m, k_min, k_2, lobe, side, g, k_neg, neg)235 print(f"{name:>10} {side:>4} {m:>5} {k_min:>5} {k_2:>5} {k_min2:>6} {f'{lobe[0]}..{lobe[1]}':>9} {N / k_min:>8.2f} {N / k_2:>7.2f} {N / k_2 / side:>10.3f} {a[k_min] / m:>10.4f} {a[k_2] / m:>9.4f} {k_neg:>5} {f'{neg[0]}..{neg[1]}':>9} {g[k_neg] / m:>10.4f}")236 print()237 results = {}238 print("run table: the state at the last step for every mask under the three named rules")239 print(f"{'mask':>10} {'rule':>13} {'B':>9} {'S':>9} {'density':>7} {'churn':>6} {'k*':>3} {'N/k*':>7} {'share':>6} {'gain':>6} {'class':>6}")240 for name, mask, side in masks:241 half, m, k_min, k_2, lobe, side, g, k_neg, neg = profiles[name]242 for rule, birth, survive in RULES:243 grid, churn = evolve(half, m, birth, survive)244 density = grid.mean()245 k, share, gain = spectrum(grid)246 cls = classify(density, churn, k, gain, side)247 results[(name, rule)] = (k, share, gain, density, churn, cls)248 b = "{}-{}".format(*thresholds(m, *birth))249 s = "{}-{}".format(*thresholds(m, *survive))250 print(f"{name:>10} {rule:>13} {b:>9} {s:>9} {density:>7.4f} {churn:>6} {k:>3} {N / k:>7.2f} {share:>6.3f} {gain:>6.1f} {cls:>6}")251 print()252 print("table 2: same mask across the three named rules, k* of every ring still")253 print(f"{'mask':>10} {'bugs':>5} {'wide':>5} {'narrow':>6} {'rings':>5}")254 rings_named = 0255 for name, mask, side in masks:256 ks = [results[(name, rule)][0] if results[(name, rule)][5] == "ring" else None for rule, _, _ in RULES]257 cells = " ".join(f"{k:>5}" if k is not None else f"{'-':>5}" for k in ks)258 kept = [k for k in ks if k is not None]259 rings_named += len(kept)260 print(f"{name:>10} {cells} {len(kept):>5}")261 print(f"ring stills under the three named rules: {rings_named} of {len(masks) * len(RULES)}; the wavelength laws are vacuous on the named rules, tested on the rule grid below")262 print()263 print("table 3: same named rule across the mask family, class of every cell against the mask's k_min and k_2")264 print(f"{'rule':>13} {'mask':>10} {'side':>4} {'class':>6} {'k*':>3} {'k_min':>5} {'k_2':>5}")265 for rule, _, _ in RULES:266 for name, mask, side in masks:267 k, share, gain, density, churn, cls = results[(name, rule)]268 half, m, k_min, k_2, lobe, side, g, k_neg, neg = profiles[name]269 print(f"{rule:>13} {name:>10} {side:>4} {cls:>6} {k if cls in ('ring', 'coarse', 'flat') else '-':>3} {k_min:>5} {k_2:>5}")270 print()271 grid_rules = [(b_lo, b_hi, s_lo, s_hi) for b_lo, b_hi, s_lo, s_hi in itertools.product(GRID_B_LO, GRID_B_HI, GRID_S_LO, GRID_S_HI)]272 print(f"rule grid: birth lower {[fmt(x) for x in GRID_B_LO]}, birth upper {[fmt(x) for x in GRID_B_HI]}, survive lower {[fmt(x) for x in GRID_S_LO]}, survive upper {[fmt(x) for x in GRID_S_HI]}, {len(grid_rules)} rules, seeds {list(SEEDS)}, {len(grid_rules) * len(masks) * len(SEEDS)} runs")273 print()274 print("census: class counts per mask over the rule grid, then the ring stills against the mask")275 print(f"{'mask':>10} {'neg rings':>9} {'dead':>4} {'full':>4} {'active':>6} {'flat':>4} {'coarse':>6} {'ring':>4} {'k* range':>9} {'k*/k_2 range':>13} {'in lobe':>7} {'nearer k_2':>10} {'|k*-k_2|<=1':>11} {'|k*-k_min|<=1':>13} {'g(k*)<0':>7} {'|k*-k_neg|<=1':>13}")276 census = {}277 coarse = {}278 totals = Counter()279 for name, mask, side in masks:280 half, m, k_min, k_2, lobe, side, g, k_neg, neg = profiles[name]281 tally = Counter()282 rings = []283 coarse[name] = Counter()284 for seed in SEEDS:285 for birth_lo, birth_hi, survive_lo, survive_hi in grid_rules:286 birth, survive = (birth_lo, birth_hi), (survive_lo, survive_hi)287 grid, churn = evolve(half, m, birth, survive, seed)288 density = grid.mean()289 k, share, gain = spectrum(grid)290 cls = classify(density, churn, k, gain, side)291 tally[cls] += 1292 if cls == "ring":293 rings.append((rule_label(birth, survive), k, share, gain, density, seed))294 if cls == "coarse":295 coarse[name][k] += 1296 census[name] = rings297 totals.update(tally)298 neg_rings = int(np.sum(g[1:] < 0))299 ks = [r[1] for r in rings]300 in_lobe = sum(lobe[0] <= k <= lobe[1] for k in ks)301 nearer = sum(abs(k - k_2) < abs(k - k_min) for k in ks)302 at_2 = sum(abs(k - k_2) <= 1 for k in ks)303 at_min = sum(abs(k - k_min) <= 1 for k in ks)304 negative = sum(g[k] < 0 for k in ks)305 at_neg = sum(abs(k - k_neg) <= 1 for k in ks)306 k_range = f"{min(ks)}..{max(ks)}" if ks else "-"307 ratio = f"{min(ks) / k_2:.2f}..{max(ks) / k_2:.2f}" if ks else "-"308 print(f"{name:>10} {neg_rings:>9} {tally['dead']:>4} {tally['full']:>4} {tally['active']:>6} {tally['flat']:>4} {tally['coarse']:>6} {tally['ring']:>4} {k_range:>9} {ratio:>13} {in_lobe:>7} {nearer:>10} {at_2:>11} {at_min:>13} {negative:>7} {at_neg:>13}")309 print(f"{'total':>10} {'-':>9} {totals['dead']:>4} {totals['full']:>4} {totals['active']:>6} {totals['flat']:>4} {totals['coarse']:>6} {totals['ring']:>4}")310 print()311 print("ring stills per mask: k* histogram over the grid rules and seeds, the coarse stills' k* histogram, and the strongest ring still")312 for name, mask, side in masks:313 rings = census[name]314 hist = Counter(r[1] for r in rings)315 line = " ".join(f"{k}x{hist[k]}" for k in sorted(hist))316 coarse_line = " ".join(f"{k}x{coarse[name][k]}" for k in sorted(coarse[name]))317 best = max(rings, key=lambda r: r[3]) if rings else None318 best_line = f"{best[0]} seed {best[5]} k* {best[1]} N/k* {N / best[1]:.2f} share {best[2]:.3f} gain {best[3]:.1f} density {best[4]:.3f}" if best else "none"319 print(f"{name:>10}: ring {line if line else '-'}; coarse {coarse_line if coarse_line else '-'}; strongest {best_line}")320 print()321 print("every ring still: mask, rule, seed, k*, N/k*, g(k*)/m, share, gain, density")322 for name, mask, side in masks:323 half, m, k_min, k_2, lobe, side, g, k_neg, neg = profiles[name]324 for label, k, share, gain, density, seed in census[name]:325 print(f"{name:>10} {label} {seed:>5} {k:>3} {N / k:>6.2f} {g[k] / m:>8.4f} {share:>6.3f} {gain:>6.1f} {density:>6.3f}")326 print()327 laws = [328 ("k* of every ring still lies in the mask's half-height first side lobe", lambda k, p: p[4][0] <= k <= p[4][1]),329 ("k* of every ring still is nearer the mask's k_2 than its k_min", lambda k, p: abs(k - p[3]) < abs(k - p[2])),330 ("k* of every ring still sits within one ring of the mask's k_2", lambda k, p: abs(k - p[3]) <= 1),331 ("k* of every ring still sits within one ring of the mask's k_min", lambda k, p: abs(k - p[2]) <= 1),332 ("the ring-mean signed DFT of the mask is negative at k* of every ring still", lambda k, p: p[6][k] < 0),333 ("k* of every ring still sits within one ring of the mask's k_neg", lambda k, p: abs(k - p[7]) <= 1),334 ("k* of every ring still lies in the mask's negative band", lambda k, p: p[8][0] <= k <= p[8][1]),335 ]336 total = sum(len(v) for v in census.values())337 for text, test in laws:338 failures = [(name, r[0], r[1]) for name in census for r in census[name] if not test(r[1], profiles[name])]339 verdict = "Verified" if total and not failures else "Refuted"340 witness = f"; first witness {failures[0][0]} {failures[0][1]} k* {failures[0][2]}, {len(failures)} of {total} fail" if failures else f"; {total} ring stills"341 print(f"law: {text}: {verdict}{witness}")342 print()343 print("per-mask verdicts of the three sharpest laws over the ring stills, with the failing k* values")344 print(f"{'mask':>10} {'rings':>5} {'in negative band':>28} {'nearer k_2 than k_min':>28} {'within one ring of k_min':>28}")345 per_mask = [346 ("negative band", lambda k, p: p[8][0] <= k <= p[8][1]),347 ("nearer k_2", lambda k, p: abs(k - p[3]) < abs(k - p[2])),348 ("at k_min", lambda k, p: abs(k - p[2]) <= 1),349 ]350 verified_masks = []351 for name, mask, side in masks:352 cells = []353 for text, test in per_mask:354 bad = sorted(set(r[1] for r in census[name] if not test(r[1], profiles[name])))355 verdict = "Verified" if census[name] and not bad else "Refuted"356 if text == "negative band" and verdict == "Verified":357 verified_masks.append(name)358 cells.append(f"{verdict + (' k* ' + ' '.join(map(str, bad)) if bad else ''):>28}")359 print(f"{name:>10} {len(census[name]):>5} " + " ".join(cells))360 print(f"negative band Verified on {len(verified_masks)} masks ({', '.join(verified_masks)}) holding {sum(len(census[n]) for n in verified_masks)} ring stills")361 spreads = []362 for name in census:363 ks = sorted(set(r[1] for r in census[name]))364 if len(ks) >= 2:365 spreads.append((name, ks[0], ks[-1]))366 fixed = "Verified" if total and not spreads else "Refuted"367 witness = f"; first witness {spreads[0][0]} k* {spreads[0][1]}..{spreads[0][2]}" if spreads else ""368 print(f"law: every ring still on one mask shares one k* across the grid rules: {fixed}{witness}")369 print()370 print("scaling: carpet family per level, k_2 and the median N/k* over side of its ring stills")371 for name in ("code 7 L2", "code 7 L3", "code 7 L4"):372 half, m, k_min, k_2, lobe, side, g, k_neg, neg = profiles[name]373 ks = [r[1] for r in census[name]]374 med = f"{np.median([N / k / side for k in ks]):.3f}" if ks else "-"375 print(f"{name:>10} side {side:>3} k_2 {k_2:>3} N/k_2/side {N / k_2 / side:.3f} ring stills {len(ks):>3} median N/k*/side {med}")376 if not riesz_ok:377 sys.exit(1)378379380if __name__ == "__main__":381 main()