measures.py
19.0 kB · python · 553 lines
1import csv2import os3import random4import time5from fractions import Fraction6from itertools import permutations7from math import comb89import numpy as np1011HERE = os.path.dirname(os.path.abspath(__file__))12DIMS = (3, 4)13MEASURES = ("s", "bs", "C", "dt", "deg", "dnf", "cnf")14KEYS = (15 ("genus",),16 ("gf2deg",),17 ("pop",),18 ("gf2deg", "pop"),19 ("fingerprint",),20 ("genus", "fingerprint"),21 ("genus", "gf2deg", "pop", "fingerprint"),22)23LETTERS = "xyzw"24SEED = 2026072525TRIALS = 5000026SAMPLE = 2000027LEVELS = (1, 2, 3, 4, 6, 8, 12, 16)28NAMED = (29 "bang dim 4, code 27",30 "bang dim 4, code 281",31 "bang dim 4, code 855",32 "bang dim 4, code 1911",33 "bang dim 4, code 7128",34)35POP16 = np.array([i.bit_count() for i in range(1 << 16)], dtype=np.uint8)3637def corner_maps(d):38 n = 1 << d39 maps = []40 for pi in permutations(range(d)):41 for flip in range(n):42 m = [0] * n43 for c in range(n):44 image = 045 for a in range(d):46 b = ((c >> (d - 1 - pi[a])) & 1) ^ ((flip >> (d - 1 - a)) & 1)47 image |= b << (d - 1 - a)48 m[c] = image49 maps.append([1 << target for target in m])50 return maps5152def catalog(d):53 n = 1 << d54 maps = corner_maps(d)55 seen = bytearray(1 << n)56 out = []57 for code in range(1 << n):58 if seen[code]:59 continue60 cells = [i for i in range(n) if (code >> i) & 1]61 orbit = set()62 for m in maps:63 moved = 064 for i in cells:65 moved |= m[i]66 orbit.add(moved)67 for member in orbit:68 seen[member] = 169 out.append((code, sorted(orbit)))70 return out7172def bit(f, x):73 return (f >> x) & 17475def differing(f, x, n):76 return f ^ ((1 << n) - 1) if bit(f, x) else f7778def cubes(d):79 n = 1 << d80 out = {}81 for fixed in range(n):82 for vals in range(n):83 if vals & fixed == vals:84 out[(fixed, vals)] = sum(1 << y for y in range(n) if y & fixed == vals)85 return out8687def constant(f, mask):88 return f & mask in (0, mask)8990def sensitivity(f, d):91 n = 1 << d92 best = 093 for x in range(n):94 m = differing(f, x, n)95 best = max(best, sum(bit(m, x ^ (1 << i)) for i in range(d)))96 return best9798def pack(blocks):99 memo = {}100101 def go(used):102 if used in memo:103 return memo[used]104 top = 0105 for b in blocks:106 if b & used == 0:107 top = max(top, 1 + go(used | b))108 memo[used] = top109 return top110111 return go(0)112113def block_sensitivity(f, d):114 n = 1 << d115 best = 0116 for x in range(n):117 m = differing(f, x, n)118 best = max(best, pack([x ^ y for y in range(n) if bit(m, y)]))119 return best120121def certificate(f, d, cube):122 n = 1 << d123 order = sorted(range(n), key=int.bit_count)124 best = 0125 for x in range(n):126 for fixed in order:127 if constant(f, cube[(fixed, x & fixed)]):128 best = max(best, fixed.bit_count())129 break130 return best131132def dt_depth(f, d, cube):133 memo = {}134135 def go(fixed, vals):136 key = (fixed, vals)137 if key in memo:138 return memo[key]139 if constant(f, cube[key]):140 out = 0141 else:142 out = d143 for i in range(d):144 if not (fixed >> i) & 1:145 low = go(fixed | 1 << i, vals)146 high = go(fixed | 1 << i, vals | 1 << i)147 out = min(out, 1 + max(low, high))148 memo[key] = out149 return out150151 return go(0, 0)152153def mobius(f, d, gf2):154 n = 1 << d155 coef = [bit(f, x) for x in range(n)]156 for i in range(d):157 step = 1 << i158 for s in range(n):159 if s & step:160 coef[s] = coef[s] ^ coef[s ^ step] if gf2 else coef[s] - coef[s ^ step]161 return coef162163def top_weight(coef):164 return max((s.bit_count() for s, c in enumerate(coef) if c), default=-1)165166def anf(f, d):167 words = []168 for s, c in enumerate(mobius(f, d, True)):169 if c:170 word = "".join(LETTERS[a] for a in range(d) if (s >> (d - 1 - a)) & 1)171 words.append((len(word), word or "1"))172 return " + ".join(word for _, word in sorted(words)) or "0"173174def min_cover(primes, need):175 best = [len(primes)]176177 def go(left, count):178 if left == 0:179 best[0] = min(best[0], count)180 return181 if count + 1 >= best[0]:182 return183 cells = [y for y in range(left.bit_length()) if bit(left, y)]184 y = min(cells, key=lambda c: sum(1 for p in primes if bit(p, c)))185 for p in primes:186 if bit(p, y):187 go(left & ~p, count + 1)188189 go(need, 0)190 return best[0]191192def cover_size(f, d, cube, target):193 n = 1 << d194 need = f if target else f ^ ((1 << n) - 1)195 if need == 0:196 return 0197 good = {key: mask for key, mask in cube.items() if mask & need == mask}198 primes = []199 for (fixed, vals), mask in good.items():200 freed = ((fixed & ~(1 << i), vals & ~(1 << i)) for i in range(d) if (fixed >> i) & 1)201 if all(key not in good for key in freed):202 primes.append(mask)203 return min_cover(primes, need)204205def all_measures(f, d, cube):206 return {207 "s": sensitivity(f, d),208 "bs": block_sensitivity(f, d),209 "C": certificate(f, d, cube),210 "dt": dt_depth(f, d, cube),211 "deg": top_weight(mobius(f, d, False)),212 "dnf": cover_size(f, d, cube, 1),213 "cnf": cover_size(f, d, cube, 0),214 }215216def fingerprint(code, d):217 coef = [0] * (d + 1)218 for i in range(1 << d):219 if bit(code, i):220 ones = i.bit_count()221 for j in range(ones + 1):222 coef[d - ones + j] += comb(ones, j) * (-1) ** (ones - j)223 return tuple(reversed(coef))224225def poly_text(coef):226 d = len(coef) - 1227 out = ""228 for j, c in enumerate(coef):229 if c == 0:230 continue231 power = d - j232 size = "" if abs(c) == 1 and power else str(abs(c))233 var = "" if power == 0 else "k" if power == 1 else f"k^{power}"234 sign = "-" if c < 0 else "+"235 out += f" {sign} {size}{var}" if out else f"{'-' if c < 0 else ''}{size}{var}"236 return out or "0"237238def render_index(d, side):239 grid = np.indices((side,) * d).reshape(d, -1) & 1240 weights = np.array([1 << (d - 1 - a) for a in range(d)])241 return (grid * weights[:, None]).sum(axis=0)242243def fitted(code, d, renders):244 m = d + 1245 table = np.array([bit(code, i) for i in range(1 << d)], dtype=np.int64)246 rows = []247 for k in range(1, m + 1):248 fill = int(table[renders[k]].sum())249 rows.append([Fraction(k ** (d - j)) for j in range(m)] + [Fraction(fill)])250 for col in range(m):251 pivot = next(r for r in range(col, m) if rows[r][col] != 0)252 rows[col], rows[pivot] = rows[pivot], rows[col]253 lead = rows[col][col]254 rows[col] = [v / lead for v in rows[col]]255 for r in range(m):256 if r != col and rows[r][col] != 0:257 factor = rows[r][col]258 rows[r] = [a - factor * b for a, b in zip(rows[r], rows[col])]259 return tuple(int(row[m]) for row in rows)260261def is_levelset(code, d):262 value = [-1] * (d + 1)263 for i in range(1 << d):264 w = i.bit_count()265 if value[w] not in (-1, bit(code, i)):266 return False267 value[w] = bit(code, i)268 return True269270def is_pin(code, d):271 cells = [i for i in range(1 << d) if bit(code, i)]272 if not cells:273 return False274 fixed = sum(1 for i in range(d) if len({bit(c, i) for c in cells}) == 1)275 return len(cells) == 1 << (d - fixed)276277def genus(orbit, d):278 if any(is_levelset(c, d) for c in orbit):279 return "iso"280 if any(is_pin(c, d) for c in orbit):281 return "axis"282 return "comp"283284def rows_of(d):285 cube = cubes(d)286 out = []287 for code, orbit in catalog(d):288 row = {289 "name": f"bang dim {d}, code {code}",290 "code": code,291 "orbit": len(orbit),292 "genus": genus(orbit, d),293 "gf2deg": top_weight(mobius(code, d, True)),294 "pop": code.bit_count(),295 "fingerprint": fingerprint(code, d),296 }297 row.update(all_measures(code, d, cube))298 out.append(row)299 return out300301def csv_cells(row):302 fp = " ".join(str(c) for c in row["fingerprint"])303 head = [row["name"], str(row["code"]), str(row["orbit"]), row["genus"], str(row["gf2deg"]), str(row["pop"]), fp]304 return head + [str(row[m]) for m in MEASURES]305306def csv_diff(rows, d):307 with open(os.path.join(HERE, f"measures_d{d}.csv"), newline="") as handle:308 stored = list(csv.reader(handle))309 compared = 0310 bad = 0311 for row, kept in zip(rows, stored[1:]):312 for mine, theirs in zip(csv_cells(row)[1:], kept[1:]):313 compared += 1314 bad += mine != theirs315 bad += abs(len(rows) - len(stored) + 1) * 13316 return compared, bad317318def key_of(row, fields):319 return tuple(row[field] for field in fields)320321def groups(rows, fields):322 out = {}323 for row in rows:324 out.setdefault(key_of(row, fields), []).append(row)325 return out326327def splits(group, measure):328 return len({row[measure] for row in group}) > 1329330def determination(rows):331 out = {}332 for measure in MEASURES:333 out[measure] = "none"334 for fields in KEYS:335 if not any(splits(group, measure) for group in groups(rows, fields).values()):336 out[measure] = " + ".join(fields)337 break338 return out339340def pin_code(d, r):341 return sum(1 << i for i in range(1 << d) if i >> (d - r) == 0)342343def popcount(a):344 return POP16[a & 0xFFFF].astype(np.int8) + POP16[a >> 16].astype(np.int8)345346def families(d):347 out = []348349 def go(remaining, chosen):350 if remaining == 0:351 out.append(tuple(chosen))352 return353 low = remaining & -remaining354 rest = remaining ^ low355 go(rest, chosen)356 sub = rest357 while True:358 go(rest ^ sub, chosen + [low | sub])359 if sub == 0:360 break361 sub = (sub - 1) & rest362363 go((1 << d) - 1, [])364 return out365366def profile(d, codes):367 n = 1 << d368 full = np.uint32((1 << n) - 1)369 cube = cubes(d)370 by_size = {}371 for family in families(d):372 by_size.setdefault(len(family), []).append(family)373 by_codim = {}374 for fixed in range(n):375 by_codim.setdefault(fixed.bit_count(), []).append(fixed)376 s = np.zeros(len(codes), np.int8)377 bs = np.zeros(len(codes), np.int8)378 C = np.zeros(len(codes), np.int8)379 for x in range(n):380 m = codes ^ (((codes >> x) & 1) * full)381 nb = np.uint32(sum(1 << (x ^ (1 << i)) for i in range(d)))382 s = np.maximum(s, popcount(m & nb))383 bx = np.zeros(len(codes), np.int8)384 for k in range(1, d + 1):385 hit = np.zeros(len(codes), bool)386 for family in by_size[k]:387 fm = np.uint32(sum(1 << (x ^ b) for b in family))388 hit |= (m & fm) == fm389 if not hit.any():390 break391 bx[hit] = k392 bs = np.maximum(bs, bx)393 cx = np.full(len(codes), d, np.int8)394 for c in range(d - 1, -1, -1):395 hit = np.zeros(len(codes), bool)396 for fixed in by_codim[c]:397 mask = np.uint32(cube[(fixed, x & fixed)])398 v = codes & mask399 hit |= (v == 0) | (v == mask)400 if not hit.any():401 break402 cx[hit] = c403 C = np.maximum(C, cx)404 return s, bs, C405406def histogram(s, bs, C):407 keys, counts = np.unique(np.stack([s, bs, C], axis=1), axis=0, return_counts=True)408 return [(tuple(int(v) for v in key), int(count)) for key, count in zip(keys, counts)]409410def print_profile(label, s, bs, C):411 print(f"{label}: designs {len(s)}, C != bs {int((C != bs).sum())}, bs - s > 1 {int((bs - s > 1).sum())}")412 for (a, b, c), count in histogram(s, bs, C):413 print(f" (s, bs, C) = ({a}, {b}, {c}): {count}")414415def table_line(row, d):416 cells = [row["name"], row["genus"], str(row["gf2deg"]), str(row["pop"])]417 if d == 4:418 cells.append(poly_text(row["fingerprint"]))419 return " | ".join(cells + [str(row[m]) for m in MEASURES])420421def counts(rows, field):422 out = {}423 for row in rows:424 out[row[field]] = out.get(row[field], 0) + 1425 return ", ".join(f"{key}: {out[key]}" for key in sorted(out))426427def influence(d):428 n = 1 << d429 codes = np.arange(1 << n, dtype=np.uint32)430 full = np.uint32((1 << n) - 1)431 total = np.zeros(len(codes), np.int64)432 top = np.zeros(len(codes), np.int8)433 for x in range(n):434 m = codes ^ (((codes >> x) & 1) * full)435 nb = np.uint32(sum(1 << (x ^ (1 << i)) for i in range(d)))436 sx = popcount(m & nb)437 total += sx438 top = np.maximum(top, sx)439 count = len(codes)440 mean = Fraction(int(total.sum()), count * n)441 var = Fraction(int((total * total).sum()), count * n * n) - mean * mean442 full_share = Fraction(int((top == d).sum()), count)443 one_edge = int((total == 2).sum())444 edges = sorted(set((total // 2).tolist()))445 return mean, var, full_share, one_edge, edges446447def sample_bits():448 rng = random.Random(SEED)449 codes = []450 for _ in range(TRIALS):451 f = 0452 for i in range(32):453 f |= rng.getrandbits(1) << i454 codes.append(f)455 return np.array(codes, dtype=np.uint32)456457def sample_words():458 rng = random.Random(SEED)459 uniform = [rng.getrandbits(32) for _ in range(SAMPLE)]460 thinned = []461 for _ in range(SAMPLE):462 p = LEVELS[rng.randint(0, len(LEVELS) - 1)]463 code = 0464 for i in range(32):465 if rng.randint(0, 31) < p:466 code |= 1 << i467 thinned.append(code)468 return np.array(uniform, dtype=np.uint32), np.array(thinned, dtype=np.uint32)469470def main():471 start = time.time()472 store = {}473 for d in DIMS:474 rows = rows_of(d)475 store[d] = rows476 print(f"D={d} classes {len(rows)}, designs {sum(row['orbit'] for row in rows)}")477 compared, bad = csv_diff(rows, d)478 print(f"D={d} csv cells compared {compared}, mismatches {bad}")479 renders = {k: render_index(d, 2 * k - 1) for k in range(1, d + 2)}480 bad = sum(fitted(row["code"], d, renders) != row["fingerprint"] for row in rows)481 print(f"D={d} fill polynomial fitted from rendered grids against closed form, mismatches {bad} of {len(rows)}")482 print("D=3 table: name | genus | gf2deg | pop | " + " | ".join(MEASURES))483 for row in store[3]:484 print(" " + table_line(row, 3))485 print(f"D=3 genus split {counts(store[3], 'genus')}")486 print(f"D=3 gf2deg histogram {counts(store[3], 'gf2deg')}")487 deeper = [f"{row['name']} (bs {row['bs']}, dt {row['dt']})" for row in store[3] if row["dt"] > row["bs"]]488 print(f"D=3 dt > bs classes {len(deeper)}: {deeper}")489 for d in DIMS:490 cube = cubes(d)491 for r in range(d + 1):492 got = all_measures(pin_code(d, r), d, cube)493 got["gf2deg"] = top_weight(mobius(pin_code(d, r), d, True))494 six = {got[m] for m in ("s", "bs", "C", "dt", "deg", "gf2deg")}495 print(f"D={d} pin of {r} axes: six measures {sorted(six)}, dnf {got['dnf']}, cnf {got['cnf']}")496 d3 = store[3]497 d4 = store[4]498 by_name = {row["name"]: row for row in d4}499 fp3 = groups(d3, ("fingerprint",))500 shared3 = [group for group in fp3.values() if len(group) > 1]501 print(f"D=3 distinct fill polynomials {len(fp3)} of {len(d3)}, shared by two or more {len(shared3)}")502 for group in shared3:503 names = ", ".join(f"{row['name']} ({row['genus']}, gf2deg {row['gf2deg']})" for row in group)504 print(f"D=3 collision {poly_text(group[0]['fingerprint'])}: {names}")505 report = determination(d3)506 print("D=3 coarsest key determining each measure: " + ", ".join(f"{m} {report[m]}" for m in MEASURES))507 fp4 = groups(d4, ("fingerprint",))508 shared4 = sum(1 for group in fp4.values() if len(group) > 1)509 print(f"D=4 distinct fill polynomials {len(fp4)} of {len(d4)}, shared by two or more {shared4}")510 report = determination(d4)511 print("D=4 coarsest key determining each measure: " + ", ".join(f"{m} {report[m]}" for m in MEASURES))512 big = [group for group in groups(d4, KEYS[-1]).values() if len(group) > 1]513 parted = sum(1 for group in big if any(splits(group, m) for m in MEASURES))514 total = sum(1 for group in big for m in MEASURES if splits(group, m))515 print(f"D=4 full-key groups of two or more {len(big)}, splitting some measure {parted}, measure-splits {total}")516 pairs = [group for group in big if len(group) == 2 and splits(group, "bs")]517 pairs.sort(key=lambda group: min(row["code"] for row in group))518 first = pairs[0]519 names = " and ".join(row["name"] for row in first)520 split_names = [m for m in MEASURES if splits(first, m)]521 agree = [m for m in MEASURES if not splits(first, m)]522 print(f"D=4 size-two full-key groups splitting bs {len(pairs)}, smallest code {first[0]['code']}: {names}")523 print(f" shared key genus {first[0]['genus']}, gf2deg {first[0]['gf2deg']}, pop {first[0]['pop']}, fill {poly_text(first[0]['fingerprint'])}")524 print(f" split {len(split_names)} of 7 measures {split_names}, agree {agree}")525 print("D=4 witness table: name | genus | gf2deg | pop | fill | " + " | ".join(MEASURES))526 for row in first:527 print(" " + table_line(row, 4))528 for d, rows in store.items():529 gaps = [row["name"] for row in rows if row["C"] != row["bs"]]530 print(f"D={d} C != bs classes: {gaps}")531 for d, rows in store.items():532 seps = [f"{row['name']} (s {row['s']}, bs {row['bs']}, orbit {row['orbit']})" for row in rows if row["s"] != row["bs"]]533 print(f"D={d} s != bs classes: {seps}")534 both = d3 + d4535 print(f"deg <= s^2 violations across both catalogs: {sum(1 for row in both if row['deg'] > row['s'] ** 2)} of {len(both)}")536 tight = [f"{row['name']} (deg {row['deg']}, s {row['s']}, orbit {row['orbit']})" for row in d4 if row["s"] >= 2 and row["deg"] == row["s"] ** 2]537 print(f"D=4 deg = s^2 with s >= 2: {tight}")538 for name in NAMED:539 print(f"{name} anf {anf(by_name[name]['code'], 4)}")540 for d in range(1, 5):541 mean, var, full_share, one_edge, edges = influence(d)542 print(f"D={d} influence mean {mean}, variance {var}, P[s(f) = D] {full_share}, designs with one bichromatic edge {one_edge}")543 print(f"D={d} realised bichromatic edge counts {edges}")544 codes = np.arange(1 << 16, dtype=np.uint32)545 print_profile("D=4 exhaustive (s, bs, C) profile by the subcube route", *profile(4, codes))546 print_profile(f"D=5 spot check, {TRIALS} uniform designs, seed {SEED}, one bit per draw", *profile(5, sample_bits()))547 uniform, thinned = sample_words()548 print_profile(f"D=5 spot check, {SAMPLE} uniform designs, seed {SEED}, one word per draw", *profile(5, uniform))549 print_profile(f"D=5 spot check, {SAMPLE} thinned designs, densities {LEVELS} of 32", *profile(5, thinned))550 print(f"domain: every design at D = 1..4, {TRIALS} + {SAMPLE} + {SAMPLE} sampled designs at D = 5, {time.time() - start:.1f} s")551552if __name__ == "__main__":553 main()