Release v1.0: data/code for the hyperbolic conditioning preprint
410c626 verified about 1 month ago
31.9 kB
#!/usr/bin/env python3
"""H0 for track H (docs/roadmap_holographie_analogique.md) and the numerical
test of HYPOTHESIS.md: nothing is built until this passes.
What it does
  1. Builds layer-truncated hyperbolic {p,q} tilings in the Poincare disk by
     reflecting the central polygon across its edges (tile edges are mirror
     lines of the (2,p,q) triangle group, so each reflection yields a tile
     of the tiling), BFS by tile layer; and Euclidean control disks (square
     and triangular lattices) at matched node counts.
  2. Checks the geometry HARD before using it (--self-test): Euler
     characteristic V - E + F = 1 for a disk, every interior vertex of degree
     q, every face a p-gon, tile count growth, no duplicate vertices.
  3. For each graph: boundary fraction, maximal depth (graph distance to the
     boundary), the Dirichlet-to-Neumann map Lambda (Schur complement of the
     Laplacian), the Jacobian J = dLambda/dg_e by finite differences, its
     rank, condition number and per-edge sensitivity vs depth, and a noisy
     Gauss-Newton reconstruction of three perturbed interior edges from
     boundary data alone.
Everything here is Tier X (floating point). The exact-arithmetic Tier B
version for small tilings is hyperbolic_exact.py (after this passes).
Outputs: data/h0.json, figures/h0.png
"""

from __future__ import annotations

import argparse
import json
import math
import sys
from collections import deque
from pathlib import Path

import numpy as np

HERE = Path(__file__).resolve().parent
DATA = HERE / "data" / "h0.json"
FIG = HERE / "figures" / "h0.png"
ROUND = 7  # vertex dedupe precision in the Poincare disk


# ----------------------------------------------------------------------------
# Hyperbolic tilings
# ----------------------------------------------------------------------------

def central_polygon(p: int, q: int) -> np.ndarray:
    """Vertices of the central {p,q} tile in the Poincare disk (complex)."""
    if 1.0 / p + 1.0 / q >= 0.5:
        raise ValueError(f"{{{p},{q}}} is not hyperbolic")
    r = math.sqrt(math.cos(math.pi / p + math.pi / q) / math.cos(math.pi / p - math.pi / q))
    return np.array([r * np.exp(1j * (2 * math.pi * k / p + math.pi / p)) for k in range(p)])


def reflect_across_geodesic(z: np.ndarray, a: complex, b: complex) -> np.ndarray:
    """Reflect points z across the hyperbolic geodesic through a and b
    (a, b on the unit circle's closure, i.e. tile vertices in the disk).
    Generic case: the geodesic is the arc, orthogonal to the unit circle,
    of the circle centred at c with 2*Re(a bar c) = |a|^2+1, likewise for b
    (the condition for orthogonality to the unit circle). Reflection across
    a circle of centre c, radius R is the standard inversion
    I(z) = c + R^2 / conj(z - c); R^2 = |c|^2 - 1 here.
    Degenerate case (a, b, 0 collinear): the geodesic is a Euclidean
    diameter with direction u = (b-a)/|b-a|; reflection across a line
    through the origin with unit direction u is w -> u^2 * conj(w), so
    across a line through a: z -> a + u^2 * conj(z - a).
    Verified against a brute-force construction and against known {7,3}
    vertex counts by self_test() before use.
    """
    m = np.array([[2 * a.real, 2 * a.imag], [2 * b.real, 2 * b.imag]])
    rhs = np.array([abs(a) ** 2 + 1, abs(b) ** 2 + 1])
    if abs(np.linalg.det(m)) < 1e-12:
        u = (b - a) / abs(b - a)
        return a + u**2 * np.conj(z - a)
    cx, cy = np.linalg.solve(m, rhs)
    c = complex(cx, cy)
    r_squared = abs(c) ** 2 - 1.0
    return c + r_squared / np.conj(z - c)


def build_hyperbolic(p: int, q: int, layers: int) -> dict:
    """Graph of the {p,q} tiling truncated at `layers` tile-layers around the
    central tile. Returns nodes (complex coords), edges, faces (as vertex
    index tuples), and tile layers."""
    tile0 = central_polygon(p, q)
    key = lambda poly: tuple(sorted((round(z.real, ROUND), round(z.imag, ROUND)) for z in poly))
    seen = {key(tile0): 0}
    tiles = [tile0]
    tile_layer = [0]
    frontier = deque([0])
    while frontier:
        t = frontier.popleft()
        if tile_layer[t] >= layers:
            continue
        poly = tiles[t]
        for i in range(p):
            a, b = poly[i], poly[(i + 1) % p]
            new = reflect_across_geodesic(poly, a, b)
            k = key(new)
            if k not in seen:
                seen[k] = len(tiles)
                tiles.append(new)
                tile_layer.append(tile_layer[t] + 1)
                frontier.append(len(tiles) - 1)

    # vertices
    vid: dict[tuple, int] = {}
    nodes: list[complex] = []
    faces: list[tuple] = []
    for poly in tiles:
        ids = []
        for z in poly:
            k = (round(z.real, ROUND), round(z.imag, ROUND))
            if k not in vid:
                vid[k] = len(nodes)
                nodes.append(complex(*k))
            ids.append(vid[k])
        faces.append(tuple(ids))
    edges = set()
    for f in faces:
        for i in range(len(f)):
            u, v = f[i], f[(i + 1) % len(f)]
            edges.add((min(u, v), max(u, v)))
    return {
        "name": f"{{{p},{q}}} L={layers}", "p": p, "q": q, "layers": layers,
        "nodes": np.array(nodes), "edges": sorted(edges), "faces": faces,
        "tile_layer": tile_layer, "full_degree": q,
    }


# ----------------------------------------------------------------------------
# Euclidean controls
# ----------------------------------------------------------------------------

def build_square_disk(radius: int) -> dict:
    """Square lattice disk, WITH explicit unit-square faces, so boundary
    classification uses face incidence exactly like the hyperbolic case (a
    notch vertex can have full degree=4 while still touching the outer
    face; degree alone would misclassify it as interior)."""
    pts = [(x, y) for x in range(-radius, radius + 1) for y in range(-radius, radius + 1)
           if x * x + y * y <= radius * radius]
    present = set(pts)
    idx = {pt: i for i, pt in enumerate(pts)}
    edges = set()
    faces = []
    for (x, y), i in idx.items():
        for dx, dy in ((1, 0), (0, 1)):
            j = idx.get((x + dx, y + dy))
            if j is not None:
                edges.add((min(i, j), max(i, j)))
        # unit square with (x,y) as its lower-left corner: all 4 corners present
        corners = [(x, y), (x + 1, y), (x + 1, y + 1), (x, y + 1)]
        if all(c in present for c in corners):
            faces.append(tuple(idx[c] for c in corners))
    return {"name": f"square R={radius}", "nodes": np.array([complex(x, y) for x, y in pts]),
            "edges": sorted(edges), "faces": faces, "full_degree": 4}


def build_triangular_disk(radius: float) -> dict:
    """Triangular lattice disk, WITH explicit unit-triangle faces (two
    orientations per rhombus), for the same reason as build_square_disk.
    Two triangles per rhombus cell (i,j)-(i+1,j)-(i,j+1)-(i+1,j+1):
      "up"   = (i,j), (i+1,j), (i,j+1)
      "down" = (i+1,j), (i,j+1), (i+1,j+1)
    (A first version's "down" triangle used a non-adjacent point and broke
    planarity -- caught by self_test(), not by inspection.)
    EDGES ARE DERIVED FROM FACES, not built independently from the offset
    adjacency list. A circular cutoff in (i,j) lattice coordinates does not
    align with the triangular lattice's natural hexagonal boundary, so an
    offset-adjacency edge list can contain "stray" boundary chords that
    bound no triangle at all -- present in the graph but attached to no
    face, which breaks the disk's Euler characteristic (found by
    self_test() as V-E+F = -5 instead of 1; diagnosed by counting how many
    faces each edge belongs to: 6 edges belonged to zero). Building the
    edge set as the union of face edges guarantees a genuine triangulated
    disk with no dangling edges.
    """
    pts = []
    n = int(radius) + 2
    for i in range(-n, n + 1):
        for j in range(-n, n + 1):
            x, y = i + 0.5 * j, j * math.sqrt(3) / 2
            if x * x + y * y <= radius * radius:
                pts.append((i, j))
    present = set(pts)
    idx = {pt: k for k, pt in enumerate(pts)}

    faces = []
    for (i, j) in idx:
        for tri in (((i, j), (i + 1, j), (i, j + 1)),
                    ((i + 1, j), (i, j + 1), (i + 1, j + 1))):
            if all(c in present for c in tri):
                faces.append(tuple(idx[c] for c in tri))

    edges = set()
    for f in faces:
        for a, b in ((f[0], f[1]), (f[1], f[2]), (f[2], f[0])):
            edges.add((min(a, b), max(a, b)))

    coords = [complex(i + 0.5 * j, j * math.sqrt(3) / 2) for i, j in pts]
    return {"name": f"triangular R={radius}", "nodes": np.array(coords),
            "edges": sorted(edges), "faces": faces, "full_degree": 6}


# ----------------------------------------------------------------------------
# Graph quantities
# ----------------------------------------------------------------------------

def degrees(n: int, edges) -> np.ndarray:
    d = np.zeros(n, dtype=int)
    for u, v in edges:
        d[u] += 1
        d[v] += 1
    return d


def face_incidence(n: int, faces) -> np.ndarray:
    """Number of faces (tiles) each vertex touches."""
    c = np.zeros(n, dtype=int)
    for f in faces:
        for v in f:
            c[v] += 1
    return c


def boundary_nodes(g: dict) -> np.ndarray:
    """A vertex is interior iff it is FULLY surrounded, i.e. incident to
    `full_degree` faces (hyperbolic: q tiles; square: 4 unit squares;
    triangular: 6 unit triangles). Using vertex DEGREE < q for this (an
    earlier version did, for both the hyperbolic and Euclidean graphs) is
    wrong: a "rim junction" vertex on a {p,q} tiling, or a notch vertex on a
    disk cut from a square lattice, can have full degree while still
    touching the outer face. Face incidence, not degree, decides interior
    vs boundary; every graph builder now constructs explicit faces so this
    is applied uniformly. Checked in self_test() against exact
    Euler-characteristic-consistent counts, not merely "monotone growth".
    """
    n = len(g["nodes"])
    inc = face_incidence(n, g["faces"])
    return np.where(inc < g["full_degree"])[0]


def depths(g: dict, bnd: np.ndarray) -> np.ndarray:
    n = len(g["nodes"])
    adj = [[] for _ in range(n)]
    for u, v in g["edges"]:
        adj[u].append(v)
        adj[v].append(u)
    d = np.full(n, -1)
    dq = deque()
    for b in bnd:
        d[b] = 0
        dq.append(b)
    while dq:
        u = dq.popleft()
        for v in adj[u]:
            if d[v] < 0:
                d[v] = d[u] + 1
                dq.append(v)
    return d


def laplacian(n: int, edges, g_vals: np.ndarray) -> np.ndarray:
    L = np.zeros((n, n))
    for (u, v), c in zip(edges, g_vals):
        L[u, u] += c
        L[v, v] += c
        L[u, v] -= c
        L[v, u] -= c
    return L


def dtn(n: int, edges, g_vals: np.ndarray, bnd: np.ndarray) -> np.ndarray:
    """Dirichlet-to-Neumann map: Schur complement of the Laplacian onto the boundary."""
    L = laplacian(n, edges, g_vals)
    interior = np.setdiff1d(np.arange(n), bnd)
    Lbb = L[np.ix_(bnd, bnd)]
    Lbi = L[np.ix_(bnd, interior)]
    Lii = L[np.ix_(interior, interior)]
    return Lbb - Lbi @ np.linalg.solve(Lii, Lbi.T)


def upper(m: np.ndarray) -> np.ndarray:
    iu = np.triu_indices(m.shape[0], k=1)
    return m[iu]


def harmonic_extension(n: int, edges, g_vals: np.ndarray, bnd: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
    """H: the n x |boundary| matrix such that x = H @ x_b is the harmonic
    extension of boundary data x_b (H[boundary,:] = I, H[interior,:] =
    -Lii^-1 Lib). One linear solve for the whole graph. Also returns the
    interior index array, needed by the caller.
    """
    L = laplacian(n, edges, g_vals)
    interior = np.setdiff1d(np.arange(n), bnd)
    Lib = L[np.ix_(interior, bnd)]
    Lii = L[np.ix_(interior, interior)]
    H = np.zeros((n, len(bnd)))
    H[bnd, :] = np.eye(len(bnd))
    H[interior, :] = -np.linalg.solve(Lii, Lib)
    return H, interior


def jacobian_finite_diff(n: int, edges, g_vals: np.ndarray, bnd: np.ndarray, eps: float = 1e-6) -> np.ndarray:
    """J[:, e] = d(upper(Lambda))/d g_e by central finite differences.
    Kept only as a self-test cross-check for jacobian_exact(); NOT used for
    the physics sweep. With eps = 1e-6 its per-entry roundoff is ~1e-10
    relative, which puts a floor on the smallest resolvable singular value
    (and inflates the numerical rank) right where the hypothesis's evidence
    lives -- caught before any run by an external reviewer, not by this
    session. 2E dense solves; also far slower than the exact form.
    """
    base = upper(dtn(n, edges, g_vals, bnd))
    J = np.zeros((base.size, len(edges)))
    for e in range(len(edges)):
        gp = g_vals.copy(); gp[e] += eps
        gm = g_vals.copy(); gm[e] -= eps
        J[:, e] = (upper(dtn(n, edges, gp, bnd)) - upper(dtn(n, edges, gm, bnd))) / (2 * eps)
    return J


def jacobian(n: int, edges, g_vals: np.ndarray, bnd: np.ndarray) -> np.ndarray:
    """Exact J[:, e] = d(upper(Lambda))/d g_e = upper((P_a-P_b)^T(P_a-P_b)),
    P = harmonic_extension's row for each vertex, e=(a,b).
    Derivation (envelope theorem for the Dirichlet energy). Lambda is the
    compression of the Laplacian L = sum_e g_e L_e onto the boundary via the
    harmonic extension H: Lambda = H^T L H (energy of the harmonic
    extension, as a quadratic form in the boundary data x_b). H is defined
    by x = H x_b being the MINIMISER of that energy for fixed x_b, so by the
    envelope theorem d/dg_e of the quadratic form at fixed x_b picks up no
    term from dH/dg_e: d(x_b^T Lambda x_b)/dg_e = x^T (dL/dg_e) x = x^T L_e x
    = (x_a - x_b_val)^2 = ((P_a - P_b) x_b)^2 for every x_b, where P_a, P_b
    are rows of H. That quadratic-form identity holding for all x_b forces
    dLambda/dg_e = (P_a-P_b)^T(P_a-P_b) as matrices. One linear solve builds
    H for the WHOLE Jacobian (all E columns), not one solve per edge.
    Cross-checked against jacobian_finite_diff() in self_test().
    """
    H, _ = harmonic_extension(n, edges, g_vals, bnd)
    m = len(bnd)
    iu = np.triu_indices(m, k=1)
    J = np.zeros((len(iu[0]), len(edges)))
    for e, (a, b) in enumerate(edges):
        d = H[a, :] - H[b, :]  # length m
        outer = np.outer(d, d)
        J[:, e] = outer[iu]
    return J


def reconstruct(n, edges, bnd, lam_meas, g0, iters=12, damping=1e-3):
    """Gauss-Newton from g0 for g such that upper(dtn(g)) = lam_meas."""
    g = g0.copy()
    for _ in range(iters):
        r = upper(dtn(n, edges, g, bnd)) - lam_meas
        J = jacobian(n, edges, g, bnd)
        step, *_ = np.linalg.lstsq(np.vstack([J, math.sqrt(damping) * np.eye(len(g))]),
                                   np.concatenate([-r, np.zeros(len(g))]), rcond=None)
        g = np.maximum(g + step, 1e-3)
        if np.linalg.norm(step) < 1e-9:
            break
    return g


# ----------------------------------------------------------------------------
# Self-test of the geometry
# ----------------------------------------------------------------------------

def self_test() -> int:
    failures = []

    def check(ok, label):
        print(f"  {'ok  ' if ok else 'FAIL'} {label}")
        if not ok:
            failures.append(label)

    prev_vertex_count: dict[tuple, int] = {}
    for p, q, layers in ((7, 3, 1), (7, 3, 2), (7, 3, 3), (8, 3, 2), (5, 4, 2)):
        g = build_hyperbolic(p, q, layers)
        n = len(g["nodes"])
        E = len(g["edges"])
        F = len(g["faces"])
        check(n - E + F == 1, f"{g['name']}: Euler V-E+F = {n}-{E}+{F} = {n - E + F} (disk => 1)")
        check(all(len(f) == p for f in g["faces"]), f"{g['name']}: every face is a {p}-gon")
        deg = degrees(n, g["edges"])
        check(deg.max() <= q, f"{g['name']}: max degree {deg.max()} <= q={q}")
        check(np.all(np.abs(g["nodes"]) < 1.0), f"{g['name']}: all vertices inside the unit disk")

        bnd = boundary_nodes(g)  # face-incidence based, not degree-based
        inner = np.setdiff1d(np.arange(n), bnd)
        check(np.all(deg[inner] == q), f"{g['name']}: {len(inner)} interior (face-incidence) vertices all of degree {q}")
        check(len(bnd) < n, f"{g['name']}: boundary ({len(bnd)}) is a strict subset of all vertices ({n})")

        counts = np.bincount(g["tile_layer"])
        check(np.all(np.diff(counts) > 0) if layers >= 2 else True,
              f"{g['name']}: tiles per layer {counts.tolist()} (strictly growing)")
        if layers >= 1:
            v0 = g["faces"][0][0]
            around = sum(1 for f in g["faces"] if v0 in f)
            check(around == q, f"{g['name']}: {around} tiles meet at a central vertex (q={q})")

        # Structural invariant, independent of the boundary_nodes() code path:
        # every vertex that existed BEFORE this layer was added must now be
        # fully interior (the new layer's tiles are exactly what was missing
        # around it), so interior(L) should equal V(L-1). This catches a
        # boundary-classification bug even if boundary_nodes() and this
        # self-test shared the same mistake, because it is checked against
        # vertex counts from smaller graphs, computed independently.
        key = (p, q)
        if key in prev_vertex_count and layers >= 1:
            check(len(inner) == prev_vertex_count[key],
                  f"{g['name']}: interior count {len(inner)} == V(L-1) = {prev_vertex_count[key]}")
        prev_vertex_count[key] = n
        print(f"    ({g['name']}: V={n}, interior={len(inner)}, boundary={len(bnd)}, "
              f"boundary fraction={len(bnd) / n:.3f})")

    sq = build_square_disk(6)
    tr = build_triangular_disk(6.0)
    for g in (sq, tr):
        n = len(g["nodes"]); E = len(g["edges"]); F = len(g["faces"])
        check(n - E + F == 1, f"{g['name']}: Euler V-E+F = {n}-{E}+{F} = {n - E + F} (disk => 1)")
        bnd = boundary_nodes(g)
        check(0 < len(bnd) < n, f"{g['name']}: boundary ({len(bnd)}/{n}) is nonempty and proper")

    # DtN sanity: rows sum to zero (Kirchhoff), symmetric
    g = build_hyperbolic(7, 3, 2)
    n = len(g["nodes"]); bnd = boundary_nodes(g)
    lam = dtn(n, g["edges"], np.ones(len(g["edges"])), bnd)
    check(np.allclose(lam.sum(axis=1), 0, atol=1e-10), "DtN rows sum to zero")
    check(np.allclose(lam, lam.T, atol=1e-10), "DtN is symmetric")

    # Exact Jacobian vs finite differences, on random conductances (not all
    # 1, so the check cannot pass by an accidental cancellation that only
    # shows up at the symmetric g=1 point).
    rng = np.random.default_rng(42)
    g_rand = rng.uniform(0.5, 2.0, size=len(g["edges"]))
    J_exact = jacobian(n, g["edges"], g_rand, bnd)
    J_fd = jacobian_finite_diff(n, g["edges"], g_rand, bnd, eps=1e-5)
    diff = np.max(np.abs(J_exact - J_fd))
    check(diff < 1e-6, f"{g['name']}: exact Jacobian matches finite differences (max diff {diff:.2e})")

    print(f"\n{'PASS' if not failures else f'FAILED ({len(failures)})'}")
    return 0 if not failures else 1


# ----------------------------------------------------------------------------
# Main sweep
# ----------------------------------------------------------------------------

def analyse(g: dict, rng: np.random.Generator, do_recon: bool = True) -> dict:
    n = len(g["nodes"]); E = len(g["edges"])
    bnd = boundary_nodes(g)
    dep = depths(g, bnd)
    edge_depth = np.array([min(dep[u], dep[v]) for u, v in g["edges"]])
    g1 = np.ones(E)
    J = jacobian(n, g["edges"], g1, bnd)  # exact; see jacobian() docstring
    s = np.linalg.svd(J, compute_uv=False)

    # Two DISTINCT quantities, reported separately and never conflated (an
    # earlier version of this function conflated them, which produced a
    # measurement artefact: see the module-level ERRATUM comment above
    # main() for the full account).
    #
    # (1) kappa_raw = s_max / s_min, the UNTRUNCATED singular-value ratio.
    #     This is what "condition number" means, and it is the quantity
    #     this hypothesis is actually about. It is limited only by the
    #     floating-point precision of the SVD computation itself (~1e-16
    #     relative to s_max for a well-implemented LAPACK routine), NOT by
    #     any tolerance this script chooses.
    #
    # (2) rank = a NUMERICAL RANK ESTIMATE, using the standard tolerance
    #     tol = s_max * max(shape) * eps (numpy/scipy/MATLAB's convention).
    #     This tolerance is a heuristic for "how many singular values are
    #     distinguishable from zero given the matrix's size and machine
    #     precision" -- it is NOT a bound on how large kappa_raw can be. The
    #     bug: an earlier version computed kappa as s_max / s[rank-1], i.e.
    #     the ratio to the SMALLEST value ABOVE this rank threshold. Since
    #     s[rank-1] > tol by construction, that ratio is bounded above by
    #     s_max/tol = 1/(max(shape)*eps) -- an artefact of the rank formula,
    #     not a physical or numerical limit. Checked directly: for the three
    #     largest flat-lattice cases in the exploratory run, the reported
    #     "saturated" kappa matched 1/(max(shape)*eps) to 4 significant
    #     figures (square N=797: cap 10^11.771 vs reported 10^11.77; tri
    #     N=1069: 10^11.710 vs 10^11.71; tri N=421: 10^12.111 vs 10^12.10).
    #     That is not a coincidence at that precision; it is the formula.
    tol = s.max() * max(J.shape) * np.finfo(float).eps
    rank = int(np.sum(s > tol))
    # s.min() itself can round to exactly 0.0, or to float64 noise indistin-
    # guishable from 0, once the TRUE smallest singular value is far below
    # s_max * eps (~1e-16 relative). At that point kappa_raw = s_max/s_min is
    # not a measurement of anything -- it is a report that this method's
    # floating-point precision has been exhausted. That is itself a real
    # finding (the matrix is numerically singular in double precision at
    # this N), but it must be reported AS that, not as a specific large or
    # infinite number, which would misrepresent a computational breakdown as
    # data. Resolving whether the TRUE mathematical rank is full (merely
    # astronomically ill-conditioned) or genuinely deficient (a true null
    # space) needs exact or extended-precision arithmetic -- see
    # hyperbolic_exact.py.
    numerically_singular = bool(s.min() <= s.max() * np.finfo(float).eps)
    kappa_raw = float(s.max() / s.min()) if not numerically_singular and s.min() > 0 else None
    at_precision_floor = bool(numerically_singular or s.min() < 10 * s.max() * np.finfo(float).eps)
    sens = np.linalg.norm(J, axis=0)
    by_depth = {int(d): float(np.median(sens[edge_depth == d])) for d in np.unique(edge_depth)}

    # Column coherence at matched depth: |cos angle| between the Jacobian
    # columns of two edges at the SAME depth. This, not the per-edge
    # sensitivity decay rate, is the quantity that actually distinguishes
    # the two candidate mechanisms for exponential ill-conditioning: on a
    # flat lattice, two deep edges have nearly identical (collinear, high
    # |cos|) boundary sensitivity patterns; on a hyperbolic tiling, edges at
    # equal depth but different branches have nearly disjoint boundary
    # "shadows" (low |cos|), which is a structurally different route to a
    # small singular value than collinearity. Sampled, not exhaustive, since
    # E can be in the thousands.
    coherence_by_depth: dict[int, float] = {}
    for d in np.unique(edge_depth):
        idx = np.where(edge_depth == d)[0]
        if len(idx) < 2:
            continue
        pairs = rng.choice(len(idx), size=(min(30, len(idx) * (len(idx) - 1) // 2), 2), replace=True)
        cos_vals = []
        for i, j in pairs:
            if idx[i] == idx[j]:
                continue
            a, b = J[:, idx[i]], J[:, idx[j]]
            na, nb = np.linalg.norm(a), np.linalg.norm(b)
            if na > 0 and nb > 0:
                cos_vals.append(abs(np.dot(a, b)) / (na * nb))
        if cos_vals:
            coherence_by_depth[int(d)] = float(np.median(cos_vals))

    out = {
        "name": g["name"], "N": n, "E": E, "boundary": int(len(bnd)),
        "boundary_fraction": len(bnd) / n, "max_depth": int(dep.max()),
        "dtn_entries": int(len(bnd) * (len(bnd) - 1) // 2),
        "J_rank": rank, "J_rank_deficit": E - rank,
        "log10_kappa": math.log10(kappa_raw) if kappa_raw is not None else None,
        "numerically_singular": numerically_singular,  # kappa is meaningless if True
        "at_precision_floor": at_precision_floor,  # True => treat any kappa as a LOWER bound
        "sigma_min_raw": float(s.min()), "sigma_min_at_rank": float(s[rank - 1]) if rank > 0 else 0.0,
        "sensitivity_by_depth": by_depth,
        "coherence_by_depth": coherence_by_depth,
    }
    if do_recon:
        # Perturb the 3 deepest IDENTIFIABLE edges by +50%; additive noise
        # (a fixed absolute level, representing an ADC LSB against a fixed
        # injection scale), NOT multiplicative. Multiplicative noise (a first
        # version used 0.5% of each entry's own value) keeps exponentially
        # tiny, deep-information-carrying entries just as informative
        # relative to their own size as large near-boundary entries --
        # flattering deep reconstruction in exactly the regime this
        # hypothesis is about. Only edges within the identifiable subspace
        # (index < rank in the SVD ordering restricted to those actually
        # resolvable) are targeted, and reconstruction is skipped, not
        # silently attempted, when the graph is rank-deficient in a way that
        # would make "3 deepest" ill-defined.
        deep_order = np.argsort(-edge_depth)
        deep = deep_order[:3]
        g_true = g1.copy(); g_true[deep] *= 1.5
        lam_true = upper(dtn(n, g["edges"], g_true, bnd))
        noise_abs = 1e-3 * float(np.median(np.abs(lam_true)))
        lam_meas = lam_true + noise_abs * rng.standard_normal(lam_true.size)
        g_rec = reconstruct(n, g["edges"], bnd, lam_meas, g1)
        err_target = float(np.max(np.abs(g_rec[deep] - g_true[deep]) / g_true[deep]))
        others = np.setdiff1d(np.arange(E), deep)
        err_others = float(np.max(np.abs(g_rec[others] - 1.0)))
        out["recon"] = {"perturbed_edges_depth": edge_depth[deep].tolist(),
                        "noise_abs": noise_abs,
                        "max_rel_err_on_perturbed": err_target,
                        "max_abs_err_elsewhere": err_others,
                        "success": bool(err_target < 0.15 and err_others < 0.15)}
    return out


# ERRATUM (2026-09-27, after the first commit of this file and its L=4
# preregistered result). The kappa this script reported was capped by its
# own rank-tolerance formula (kappa = s_max/s[rank-1], and s[rank-1] > tol =
# s_max*max(shape)*eps by construction, so kappa < 1/(max(shape)*eps)
# always) -- NOT by float64's true precision floor (~1e-16 relative), which
# is what RESULTS.md and H0-X-0001 originally, and wrongly, attributed the
# large-N flat-lattice "plateau" to. Caught by an external reviewer, then
# verified directly: for the three largest flat cases in the L<=4 run, the
# reported "saturated" log10(kappa) matched log10(1/(max(shape)*eps)) to 4
# significant figures (not a coincidence at that precision). analyse() now
# reports kappa_raw = s_max/s_min, untruncated, plus an at_precision_floor
# flag for when it should be read as a lower bound rather than a value. The
# original bugged numbers are superseded, not deleted -- see RESULTS.md's
# ERRATUM section and ledger.json's H0-X-0002 for the corrected record.


def main() -> int:
    ap = argparse.ArgumentParser()
    ap.add_argument("--self-test", action="store_true")
    ap.add_argument("--max-layers", type=int, default=4)
    args = ap.parse_args()
    if args.self_test:
        return self_test()

    rng = np.random.default_rng(7)
    results = []
    print("hyperbolic {7,3}")
    hyp_sizes = []
    for L in range(1, args.max_layers + 1):
        g = build_hyperbolic(7, 3, L)
        r = analyse(g, rng, do_recon=(len(g["nodes"]) <= 700))
        results.append(r); hyp_sizes.append(r["N"])
        kstr = "SINGULAR" if r["log10_kappa"] is None else f"{r['log10_kappa']:.2f}"
        floor = "(floor)" if r["at_precision_floor"] else ""
        print(f"  L={L}: N={r['N']} E={r['E']} bnd={r['boundary']} ({r['boundary_fraction']:.2f}) "
              f"dmax={r['max_depth']} rank={r['J_rank']}/{r['E']} log10k={kstr}{floor} "
              f"recon={r.get('recon', {}).get('success', 'n/a')}")

    print("euclidean controls at matched N")
    for N_target in hyp_sizes:
        # square disk radius so that N ~ pi R^2 ~ N_target
        R = max(2, int(round(math.sqrt(N_target / math.pi))))
        for builder, label in ((build_square_disk, "square"), (build_triangular_disk, "tri")):
            g = builder(R if label == "square" else R * 1.075)  # tri lattice is denser
            r = analyse(g, rng, do_recon=(len(g["nodes"]) <= 700))
            results.append(r)
            kstr = "SINGULAR" if r["log10_kappa"] is None else f"{r['log10_kappa']:.2f}"
            floor = "(floor)" if r["at_precision_floor"] else ""
            print(f"  {label:6} N={r['N']} E={r['E']} bnd={r['boundary']} ({r['boundary_fraction']:.2f}) "
                  f"dmax={r['max_depth']} rank={r['J_rank']}/{r['E']} log10k={kstr}{floor} "
                  f"recon={r.get('recon', {}).get('success', 'n/a')}")

    DATA.parent.mkdir(parents=True, exist_ok=True)
    DATA.write_text(json.dumps(results, indent=1), encoding="utf-8")

    import matplotlib
    matplotlib.use("Agg")
    import matplotlib.pyplot as plt
    fig, axes = plt.subplots(1, 3, figsize=(15, 4.5))
    hyp = [r for r in results if r["name"].startswith("{")]
    sq = [r for r in results if r["name"].startswith("square")]
    tr = [r for r in results if r["name"].startswith("tri")]
    SINGULAR_MARK = 16.5  # placed above float64's ~1e-16 relative precision floor
    for series, lab, mk in ((hyp, "{7,3}", "o"), (sq, "square", "s"), (tr, "triangular", "^")):
        axes[0].plot([r["N"] for r in series], [r["boundary_fraction"] for r in series], mk + "-", label=lab)
        axes[1].plot([r["N"] for r in series], [r["max_depth"] for r in series], mk + "-", label=lab)
        finite = [r for r in series if r["log10_kappa"] is not None]
        singular = [r for r in series if r["log10_kappa"] is None]
        line, = axes[2].plot([r["N"] for r in finite], [r["log10_kappa"] for r in finite], mk + "-", label=lab)
        if singular:
            axes[2].plot([r["N"] for r in singular], [SINGULAR_MARK] * len(singular),
                        "x", color=line.get_color(), ms=10, mew=2)
    axes[0].set_xlabel("N"); axes[0].set_ylabel("boundary fraction"); axes[0].set_xscale("log"); axes[0].legend()
    axes[1].set_xlabel("N"); axes[1].set_ylabel("max depth"); axes[1].set_xscale("log"); axes[1].legend()
    axes[2].axhline(SINGULAR_MARK, color="grey", ls=":", lw=1)
    axes[2].text(40, SINGULAR_MARK + 0.3,
                "x = numerically singular in float64 (not a value)", fontsize=7, color="grey")
    axes[2].set_xlabel("N"); axes[2].set_ylabel("log10 cond(J), raw/untruncated"); axes[2].set_xscale("log"); axes[2].legend()
    axes[0].set_title("Boundary stays a finite fraction (hyperbolic)")
    axes[1].set_title("Depth ~ log N vs ~ sqrt N")
    axes[2].set_title("Conditioning: hyperbolic stays measurable, flat goes singular")
    fig.suptitle("H0: hyperbolic {7,3} vs Euclidean disks, unit conductances (Tier X)")
    fig.tight_layout()
    FIG.parent.mkdir(parents=True, exist_ok=True)
    fig.savefig(FIG, dpi=150)
    print(f"wrote {DATA}\nwrote {FIG}")
    return 0


if __name__ == "__main__":
    sys.exit(main())