#!/usr/bin/env python3 """Reproduce the finite-group calculations and error constants. Run from any directory: python3 code/reproduce.py Core state sums, projections, and blocking checks use exact rational arithmetic. NumPy and Matplotlib are used only for spectra, graph constants, and plots. All output paths are relative to this project. """ from __future__ import annotations import csv import itertools import json import math import os import platform import tempfile from fractions import Fraction as F from pathlib import Path import numpy as np os.environ.setdefault("MPLCONFIGDIR", str(Path(tempfile.gettempdir()) / "spin_foam_matplotlib")) import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt ROOT = Path(__file__).resolve().parents[2] DATA, FIGURES = ROOT / "anc" / "data", ROOT / "figures" DATA.mkdir(exist_ok=True) FIGURES.mkdir(exist_ok=True) def mm(a, b): return [[sum((x * y for x, y in zip(row, col)), F(0)) for col in zip(*b)] for row in a] def transpose(a): return [list(x) for x in zip(*a)] def scale(a, t): return [[t * x for x in row] for row in a] def subtract(a, b): return [[x - y for x, y in zip(ar, br)] for ar, br in zip(a, b)] def identity(n): return [[F(i == j) for j in range(n)] for i in range(n)] def power(a, n): result = identity(len(a)) for _ in range(n): result = mm(result, a) return result def diag(values): return [[F(x) if i == j else F(0) for j in range(len(values))] for i, x in enumerate(values)] def asfloat(a): return np.array([[float(x) for x in row] for row in a]) def save_csv(name, fieldnames, rows): with (DATA / name).open("w", newline="", encoding="utf-8") as stream: writer = csv.DictWriter(stream, fieldnames=fieldnames) writer.writeheader() writer.writerows(rows) def compositions(n): if n == 0: yield () else: for first in range(1, n + 1): for rest in compositions(n - first): yield (first,) + rest def finite_group_checks(): """Directly enumerate nonabelian as well as abelian BF state sums.""" cyclic = lambda n: (list(range(n)), lambda a, b: (a + b) % n, lambda a: (-a) % n, 0) klein = (list(itertools.product(range(2), repeat=2)), lambda a, b: tuple(x ^ y for x, y in zip(a, b)), lambda a: a, (0, 0)) perms = list(itertools.permutations(range(3))) smul = lambda a, b: tuple(a[b[i]] for i in range(3)) sinv = lambda a: tuple(a.index(i) for i in range(3)) groups = [("Z2", cyclic(2)), ("Z3", cyclic(3)), ("Z2xZ2", klein), ("S3", (perms, smul, sinv, (0, 1, 2)))] rows = [] count = 0 for name, (group, mul, inv, unit) in groups: q = len(group) delta = lambda x: q if x == unit else 0 # A polygon with n boundary edges has integral one. for n in range(1, 5): total = 0 for edges in itertools.product(group, repeat=n): hol = unit for edge in edges: hol = mul(hol, edge) total += delta(hol) value = F(total, q ** n) assert value == 1 rows.append(dict(group=name, test=f"disk_{n}_edges", value=value)) count += 1 # Two BF faces glued along one edge: Delta * Delta = Delta. for a, b in itertools.product(group, repeat=2): value = F(sum(delta(mul(a, h)) * delta(mul(inv(h), b)) for h in group), q) assert value == delta(mul(a, b)) count += 1 # Three triangles meeting at one internal vertex; fixed boundary a,b,c. # Every face receives q*1_flat and every internal edge receives 1/q. for a, b, c in itertools.product(group, repeat=3): flat_assignments = 0 for r0, r1, r2 in itertools.product(group, repeat=3): faces = (mul(mul(r0, a), inv(r1)), mul(mul(r1, b), inv(r2)), mul(mul(r2, c), inv(r0))) flat_assignments += all(h == unit for h in faces) kernel = F(flat_assignments * q ** 3, q ** 3) assert kernel == delta(mul(mul(a, b), c)) count += 1 sphere = F(sum(delta(g) ** 2 for g in group), q) assert sphere == q rows.append(dict(group=name, test="sphere_from_two_disks", value=sphere)) rows.append(dict(group=name, test="groupoid_sphere", value=F(1, q))) count += 1 save_csv("bf_normalization.csv", ["group", "test", "value"], rows) return count def transfer_checks(alpha, beta, tag): """Compute from the group convolution matrix before using Fourier formulas.""" group = list(itertools.product(range(2), repeat=2)) chi1 = [F((-1) ** x) for x, y in group] chi2 = [F((-1) ** y) for x, y in group] characters = [[F(1), chi1[i], chi2[i], chi1[i] * chi2[i]] for i in range(4)] weights = [1 + alpha * chi1[i] + beta * chi2[i] for i in range(4)] assert all(w > 0 for w in weights) and sum(weights) == 4 index = {g: i for i, g in enumerate(group)} matrix = [[weights[index[(x[0] ^ y[0], x[1] ^ y[1])]] / 4 for y in group] for x in group] fourier = scale(mm(mm(transpose(characters), matrix), characters), F(1, 4)) assert fourier == diag([1, alpha, beta, 0]) projection = [[F(1, 4) + (chi1[i] + chi2[i]) * (chi1[j] + chi2[j]) / 8 for j in range(4)] for i in range(4)] assert mm(projection, projection) == projection assert transpose(projection) == projection retained = mm(mm(projection, matrix), projection) complement = subtract(identity(4), projection) gap = subtract(mm(mm(projection, power(matrix, 2)), projection), power(retained, 2)) leakage = mm(mm(complement, matrix), projection) assert gap == mm(transpose(leakage), leakage) delta = (alpha - beta) ** 2 / 4 gap_eigenvalues = np.linalg.eigvalsh(asfloat(gap)) assert np.allclose(gap_eigenvalues, [0, 0, 0, float(delta)], atol=1e-14) if beta == -alpha: assert power(retained, 2) == retained # Check the analytic compression formula by applying its group-space projector. powers = {n: mm(mm(projection, power(matrix, n)), projection) for n in range(1, 13)} mode_projection = [[(chi1[i] + chi2[i]) * (chi1[j] + chi2[j]) / 8 for j in range(4)] for i in range(4)] vacuum_projection = [[F(1, 4)] * 4 for _ in range(4)] compressed_coefficient = lambda n: (alpha ** n + beta ** n) / 2 for n in powers: expected = [[vacuum_projection[i][j] + compressed_coefficient(n) * mode_projection[i][j] for j in range(4)] for i in range(4)] assert powers[n] == expected residual_rows = [] for m, n in itertools.product(range(1, 7), repeat=2): residual = subtract(powers[m + n], mm(powers[m], powers[n])) coefficient = compressed_coefficient(m + n) - compressed_coefficient(m) * compressed_coefficient(n) assert residual == scale(mode_projection, coefficient) assert abs(coefficient) <= m * n * delta if beta == -alpha: assert coefficient == (alpha ** (m + n) if m % 2 and n % 2 else 0) residual_rows.append(dict(model=tag, m=m, n=n, retained_residual=coefficient, rigorous_bound=m * n * delta)) # Compare actual matrix products for every ordered partition through n=10. schedule_rows = [] schedule_count = 0 for n in range(1, 11): largest_error = F(0) for parts in compositions(n): product_matrix = identity(4) for part in parts: product_matrix = mm(product_matrix, powers[part]) error = subtract(powers[n], product_matrix) coefficient = compressed_coefficient(n) - math.prod(compressed_coefficient(x) for x in parts) assert error == scale(mode_projection, coefficient) bound = F(n * n - sum(x * x for x in parts), 2) * delta assert abs(coefficient) <= bound largest_error = max(largest_error, abs(coefficient)) schedule_count += 1 schedule_rows.append(dict(model=tag, length=n, schedules=2 ** (n - 1), largest_error=largest_error, common_bound=F(n * (n - 1), 2) * delta)) details = dict(model=tag, alpha=str(alpha), beta=str(beta), plaquette_weights=list(map(str, weights)), transfer_eigenvalues=list(map(str, [1, alpha, beta, 0])), retained_one_step=list(map(str, [1, compressed_coefficient(1)])), retained_two_step=list(map(str, [1, compressed_coefficient(2)])), residual_two_step=str(delta), retained_one_step_is_fixed=(power(retained, 2) == retained), ordered_schedules_checked=schedule_count) return details, residual_rows, schedule_rows def matched_model_checks(): """Construct the two matched positive models in group value coordinates.""" group = list(itertools.product(range(2), repeat=2)) index = {g: i for i, g in enumerate(group)} chi1 = [F((-1) ** x) for x, y in group] chi2 = [F((-1) ** y) for x, y in group] characters = [[F(1), chi1[i], chi2[i], chi1[i] * chi2[i]] for i in range(4)] mode_projection = [[(chi1[i] + chi2[i]) * (chi1[j] + chi2[j]) / 8 for j in range(4)] for i in range(4)] projection = [[F(1, 4) + mode_projection[i][j] for j in range(4)] for i in range(4)] retained, rows = [], [] for tag, alpha, beta in [("A", F(3, 8), F(1, 8)), ("B", F(1, 4), F(1, 4))]: weights = [1 + alpha * chi1[i] + beta * chi2[i] for i in range(4)] assert all(w > 0 for w in weights) and sum(weights) == 4 matrix = [[weights[index[(x[0] ^ y[0], x[1] ^ y[1])]] / 4 for y in group] for x in group] fourier = scale(mm(mm(transpose(characters), matrix), characters), F(1, 4)) assert fourier == diag([1, alpha, beta, 0]) assert all(0 <= fourier[i][i] <= 1 for i in range(4)) one_step = mm(mm(projection, matrix), projection) two_step = mm(mm(projection, power(matrix, 2)), projection) delta = subtract(two_step, power(one_step, 2)) # u=(chi1+chi2)/sqrt(2); its matrix element is rational because # the two factors of sqrt(2) cancel in the Haar inner product. u_scaled = [[chi1[i] + chi2[i]] for i in range(4)] one_element = scale(mm(mm(transpose(u_scaled), one_step), u_scaled), F(1, 8))[0][0] two_element = scale(mm(mm(transpose(u_scaled), two_step), u_scaled), F(1, 8))[0][0] delta_element = scale(mm(mm(transpose(u_scaled), delta), u_scaled), F(1, 8))[0][0] assert one_element == F(1, 4) assert two_element == (F(5, 64) if tag == "A" else F(1, 16)) assert delta_element == (F(1, 64) if tag == "A" else 0) assert delta == scale(mode_projection, delta_element) retained.append((one_step, two_step)) rows.append(dict(model=tag, alpha=alpha, beta=beta, vacuum_one_step=1, vacuum_two_step=1, retained_one_step=one_element, retained_two_step=two_element, residual_two_step=delta_element)) assert retained[0][0] == retained[1][0] assert subtract(retained[0][1], retained[1][1]) == scale(mode_projection, F(1, 64)) save_csv("matched_models.csv", list(rows[0]), rows) return {"models_checked": len(rows), "arithmetic": "exact rational", "one_step_operators_equal": True, "two_step_difference": "1/64", "residuals": {"A": "1/64", "B": "0"}} def error_constants(): """Plot exact size dependence, with independent finite-matrix checks. Panel (a) answers how the choice of norm changes the measured residual as independent channels are added: the operator norm is constant while the Hilbert-Schmidt norm grows. Panel (b) answers how the comparison graph changes the linear correction constant: a path and a complete graph have opposite size dependence under the same unweighted edge norm. Curves evaluate the paper's exact formulas at integer channel/vertex counts; connecting lines only aid reading. The two panels compare scales. """ alpha = F(1, 4) rows = [] for k in range(1, 65): hs_squared = (1 + alpha ** 4) ** k - 1 rows.append(dict(channels=k, retained_dimension=2 ** k, operator_residual=str(alpha ** 2), hs_residual_squared=str(hs_squared), hs_residual=math.sqrt(float(hs_squared)))) # Independently form the Kronecker operators for the small cases. if k <= 7: two = np.array([[1.0]]) one = np.array([[1.0]]) for _ in range(k): two = np.kron(two, np.diag([1, float(alpha ** 2)])) one = np.kron(one, np.diag([1, 0])) difference = two - one assert math.isclose(np.linalg.norm(difference, 2), float(alpha ** 2), rel_tol=1e-12) assert math.isclose(np.linalg.norm(difference, "fro") ** 2, float(hs_squared), rel_tol=1e-12) save_csv("interface_constants.csv", list(rows[0]), rows) graph_rows = [] for n in range(2, 65): path_incidence = np.zeros((n - 1, n)) for i in range(n - 1): path_incidence[i, i] = -1 path_incidence[i, i + 1] = 1 spectrum = np.linalg.eigvalsh(path_incidence.T @ path_incidence) gap = 2 - 2 * math.cos(math.pi / n) assert math.isclose(spectrum[1], gap, rel_tol=1e-10) complete_edges = list(itertools.combinations(range(n), 2)) complete_incidence = np.zeros((len(complete_edges), n)) for edge, (i, j) in enumerate(complete_edges): complete_incidence[edge, i] = -1 complete_incidence[edge, j] = 1 complete_laplacian = complete_incidence.T @ complete_incidence assert np.array_equal(complete_laplacian, n * np.eye(n) - np.ones((n, n))) complete_spectrum = np.linalg.eigvalsh(complete_laplacian) assert math.isclose(complete_spectrum[0], 0, abs_tol=1e-12) assert np.allclose(complete_spectrum[1:], n, rtol=1e-12, atol=1e-12) # Standard unweighted Euclidean sum over all tested edges. graph_rows.append(dict(vertices=n, path_edges=n - 1, complete_edges=n * (n - 1) // 2, path_gap=gap, complete_gap=n, path_correction_constant=1 / math.sqrt(gap), complete_correction_constant=1 / math.sqrt(n))) save_csv("test_graph_constants.csv", list(graph_rows[0]), graph_rows) # The article has 25 mm margins on A4 and uses T1-encoded Latin Modern. # At this exact 160 mm width, 10.95 TeX points match its 11 pt body text. plt.rcParams.update({"text.usetex": True, "text.latex.preamble": r"\usepackage[T1]{fontenc}\usepackage{lmodern}", "font.family": "serif", "font.serif": ["Computer Modern Roman"], "font.size": 10.95, "axes.titlesize": 10.95, "axes.labelsize": 10.95, "xtick.labelsize": 10.95, "ytick.labelsize": 10.95, "legend.fontsize": 10.95, "axes.spines.top": False, "axes.spines.right": False, "pdf.fonttype": 42, "ps.fonttype": 42}) fig, axes = plt.subplots(1, 2, figsize=(160 / 25.4, 3.05), layout="constrained") fig.get_layout_engine().set(w_pad=0.025, h_pad=0.045, wspace=0.08) k = np.array([r["channels"] for r in rows]) axes[0].plot(k, [r["hs_residual"] for r in rows], color="#16697a", label="Hilbert--Schmidt") axes[0].axhline(float(alpha ** 2), color="#c66b30", linestyle="--", label="Operator") axes[0].set(xlabel=r"Interface channels $k$", ylabel="Two-step residual norm", xlim=(1, 64), title=r"(a) Independent channels, $a=1/4$") axes[0].legend(frameon=False, loc="upper left") n = np.array([r["vertices"] for r in graph_rows]) axes[1].plot(n, [r["path_correction_constant"] for r in graph_rows], color="#16697a", label="Path") axes[1].plot(n, [r["complete_correction_constant"] for r in graph_rows], color="#c66b30", linestyle="--", label="Complete graph") axes[1].set(xlabel=r"Vertices $N$", ylabel=r"Correction constant $1/\sqrt{\lambda_2}$", xlim=(2, 64), yscale="log", title="(b) Unweighted edge residuals") axes[1].legend(frameon=False, loc="center right") for ax in axes: ax.grid(alpha=0.2, linewidth=0.5) fig.savefig(FIGURES / "error_constants.pdf", metadata={"CreationDate": None}) plt.close(fig) return {"path_spectra_checked": len(graph_rows), "complete_graph_spectra_checked": len(graph_rows), "vertex_range": [2, 64], "independent_channel_matrices_checked": 7} def sharpness_checks(): """Check the sharpness witness from its two-by-two matrix, exactly.""" rows = [] partitions = [(1, 1), (2, 3), (1, 2, 3), (2, 2, 2), (1,) * 6] for parts in partitions: n = sum(parts) constant = F(n * n - sum(x * x for x in parts), 2) for denominator in (100, 1000, 10000): epsilon = F(1, denominator) matrix = [[1 - epsilon, epsilon], [epsilon, 1 - epsilon]] compressed = {j: power(matrix, j)[0][0] for j in range(1, n + 1)} delta = compressed[2] - compressed[1] ** 2 assert delta == epsilon ** 2 error = compressed[n] - math.prod(compressed[j] for j in parts) ratio = error / delta assert 0 <= ratio <= constant assert abs(ratio - constant) <= 256 * epsilon rows.append(dict(partition="+".join(map(str, parts)), epsilon=str(epsilon), residual=str(delta), schedule_error=str(error), ratio=float(ratio), optimal_constant=str(constant))) save_csv("sharpness_witness.csv", list(rows[0]), rows) return len(rows) def main(): count = finite_group_checks() all_residuals, all_schedules, models = [], [], [] for alpha, beta, tag in [(F(1, 4), F(-1, 4), "fixed_retained_operator"), (F(3, 8), F(1, 8), "positive_semidefinite_transfer")]: details, residuals, schedules = transfer_checks(alpha, beta, tag) models.append(details) all_residuals += residuals all_schedules += schedules save_csv("blocking_residuals.csv", list(all_residuals[0]), all_residuals) save_csv("blocking_schedules.csv", list(all_schedules[0]), all_schedules) matched_models = matched_model_checks() graph_checks = error_constants() sharpness_count = sharpness_checks() report = {"status": "all_checks_passed", "arithmetic": "exact rational for finite state sums and schedules", "bf_equalities_checked": count, "groups": ["Z2", "Z3", "Z2xZ2", "S3"], "models": models, "residual_pairs_checked": len(all_residuals), "ordered_schedules_checked": sum(x["ordered_schedules_checked"] for x in models), "interface_channels_plotted": 64, "test_graph_max_vertices": 64, "matched_models": matched_models, "graph_and_channel_checks": graph_checks, "sharpness_witnesses_checked": sharpness_count, "python_version": platform.python_version(), "numpy_version": np.__version__, "matplotlib_version": matplotlib.__version__} (DATA / "verification.json").write_text(json.dumps(report, indent=2) + "\n", encoding="utf-8") print(json.dumps(report, indent=2)) if __name__ == "__main__": main()