codebook.py
23.4 kB · python · 568 lines
1import math2import os3import zlib45import numpy as np6from scipy.signal import correlate78SIDE = 2439LEVEL = 510SEED_LIFE = 172911SEED_RANDOM = 812812SEED_NOISE = 49613SEED_CHECK = 656114LIFE_SIDE = 1615LIFE_GENS = 616LIFE_SEEDS = 817LIFE_DENSITY = 0.318CAP = 600000019CARPET2 = 720CARPET3 = 49521LETTERS = ((7, "carpet"), (14, "net"), (3, "htree"), (5, "vtree"), (9, "void"))22MAGIC_SIDES = ((3, 5), (5, 3), (3, 7), (7, 3))23RULES = (24 ((3,), (2, 3)),25 ((3, 6), (2, 3)),26 ((1, 3, 5, 7), (1, 3, 5, 7)),27 ((1, 3, 5, 7), (0, 2, 4, 6, 8)),28 ((0, 2, 4, 6, 8), (1, 3, 5, 7)),29 ((0, 2, 4, 6, 8), (0, 2, 4, 6, 8)),30 ((2,), ()),31 ((3, 4), (3, 4)),32)33BAYER = np.array([[0, 8, 2, 10], [12, 4, 14, 6], [3, 11, 1, 9], [15, 7, 13, 5]])34CLASSES = ("kronecker", "magic", "life")35HERE = os.path.dirname(os.path.abspath(__file__))36TREE = os.path.normpath(os.path.join(HERE, "..", "..", ".."))3738# TILES3940def tile(code, base, side):41 i = np.arange(side) % base42 bit = base * i[:, None] + i[None, :]43 return ((code >> bit) & 1).astype(np.uint8)4445def kron(a, b):46 return np.kron(a, b).astype(np.uint8)4748def power(a, level):49 out = a50 for _ in range(level - 1):51 out = kron(out, a)52 return out5354def square_orbit(code):55 g = tile(code, 3, 3)56 seen = set()57 for flip in (g, g[:, ::-1]):58 h = flip59 for _ in range(4):60 h = np.rot90(h)61 seen.add(int(sum(int(h[r, c]) << (3 * r + c) for r in range(3) for c in range(3))))62 return seen6364def base3_reps():65 seen = set()66 reps = []67 for code in range(1, 512):68 if code in seen:69 continue70 orbit = square_orbit(code)71 seen |= orbit72 reps.append(min(orbit))73 return sorted(reps)7475def magic_word(codes, sides):76 out = tile(codes[0], 2, sides[0])77 for code, side in zip(codes[1:], sides[1:]):78 out = kron(out, tile(code, 2, side))79 return out8081# LIFE8283def life_step(grid, birth, survive):84 n = np.zeros(grid.shape, dtype=np.int8)85 for dr in (-1, 0, 1):86 for dc in (-1, 0, 1):87 if dr or dc:88 n += np.roll(np.roll(grid, dr, axis=0), dc, axis=1)89 born = np.isin(n, birth) & (grid == 0)90 kept = np.isin(n, survive) & (grid == 1)91 return (born | kept).astype(np.uint8)9293def rule_name(birth, survive):94 return "rule birth [%s], survive [%s], wrap" % (" ".join(str(d) for d in birth), " ".join(str(d) for d in survive))9596def life_seeds():97 rng = np.random.default_rng(SEED_LIFE)98 out = [np.zeros((LIFE_SIDE, LIFE_SIDE), dtype=np.uint8)]99 out[0][LIFE_SIDE // 2, LIFE_SIDE // 2] = 1100 for _ in range(LIFE_SEEDS - 1):101 out.append((rng.random((LIFE_SIDE, LIFE_SIDE)) < LIFE_DENSITY).astype(np.uint8))102 return out103104# THE CATALOG105106def build_catalog(depth):107 reps = base3_reps()108 raw = []109 for code in range(16):110 raw.append(("kronecker", "bang dim 2, code %d, level 1" % code, tile(code, 2, 2)))111 for code in range(16):112 raw.append(("kronecker", "bang dim 2, code %d, level 2" % code, power(tile(code, 2, 2), 2)))113 for level in range(1, depth + 1):114 for code in reps:115 raw.append(("kronecker", "bang dim 2, base 3, code %d, level %d" % (code, level),116 power(tile(code, 3, 3), level)))117 for sides in MAGIC_SIDES:118 for ca, na in LETTERS:119 for cb, nb in LETTERS:120 raw.append(("magic", "%s(%d), %s(%d)" % (na, sides[0], nb, sides[1]),121 magic_word((ca, cb), sides)))122 for birth, survive in RULES:123 name = rule_name(birth, survive)124 for si, seed in enumerate(life_seeds()):125 grid = seed126 for gen in range(1, LIFE_GENS + 1):127 grid = life_step(grid, birth, survive)128 raw.append(("life", "%s s%d g%d" % (name, si, gen), grid.copy()))129 atoms = []130 index = {}131 collisions = 0132 for cls, label, arr in raw:133 if not arr.any():134 continue135 key = (arr.shape, arr.tobytes())136 if key in index:137 collisions += 1138 continue139 index[key] = len(atoms)140 atoms.append({"cls": cls, "label": label, "arr": arr, "h": int(arr.shape[0]),141 "w": int(arr.shape[1]), "ones": int(arr.sum())})142 return atoms, len(raw), collisions, len(reps)143144# CORPORA145146def corpus_tree():147 return power(tile(CARPET2, 2, 3), LEVEL)148149def corpus_text():150 need = (SIDE * SIDE + 7) // 8151 with open(os.path.join(TREE, "README.md"), "rb") as fh:152 data = fh.read()153 assert len(data) >= need154 bits = np.unpackbits(np.frombuffer(data[:need], dtype=np.uint8))155 return bits[: SIDE * SIDE].reshape(SIDE, SIDE).copy()156157def corpus_random():158 return (np.random.default_rng(SEED_RANDOM).random((SIDE, SIDE)) < 0.5).astype(np.uint8)159160def corpus_halftone():161 axis = (np.arange(SIDE) - (SIDE - 1) / 2) / ((SIDE - 1) / 2)162 grey = np.clip(1.0 - np.sqrt(axis[:, None] ** 2 + axis[None, :] ** 2), 0.0, 1.0)163 thresh = (BAYER[np.arange(SIDE)[:, None] % 4, np.arange(SIDE)[None, :] % 4] + 0.5) / 16.0164 return (grey > thresh).astype(np.uint8)165166# THE SCORE167168TABLES = {}169170def log2c_table(n):171 if n not in TABLES:172 if n <= 512:173 TABLES[n] = np.array([math.log2(math.comb(n, d)) for d in range(n + 1)])174 else:175 i = np.arange(1, n + 1, dtype=np.float64)176 TABLES[n] = np.concatenate(([0.0], np.cumsum(np.log2((n - i + 1) / i))))177 return TABLES[n]178179def log2c(n, k):180 if n <= 512:181 return math.log2(math.comb(n, k))182 if k == 0 or k == n:183 return 0.0184 return (math.lgamma(n + 1) - math.lgamma(k + 1) - math.lgamma(n - k + 1)) / math.log(2.0)185186def window_sums(x, h, w):187 ii = np.zeros((x.shape[0] + 1, x.shape[1] + 1), dtype=np.int64)188 ii[1:, 1:] = np.cumsum(np.cumsum(x.astype(np.int64), axis=0), axis=1)189 return ii[h:, w:] - ii[:-h, w:] - ii[h:, :-w] + ii[:-h, :-w]190191def mismatch(x, atom):192 corr = correlate(x.astype(np.float64), atom["arr"].astype(np.float64), mode="valid", method="fft")193 assert np.max(np.abs(corr - np.rint(corr))) < 1e-6194 return window_sums(x, atom["h"], atom["w"]) + atom["ones"] - 2 * np.rint(corr).astype(np.int64)195196def check_mismatch(x, atoms):197 rng = np.random.default_rng(SEED_CHECK)198 bad = 0199 for _ in range(200):200 atom = atoms[int(rng.integers(len(atoms)))]201 h, w = atom["h"], atom["w"]202 if h > x.shape[0] or w > x.shape[1]:203 continue204 r = int(rng.integers(x.shape[0] - h + 1))205 c = int(rng.integers(x.shape[1] - w + 1))206 direct = int((x[r:r + h, c:c + w] ^ atom["arr"]).sum())207 bad += direct != int(mismatch(x, atom)[r, c])208 return bad209210def floor_cells(log2a, positions):211 n = 1212 while n - (log2a + math.log2(positions) + math.log2(n + 1)) <= 0:213 n += 1214 return n215216# THE PURSUIT217218def candidates(x, atoms, tables):219 height, width = x.shape220 log2a = math.log2(len(atoms))221 sav_all, who_all, pos_all, cell_all = [], [], [], []222 live = 0223 for ai, atom in enumerate(atoms):224 h, w = atom["h"], atom["w"]225 if h > height or w > width:226 continue227 nw = h * w228 cols = width - w + 1229 positions = (height - h + 1) * cols230 base = log2a + math.log2(positions) + math.log2(nw + 1)231 if nw - base <= 0:232 continue233 live += 1234 sav = (nw - base) - tables[nw][mismatch(x, atom)]235 hit = np.nonzero(sav > 0)236 if hit[0].size == 0:237 continue238 sav_all.append(sav[hit].astype(np.float32))239 who_all.append(np.full(hit[0].size, ai, dtype=np.int32))240 pos_all.append((hit[0] * cols + hit[1]).astype(np.int32))241 cell_all.append(np.full(hit[0].size, nw, dtype=np.int32))242 if not sav_all:243 return np.zeros(0, np.float32), np.zeros(0, np.int32), np.zeros(0, np.int32), np.zeros(0, np.int32), live244 return (np.concatenate(sav_all), np.concatenate(who_all), np.concatenate(pos_all),245 np.concatenate(cell_all), live)246247def sweep(x, atoms, sav, who, pos, cell, key):248 height, width = x.shape249 score = sav if key == "flat" else sav / cell250 order = np.argsort(-score, kind="stable")251 capped = order.size > CAP252 if capped:253 order = order[:CAP]254 covered = np.zeros((height, width), dtype=bool)255 flat = bytearray(height * width)256 out = []257 for k in order:258 ai = int(who[k])259 atom = atoms[ai]260 h, w = atom["h"], atom["w"]261 cols = width - w + 1262 r, c = divmod(int(pos[k]), cols)263 if flat[r * width + c]:264 continue265 if covered[r:r + h, c:c + w].any():266 continue267 covered[r:r + h, c:c + w] = True268 for rr in range(r, r + h):269 flat[rr * width + c:rr * width + c + w] = b"\x01" * w270 out.append((ai, r, c))271 return out, capped272273def describe(x, atoms, placements):274 height, width = x.shape275 n = height * width276 log2a = math.log2(len(atoms))277 covered = np.zeros((height, width), dtype=bool)278 rebuilt = np.zeros((height, width), dtype=np.uint8)279 per_class = dict.fromkeys(CLASSES, 0.0)280 used = {c: {} for c in CLASSES}281 cost = 0.0282 for ai, r, c in placements:283 atom = atoms[ai]284 h, w = atom["h"], atom["w"]285 nw = h * w286 mask = x[r:r + h, c:c + w] ^ atom["arr"]287 d = int(mask.sum())288 this = (log2a + math.log2((height - h + 1) * (width - w + 1))289 + math.log2(nw + 1) + log2c(nw, d))290 cost += this291 per_class[atom["cls"]] += nw - this292 row = used[atom["cls"]].setdefault(ai, [0, 0.0])293 row[0] += 1294 row[1] += nw - this295 covered[r:r + h, c:c + w] = True296 rebuilt[r:r + h, c:c + w] = atom["arr"] ^ mask297 rebuilt[~covered] = x[~covered]298 assert np.array_equal(rebuilt, x)299 u = int((~covered).sum())300 u1 = int(x[~covered].sum())301 enum = math.log2(u + 1) + log2c(u, u1) if u else 0.0302 residual = 1.0 + min(float(u), enum)303 header = math.log2(n + 1)304 return {"bits": header + cost + residual, "header": header, "residual": residual,305 "uncovered": u, "credit": float(u) - residual, "per_class": per_class,306 "used": used, "placements": len(placements)}307308def deflate_bits(x):309 co = zlib.compressobj(9, zlib.DEFLATED, -15)310 return 8 * len(co.compress(np.packbits(x.ravel()).tobytes()) + co.flush())311312def run(x, atoms, tables, key="dense"):313 sav, who, pos, cell, live = candidates(x, atoms, tables)314 placements, capped = sweep(x, atoms, sav, who, pos, cell, key)315 full = describe(x, atoms, placements)316 bare = describe(x, atoms, [])317 assert abs(sum(full["per_class"].values()) - full["header"] + full["credit"]318 - (x.size - full["bits"])) < 1e-6319 won = full["bits"] < bare["bits"]320 rep = dict(full if won else bare)321 rep["bits"] += 1.0322 rep.update({"mode": "atoms" if won else "bare",323 "pursued": len(placements), "pursuit": full["bits"] + 1.0,324 "bare": bare["bits"] + 1.0,325 "live": live, "cands": int(sav.size), "capped": capped,326 "raw": x.size, "zlib": deflate_bits(x), "ones": int(x.sum())})327 return rep328329# THE REARRANGEMENT330331def rearrange(a, p, q):332 m, n = a.shape333 assert m % p == 0 and n % q == 0334 return a.reshape(p, m // p, q, n // q).transpose(0, 2, 1, 3).reshape(p * q, (m // p) * (n // q))335336def kron_spectrum(a, p, q):337 return np.linalg.svd(rearrange(a.astype(np.float64), p, q), compute_uv=False)338339def nearest_factors(a, p, q):340 u, s, vt = np.linalg.svd(rearrange(a.astype(np.float64), p, q), full_matrices=False)341 outer, inner = s[0] * u[:, 0], s[0] * vt[0]342 if outer[np.argmax(np.abs(outer))] < 0:343 outer, inner = -outer, -inner344 return outer.reshape(p, q), inner.reshape(a.shape[0] // p, a.shape[1] // q), s345346def block_means(a, p, q):347 m, n = a.shape348 return a.reshape(p, m // p, q, n // q).mean(axis=(1, 3))349350def threshold(m):351 return (m > 0.5 * m.max()).astype(np.uint8)352353def code_of(t):354 side = t.shape[0]355 return int(sum(int(t[r, c]) << (side * r + c) for r in range(side) for c in range(side)))356357def peel(a, sides):358 out, rest = [], a359 for side in sides:360 outer, inner, s = nearest_factors(rest, side, side)361 out.append((code_of(threshold(outer)), s[1] / s[0]))362 rest = threshold(inner)363 return out, rest364365# THE RUN366367def table(rows):368 print("%-12s %8s %8s %8s %9s %9s %9s" % ("corpus", "ones", "raw", "zlib", "codebook", "vs raw", "vs zlib"))369 for name, r in rows:370 bits = math.ceil(r["bits"])371 print("%-12s %8d %8d %8d %9d %9d %9d"372 % (name, r["ones"], r["raw"], r["zlib"], bits, r["raw"] - bits, r["zlib"] - bits))373374def classes(rows):375 print("%-12s %6s %5s %5s %5s %10s %10s %10s %9s"376 % ("corpus", "place", "kron", "magic", "life", "bits K", "bits M", "bits L", "residual"))377 for name, r in rows:378 u = r["used"]379 print("%-12s %6d %5d %5d %5d %10.0f %10.0f %10.0f %9.0f"380 % (name, r["placements"], len(u["kronecker"]), len(u["magic"]), len(u["life"]),381 r["per_class"]["kronecker"], r["per_class"]["magic"], r["per_class"]["life"],382 r["credit"]))383384def top_atoms(atoms, rep, k):385 rows = []386 for cls in CLASSES:387 for ai, (n, bits) in rep["used"][cls].items():388 rows.append((bits, n, cls, atoms[ai]))389 rows.sort(key=lambda t: -t[0])390 for bits, n, cls, atom in rows[:k]:391 print(" %-8s %-34s %3d x %-3d fill %5d %4d placements %8.0f bits"392 % (cls, atom["label"], atom["h"], atom["w"], atom["ones"], n, bits))393394def main():395 atoms, built, collisions, reps = build_catalog(2)396 assert reps == 101397 log2a = math.log2(len(atoms))398 sizes = sorted({a["h"] * a["w"] for a in atoms})399 tables = {s: log2c_table(s) for s in sizes}400 counts = {}401 for a in atoms:402 counts[(a["cls"], a["h"])] = counts.get((a["cls"], a["h"]), 0) + 1403404 print("== THE CATALOG ==")405 print("%d atoms built, %d duplicates dropped, %d empty frames dropped, %d kept"406 % (built, collisions, built - collisions - len(atoms), len(atoms)))407 for key in sorted(counts):408 print(" %-9s %3d x %-3d %4d atoms" % (key[0], key[1], key[1], counts[key]))409 print("the base-3 class is the %d nonempty square-group orbits of the 511 nonempty plane codes" % reps)410 print("naming one atom costs log2(%d) = %.4f bits, naming a position at most log2(%d) = %.4f"411 % (len(atoms), log2a, SIDE * SIDE, math.log2(SIDE * SIDE)))412 fl = floor_cells(log2a, SIDE * SIDE)413 dead = sum(1 for a in atoms if a["h"] * a["w"] < fl)414 print("an atom of n cells pays only when n > log2(%d) + log2(%d) + log2(n+1), so no atom below %d cells can ever be placed"415 % (len(atoms), SIDE * SIDE, fl))416 print("%d of the %d atoms sit below that floor: they can never pay and still charge every other atom their share of the name"417 % (dead, len(atoms)))418 print()419420 corpora = [421 ("tree render", "bang dim 2, code 7 at level %d, side %d" % (LEVEL, SIDE), corpus_tree()),422 ("text", "the first %d bytes of research/README.md, unpacked MSB first" % ((SIDE * SIDE + 7) // 8), corpus_text()),423 ("random", "uniform bits from numpy default_rng(%d)" % SEED_RANDOM, corpus_random()),424 ("halftone", "the radial gradient 1 - r ordered-dithered by the 4 x 4 Bayer matrix", corpus_halftone()),425 ]426 bad = sum(check_mismatch(x, atoms) for _, _, x in corpora)427 print("the correlation mismatch count agrees with a direct window comparison on 800 sampled placements, %d failures" % bad)428 assert bad == 0429 print()430431 rows = [(name, run(x, atoms, tables)) for name, _, x in corpora]432 print("== THE TABLE ==")433 table(rows)434 print()435 classes(rows)436 print()437 for (name, source, x), (_, r) in zip(corpora, rows):438 print("%-12s %s" % (name, source))439 print(" %d of %d cells covered, %d uncovered, reconstruction exact; %d atoms could be placed, %d candidate placements scored positive, cap hit %s"440 % (r["raw"] - r["uncovered"], r["raw"], r["uncovered"], r["live"], r["cands"], r["capped"]))441 print(" the pursuit takes %d placements and describes the corpus in %.0f bits, the placement-free description takes %.0f, the encoder emits the %s mode"442 % (r["pursued"], r["pursuit"], r["bare"], r["mode"]))443 print(" saving = %.0f atom bits - %.1f header + %.0f residual credit = %.0f"444 % (sum(r["per_class"].values()), r["header"], r["credit"], r["raw"] - r["bits"]))445 for cls in CLASSES:446 if r["used"][cls]:447 fills = [atoms[ai]["ones"] / (atoms[ai]["h"] * atoms[ai]["w"]) for ai in r["used"][cls]]448 print(" %-9s %2d distinct atoms placed, their densities running %.4f to %.4f"449 % (cls, len(fills), min(fills), max(fills)))450 top_atoms(atoms, r, 4)451 print()452453 print("== THE CONTROL ==")454 rnd = dict(rows)["random"]455 print("on uniform bits the codebook spends %.0f bits against a raw %d and a zlib %d, a saving of %.0f, with %d placements"456 % (rnd["bits"], rnd["raw"], rnd["zlib"], rnd["raw"] - rnd["bits"], rnd["pursued"]))457 print("it reads under zlib there only because deflate expands an incompressible stream by %d bits; neither compresses" % (rnd["zlib"] - rnd["raw"]))458 assert rnd["pursued"] == 0 and rnd["raw"] - rnd["bits"] <= 0459 txt = dict(rows)["text"]460 print("on text the codebook places nothing either, so its %.0f bit saving is the residual entropy code and not an atom"461 % (txt["raw"] - txt["bits"]))462 assert txt["pursued"] == 0463 print()464465 print("== THE GREEDY KEY ==")466 x = corpus_tree()467 sav, who, pos, cell, _ = candidates(x, atoms, tables)468 for key in ("dense", "flat"):469 pl, _ = sweep(x, atoms, sav, who, pos, cell, key)470 rep = describe(x, atoms, pl)471 print(" key %-5s: %4d placements, codebook %6.0f bits, saving %6.0f, K/M/L bits %6.0f %6.0f %6.0f"472 % (key, rep["placements"], rep["bits"] + 1.0, x.size - rep["bits"] - 1.0,473 rep["per_class"]["kronecker"], rep["per_class"]["magic"], rep["per_class"]["life"]))474 print("the table above uses the dense key, saving per cell; the flat key, saving per placement, spends the plane on large loose atoms")475 print()476477 print("== THE DEPTH LADDER ==")478 print("%-7s %7s %7s %9s %9s %9s %7s %10s %6s %8s %6s %8s"479 % ("depth", "atoms", "log2 A", "codebook", "vs raw", "vs zlib", "place", "kron bits",480 "r pl", "r save", "t pl", "t save"))481 for depth in (2, 3, 4, 5):482 cat, _, _, _ = build_catalog(depth)483 tab = {s: log2c_table(s) for s in sorted({a["h"] * a["w"] for a in cat})}484 rep = run(corpus_tree(), cat, tab)485 rnd2 = run(corpus_random(), cat, tab)486 txt2 = run(corpus_text(), cat, tab)487 assert rnd2["pursued"] == 0 and rnd2["raw"] - rnd2["bits"] <= 0488 bits = math.ceil(rep["bits"])489 print("%-7d %7d %7.4f %9d %9d %9d %7d %10.0f %6d %8.0f %6d %8.0f"490 % (depth, len(cat), math.log2(len(cat)), bits, rep["raw"] - bits, rep["zlib"] - bits,491 rep["placements"], rep["per_class"]["kronecker"],492 rnd2["pursued"], rnd2["raw"] - rnd2["bits"],493 txt2["pursued"], txt2["raw"] - txt2["bits"]))494 print("the random control takes zero placements at every depth, so a deeper catalog never manufactures a saving where there is none")495 print("text takes a handful from depth 3 on, each worth under two bits and each lowering the corpus total, so the encoder drops them and emits the placement-free mode")496 print("the level-%d tile of the tree corpus is itself an atom at depth %d, so the last row is recognition and not compression" % (LEVEL, LEVEL))497 print()498499 print("== THE KRONECKER SPECTRUM ==")500 carpet = corpus_tree()501 for j in (1, 2, 3, 4):502 side = 3 ** j503 outer, inner, s = nearest_factors(carpet, side, side)504 assert np.array_equal(threshold(outer), power(tile(CARPET2, 2, 3), j))505 assert np.array_equal(threshold(inner), power(tile(CARPET2, 2, 3), LEVEL - j))506 print(" split %3d x %-3d sigma_1 = %.4f, sigma_2/sigma_1 = %.3e; the rank-one factors are carpet levels %d and %d exactly"507 % (side, side, s[0], s[1] / s[0], j, LEVEL - j))508 got = code_of(threshold(nearest_factors(carpet, 3, 3)[0]))509 assert got == CARPET3510 print(" at the 3 x 3 split the recovered base-3 code is %d, the carpet's own level-1 code" % got)511 print()512513 word = (7, 9, 14, 3, 5)514 names = dict(LETTERS)515 read, rest = peel(magic_word(word, (3,) * 5), (3,) * 4)516 got = [c for c, _ in read] + [code_of(rest)]517 want = [code_of(tile(c, 2, 3)) for c in word]518 assert got == want519 print("== THE PEEL ==")520 print(" the magic word %s renders at side %d and the rearranged SVD peels it letter by letter"521 % (", ".join("%s(3)" % names[c] for c in word), SIDE))522 print(" recovered base-3 codes %s, the word's own %s, ratios sigma_2/sigma_1 %s"523 % (got, want, ", ".join("%.1e" % r for _, r in read)))524 print()525526 print("== THE NOISE DIAL ==")527 rng = np.random.default_rng(SEED_NOISE)528 trials = 40529 svd, mean = [], []530 for k in range(0, 51):531 p = k / 100.0532 a = b = 0533 for _ in range(trials):534 noisy = (carpet ^ (rng.random(carpet.shape) < p)).astype(np.uint8)535 a += code_of(threshold(nearest_factors(noisy, 3, 3)[0])) == CARPET3536 b += code_of(threshold(block_means(noisy, 3, 3))) == CARPET3537 svd.append((p, a))538 mean.append((p, b))539 print(" flipping a fraction p of the %d bits, %d seeds per p on the grid 0.00 to 0.50 in steps of 0.01" % (SIDE * SIDE, trials))540 for name, curve in (("rank one", svd), ("block mean", mean)):541 clean = max(p for p, ok in curve if all(o == trials for q, o in curve if q <= p))542 first = min((p for p, ok in curve if ok < trials), default=None)543 half = min((p for p, ok in curve if ok * 2 < trials), default=None)544 none = min((p for p, ok in curve if ok == 0), default=None)545 print(" %-11s every seed to p = %.2f, first failure p = %.2f, half the seeds by p = %.2f, none from p = %.2f"546 % (name, clean, first, half, none))547 print(" %s" % " ".join("%.2f:%d" % (p, ok) for p, ok in curve if 0.20 <= p <= 0.35))548 fill = (8.0 / 9.0) ** (LEVEL - 1)549 star = fill / (1.0 + 2.0 * fill)550 half_mean = min(p for p, ok in mean if ok * 2 < trials)551 print(" the block-mean detector's expected margin vanishes at f/(1 + 2f) = %.6f with f = (8/9)^%d = %.6f the density of a filled corner block,"552 % (star, LEVEL - 1, fill))553 print(" a closed form read off the fill law and not off this sweep; the sweep keeps more than half its seeds at p = %.2f and loses them by p = %.2f,"554 % (half_mean - 0.01, half_mean))555 print(" one grid step above the closed form because the threshold divides the largest of eight noisy filled corners rather than their mean")556 assert abs(half_mean - star) <= 0.02557 print(" the rank-one factor of the rearrangement holds %.2f further and falls off a cliff two grid steps wide"558 % (min(p for p, ok in svd if ok * 2 < trials) - half_mean))559 for p in (0.00, 0.10, 0.20, 0.30, 0.40, 0.50):560 noisy = (carpet ^ (rng.random(carpet.shape) < p)).astype(np.uint8)561 s = kron_spectrum(noisy, 3, 3)562 print(" p = %.2f: sigma_2/sigma_1 = %.4f, recovered code %d against %d"563 % (p, s[1] / s[0], code_of(threshold(nearest_factors(noisy, 3, 3)[0])), CARPET3))564 print(" the spectrum itself degrades smoothly and reads nothing, while the thresholded rank-one factor stays exact far past it")565 print()566 print("every assertion passed")567568main()