"""The curvature weight of the Kuhn lattice: the measure form D of Theorem 5.4, with the checks of Theorems 5.3 and 5.5.

Paper: "Quantum Gravity from the Recognition Ledger: A Finite One-Loop Path Sum, Einstein Gravity at Long Wavelengths, and the
Massless-Graviton Problem", J. Washburn (2026), version 1.1, https://doi.org/10.5281/zenodo.23813297. Section and theorem numbers
below are that version's. "The earlier paper" is "Linearized Einstein gravity from Lorentzian Regge calculus on the Kuhn lattice, and the
divergence of its one-loop graviton limit", J. Washburn (2026), https://doi.org/10.5281/zenodo.23764491.

What it computes
  On the flat Kuhn lattice at kappa = 1, where an edge of class d in {0,1}^4, d != 0 (the fifteen classes), has squared length
  |d|(1 - |d|), the form through which the ledger's cost of curvature enters the path sum at Gaussian order (Theorem 5.4):
      D = 1/2 sum_h (A_h / n_h) delta'_h delta'_h^T,
  with h running over the hinges (triangles), A_h > 0 the flat area of h, delta'_h the derivative at the flat lattice of the
  deficit angle at h with respect to the squared edge lengths, and n_h the channel count. The count used is the dual cycle's:
  n_h is the number of four-simplices that contain h, 4 or 6, which the paper states is a definition. The formula holds under
  the paper's two identifications (Section 5.3) and the scaling that its phase already uses (Table 2): (I1) the cost of a
  hinge record reads its area as the action does: it adds over the pieces of the hinge, does not decrease with area, and a face
  of area lambda_rec^2 carries one cycle of the hinge's bridge (a model), with coupling 1 and the Regge deficit, so that it
  costs delta_h^2 / (2 n_h) at second order; and (I2) the Born weight of a geometry is the Gibbs weight exp(-C) of its ledger
  cost. The same geometric data, the flat hinge areas and the deficit derivatives, give the action's flat Hessian, which needs
  no identification: by Schlaefli's identity the Regge action S = sum_h A_h delta_h has the Hessian sum_h A'_h delta'_h^T at
  the flat lattice (Section 5.2). The program builds both, as Bloch matrices on the
  fifteen edge classes, from the exact flat dihedral angles and their derivatives, which it evaluates analytically at
  40 digits. The Bloch matrices are evaluated in double precision. Every number it prints is a numerical evaluation.

Checks (printed, each with PASS or FAIL)
  * The hinges: 50 shapes per cell, valences 4 and 6 (Section 5.3) adding to 240 = 24 simplices x 10 triangles; flat areas
    1/2, 1, 3/2 and sqrt3 (Section 6.4); weights A_h/n_h taking the values 1/12, 1/4 and sqrt3/4 (stated after Theorem 5.5).
    The flat dihedral angles add to 2 pi around every hinge, so every flat deficit is zero.
  * The action: the Hessian from Schlaefli's identity is Hermitian, annihilates the vertex moves (the earlier paper's
    Section 4.2), has the ten constant strains as the kernel at k = 0, and at each axis corner pi e_c annihilates the doubler
    v_c of the earlier paper's Proposition 6.1 (v_c(d) = 2^{|d|-3} if d_c = 0, else 0; |v_c|^2 = 31/16). With --compare and
    the earlier paper's program kuhn_bloch.py in the same folder, it also compares this Hessian with the one assembled there
    from the earlier paper's printed cell Hessian (Appendix B).
  * Theorem 5.3: the linearized deficits vanish on the constant strains and on the vertex moves; the signed census
    sum_h sigma_c(h) delta'_h of the doubler field (-1)^{x_c} 4 v_c equals N^4 (6 + 8 sqrt3) for each axis c.
  * Theorem 5.5 and the numbers after it: D vanishes on constant strains and vertex moves; on the complement of the vertex
    moves six eigenvalues of D(k) scale like |k|^4 while the other five stay gapped on long waves, the smallest near 4.92 at the
    centre; D is positive on
    each doubler, with damping coefficient v_c^dagger D v_c / |v_c|^2 = (10 + 3 sqrt3)/124 at all four axis corners, the value
    the paper states after Theorem 5.5 (here to double precision); with a = 1 (any a != 0
    serves, Theorem 5.2) the corner precision D - i a A is invertible, where A = B + Gamma Gamma^dagger is the earlier paper's
    gauge-fixed operator, which is singular there.
  * The remark after Theorem 5.6: at the four axis corners and at sampled waves the Bloch map of the linearized deficits has
    rank eleven, so its kernel is the four vertex moves.

Output
  With --json FILE, writes the data of the figures on the gravity pages: along the axis k = s pi e_0, 0 < s <= 1, the eleven
  eigenvalues of D(k) on the complement of the vertex moves; and, for the mode on that complement on which the action pushes
  back least, the action's eigenvalue and the weight's value on the same mode.

Requires Python 3.9+, numpy and mpmath. Run: python3 curvature_weight.py   [--json weight.json]   [--compare]
License: MIT. Copyright (c) 2026 Recognition Physics Institute.
"""
import argparse
import itertools
import json
import math
import os
import sys

import mpmath as mp
import numpy as np

mp.mp.dps = 40
PAIRS = [(0, 1), (0, 2), (0, 3), (0, 4), (1, 2), (1, 3), (1, 4), (2, 3), (2, 4), (3, 4)]     # edges of a simplex chain
TRIS = list(itertools.combinations(range(5), 3))                                            # its ten triangles
CLASSES = np.array([[(m >> a) & 1 for a in range(4)] for m in range(1, 16)], dtype=float)   # index = sum_a d_a 2^a - 1
SQRT3 = math.sqrt(3)


def class_index(d):
    return int(sum(int(d[a]) << a for a in range(4))) - 1


def flat_sq(n):
    """Flat squared length at kappa = 1 of an edge spanning n postings, |d|(1 - kappa |d|)."""
    return mp.mpf(int(n) * (1 - int(n)))


# One simplex along the chain v_0, ..., v_4 with v_j - v_i of size j - i. Squared lengths l_ij and the Gram matrix at v_0.
L0 = {(i, j): flat_sq(j - i) for i, j in PAIRS}


def gram(ell):
    G = mp.matrix(4, 4)
    for a in range(1, 5):
        for b in range(1, 5):
            G[a - 1, b - 1] = (ell[(0, a)] + ell[(0, b)] - (0 if a == b else ell[(min(a, b), max(a, b))])) / 2
    return G


def dual(Gi):
    """5 x 5 matrix of inner products of the gradients of the barycentric coordinates (rows sum to zero)."""
    N = mp.matrix(5, 5)
    for i in range(5):
        for j in range(5):
            if i and j:
                N[i, j] = Gi[i - 1, j - 1]
            elif i:
                N[i, j] = -sum(Gi[i - 1, b] for b in range(4))
            elif j:
                N[i, j] = -sum(Gi[a, j - 1] for a in range(4))
            else:
                N[i, j] = sum(Gi[a, b] for a in range(4) for b in range(4))
    return N


def dgram(edge):
    """Derivative of the Gram matrix with respect to the squared length of one edge (a constant matrix)."""
    E = mp.matrix(4, 4)
    i, j = edge
    if i == 0:
        for b in range(1, 5):
            E[j - 1, b - 1] += mp.mpf(1) / 2
            E[b - 1, j - 1] += mp.mpf(1) / 2
    else:
        E[i - 1, j - 1] -= mp.mpf(1) / 2
        E[j - 1, i - 1] -= mp.mpf(1) / 2
    return E


def angles_and_gradients():
    """Flat dihedral angles theta_t at the ten triangles, and d theta_t / d l_e, analytically: theta = arccos(-N_ij / sqrt(N_ii N_jj))
    with {i, j} the two vertices off t, and dN from d(G^-1) = -G^-1 dG G^-1."""
    Gi = gram(L0) ** -1
    N = dual(Gi)
    theta, grad = [], []
    dN = {e: dual(-Gi * dgram(e) * Gi) for e in PAIRS}
    for t in TRIS:
        i, j = [v for v in range(5) if v not in t]
        r = mp.sqrt(N[i, i] * N[j, j])
        c = -N[i, j] / r
        th = mp.acos(c)
        theta.append(th)
        row = []
        for e in PAIRS:
            d = dN[e]
            dc = -d[i, j] / r + N[i, j] / (2 * r) * (d[i, i] / N[i, i] + d[j, j] / N[j, j])
            row.append(-dc / mp.sin(th))
        grad.append(row)
    return theta, grad


def tri_area_and_gradient(l1, l2, l3):
    """Area of a timelike triangle from its squared lengths, 16 A^2 = l1^2 + l2^2 + l3^2 - 2(l1 l2 + l2 l3 + l3 l1), and dA/dl."""
    q = l1 * l1 + l2 * l2 + l3 * l3 - 2 * (l1 * l2 + l2 * l3 + l3 * l1)
    A = mp.sqrt(q) / 4
    return A, [(2 * l1 - 2 * l2 - 2 * l3) / (32 * A), (2 * l2 - 2 * l1 - 2 * l3) / (32 * A), (2 * l3 - 2 * l1 - 2 * l2) / (32 * A)]


THETA, DTHETA = angles_and_gradients()
DTH = np.array([[float(x) for x in row] for row in DTHETA])

# Hinges. A hinge is a triangle x, x + d1, x + d1 + d2 with disjoint nonempty d1, d2. INC[(d1, d2)] lists its incidences in the
# 24 simplices of one cell: (triangle index, chain vertices, position of the hinge's base vertex in the chain).
INC = {}
for perm in itertools.permutations(range(4)):
    V = [np.zeros(4, dtype=int)]
    for a in perm:
        v = V[-1].copy()
        v[a] = 1
        V.append(v)
    for ti, (i, j, k) in enumerate(TRIS):
        INC.setdefault((tuple(V[j] - V[i]), tuple(V[k] - V[j])), []).append((ti, V, i))
SHAPES = sorted(INC)
VALENCE = np.array([len(INC[s]) for s in SHAPES])
AREA, AREA_GRAD = [], []
for s in SHAPES:
    p, r = sum(s[0]), sum(s[1])
    A, g = tri_area_and_gradient(flat_sq(p), flat_sq(r), flat_sq(p + r))
    AREA.append(A)
    AREA_GRAD.append([float(x) for x in g])
WEIGHT = np.array([float(A) for A in AREA]) / VALENCE


def _deficit_terms():
    """(hinge shape, edge class, -d theta/d l, edge base minus hinge base) for every edge of every simplex at every hinge."""
    rows = []
    for si, s in enumerate(SHAPES):
        for ti, V, i in INC[s]:
            for e, (a, b) in enumerate(PAIRS):
                rows.append((si, class_index(V[b] - V[a]), -DTH[ti, e], V[a] - V[i]))
    return (np.array([r[0] for r in rows]), np.array([r[1] for r in rows]), np.array([r[2] for r in rows]),
            np.array([r[3] for r in rows], dtype=float))


def _area_terms():
    """(hinge shape, edge class, d A/d l, edge base minus hinge base) for the three edges of each hinge."""
    rows = []
    for si, s in enumerate(SHAPES):
        ti, V, i = INC[s][0]
        a, b, c = TRIS[ti]
        for (u, w), g in zip([(a, b), (b, c), (a, c)], AREA_GRAD[si]):
            rows.append((si, class_index(V[w] - V[u]), g, V[u] - V[i]))
    return (np.array([r[0] for r in rows]), np.array([r[1] for r in rows]), np.array([r[2] for r in rows]),
            np.array([r[3] for r in rows], dtype=float))


_DS, _DC, _DW, _DO = _deficit_terms()
_AS, _AC, _AW, _AO = _area_terms()


def deficit_rows(k):
    """Bloch rows of the linearized deficits: delta'_h(u) = sum_d R_s(k)_d u_d for a hinge of shape s at the origin and the
    plane wave u_e = u_d exp(i k.b_e) on edges of class d based at b_e. delta_h = 2 pi - sum of the dihedral angles at h."""
    R = np.zeros((len(SHAPES), 15), dtype=complex)
    np.add.at(R, (_DS, _DC), _DW * np.exp(1j * (_DO @ np.asarray(k, dtype=float))))
    return R


def area_rows(k):
    """Bloch rows of the hinge areas, which depend only on the hinge's own three edges."""
    R = np.zeros((len(SHAPES), 15), dtype=complex)
    np.add.at(R, (_AS, _AC), _AW * np.exp(1j * (_AO @ np.asarray(k, dtype=float))))
    return R


def D(k):
    """Bloch matrix of the derived weight, per lattice site (Theorem 5.4, dual-cycle counts)."""
    R = deficit_rows(k)
    return 0.5 * (R.conj().T * WEIGHT) @ R


def B(k):
    """Bloch matrix of the flat Regge Hessian from Schlaefli's identity, sum_h A'_h delta'_h^T, per lattice site."""
    return area_rows(k).conj().T @ deficit_rows(k)


def Gamma(k):
    """Vertex-move columns, (Gamma(k) eta)_d = 2 (exp(i k.d) - 1) d.eta (the earlier paper's Section 4.2)."""
    return 2.0 * (np.exp(1j * (CLASSES @ np.asarray(k, dtype=float))) - 1.0)[:, None] * CLASSES


def A(k):
    """The earlier paper's gauge-fixed operator B + Gamma Gamma^dagger."""
    G = Gamma(k)
    return B(k) + G @ G.conj().T


def doubler(c):
    """The doubler v_c of the earlier paper's Proposition 6.1 at the axis corner pi e_c."""
    return np.array([2.0 ** (d.sum() - 3) if d[c] == 0 else 0.0 for d in CLASSES])


def complement(k):
    """Orthonormal basis of the complement of the vertex moves at a wave k != 0."""
    Q, _ = np.linalg.qr(Gamma(k))
    U, _, _ = np.linalg.svd(np.eye(15) - Q @ Q.conj().T)
    return U[:, :11]


def on_complement(k):
    """Eigenvalues of D(k) on the complement of the vertex moves, and the action's softest mode there: the eigenvalue of A of
    least size and the weight's value on its eigenvector. The complement is invariant under A, since A Gamma = Gamma Gamma^dagger
    Gamma."""
    C = complement(k)
    Dc = C.conj().T @ D(k) @ C
    Ac = C.conj().T @ A(k) @ C
    w, V = np.linalg.eigh((Ac + Ac.conj().T) / 2)
    i = int(np.argmin(np.abs(w)))
    u = V[:, i]
    return np.linalg.eigvalsh((Dc + Dc.conj().T) / 2), float(w[i]), float((u.conj() @ Dc @ u).real)


STRAINS = [np.array([d @ H @ d for d in CLASSES]) for H in
           [np.outer(np.eye(4)[a], np.eye(4)[b]) + np.outer(np.eye(4)[b], np.eye(4)[a]) for a in range(4) for b in range(a, 4)]]
WAVES = [np.array([0.3, -0.2, 0.7, 0.1]), np.array([1.0, 2.0, -0.5, 0.25]), np.array([2.9, 0.4, -1.3, 0.8])]
results = []


def report(name, ok, detail):
    results.append(ok)
    print(f"[{'PASS' if ok else 'FAIL'}] {name}: {detail}")


def run_checks(compare):
    # The hinges (Section 5.3).
    areas = sorted({round(float(a), 12) for a in AREA})
    wvals = sorted({round(float(x), 9) for x in WEIGHT})
    wexact = [1 / 12, 1 / 4, SQRT3 / 4]
    report('Hinges', len(SHAPES) == 50 and set(VALENCE) == {4, 6} and VALENCE.sum() == 240
           and np.allclose(areas, [0.5, 1, 1.5, SQRT3]) and len(wvals) == 3
           and all(min(abs(float(w) - e) for e in wexact) < 1e-14 for w in WEIGHT),
           f'{len(SHAPES)} shapes per cell, valences {sorted(set(VALENCE.tolist()))} adding to {VALENCE.sum()} = 24 x 10; areas '
           f'{", ".join(f"{a:.4f}" for a in areas)}; A_h/n_h takes the values {", ".join(f"{w:.6f}" for w in wvals)} = 1/12, 1/4, sqrt3/4')
    worst = max(abs(sum(THETA[ti] for ti, _, _ in INC[s]) - 2 * mp.pi) for s in SHAPES)
    report('Flat lattice', worst < mp.mpf(10) ** -35,
           f'dihedral angles {sorted({int(mp.nint(mp.degrees(t))) for t in THETA})} degrees; around every hinge they add to 2 pi '
           f'(largest deviation {mp.nstr(worst, 3)}), so every flat deficit vanishes')

    # The action's Hessian from Schlaefli's identity.
    herm = max(np.abs(B(k) - B(k).conj().T).max() for k in WAVES)
    bianchi = max(np.abs(B(k) @ Gamma(k)).max() for k in WAVES)
    ev0 = np.linalg.eigvalsh(B(np.zeros(4)).real)
    S = np.array(STRAINS).T
    strain_kernel = np.abs(B(np.zeros(4)) @ S).max()
    corner = []
    for c in range(4):
        k = np.pi * np.eye(4)[c]
        v = doubler(c)
        corner.append((np.abs(B(k) @ v).max(), np.abs(A(k) @ v).max(), np.abs(Gamma(k).conj().T @ v).max(),
                       int(np.sum(np.abs(np.linalg.eigvalsh(A(k))) < 1e-9))))
    report('Action, from Schlaefli', herm < 1e-13 and bianchi < 1e-13 and int(np.sum(np.abs(ev0) < 1e-9)) == 10 and strain_kernel < 1e-13
           and max(max(x[:3]) for x in corner) < 1e-13 and all(x[3] == 1 for x in corner) and abs(doubler(0) @ doubler(0) - 31 / 16) < 1e-15,
           f'B Hermitian ({herm:.1e}), B Gamma = 0 ({bianchi:.1e}); kernel of B(0) = the ten constant strains; at each axis corner '
           f'B v_c = A v_c = Gamma^dagger v_c = 0 and A has a one-dimensional kernel; |v_c|^2 = 31/16 '
           f'(earlier paper, Section 4.2 and Proposition 6.1)')
    if compare:
        sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
        import kuhn_bloch
        diff = max(np.abs(B(k) - kuhn_bloch.B(k)).max() for k in WAVES + [np.pi * np.eye(4)[0]])
        report('Action, against Appendix B', diff < 1e-13,
               f'largest difference from the Hessian assembled by kuhn_bloch.py from the printed cell Hessian: {diff:.1e}')

    # Theorem 5.3.
    flat = np.abs(deficit_rows(np.zeros(4)) @ S).max()
    moves = max(np.abs(deficit_rows(k) @ Gamma(k)).max() for k in WAVES)
    census = []
    for c in range(4):
        R = deficit_rows(np.pi * np.eye(4)[c])
        v4 = 4 * doubler(c)
        tot = sum((1 + (-1) ** s[0][c] + (-1) ** (s[0][c] + s[1][c])) * (R[si] @ v4) for si, s in enumerate(SHAPES))
        census.append(tot.real)
    exact = 6 + 8 * SQRT3
    report('Theorem 5.3', flat < 1e-13 and moves < 1e-13 and max(abs(x - exact) for x in census) < 1e-12,
           f'(i) linearized deficits vanish on constant strains ({flat:.1e}) and vertex moves ({moves:.1e}); (ii) the signed census '
           f'of (-1)^x_c 4 v_c is N^4 times {", ".join(f"{x:.12f}" for x in census)} for c = 0..3; 6 + 8 sqrt3 = {exact:.12f}')

    # Theorem 5.5 and the numbers after it.
    d0 = np.linalg.eigvalsh(D(np.zeros(4)).real)
    on_strains = np.abs(D(np.zeros(4)) @ S).max()
    on_moves = max(np.abs(D(k) @ Gamma(k)).max() for k in WAVES)
    report('Theorem 5.5 (i)', on_strains < 1e-13 and on_moves < 1e-13 and int(np.sum(np.abs(d0) < 1e-9)) == 10,
           f'D vanishes on the constant strains ({on_strains:.1e}) and the vertex moves ({on_moves:.1e}); D(0) has exactly ten zero '
           f'eigenvalues')
    lines, ok4, low = [], True, []
    for name, e in [('axis', np.array([1.0, 0, 0, 0])), ('face diagonal', np.array([1.0, 1, 0, 0]) / math.sqrt(2)),
                    ('(1, 2, -1, 0.5)', np.array([1.0, 2, -1, 0.5]) / 2.5)]:
        r = []
        for q in (0.1, 0.05, 0.025):
            ev, _, _ = on_complement(q * e)
            r.append((ev[0] / q ** 4, ev[5] / q ** 4, ev[6]))
        ok4 &= all(abs(r[2][j] / r[1][j] - 1) < 0.01 for j in (0, 1))
        low.append(min(x[2] for x in r))
        lines.append(f'{name}: ' + ', '.join(f'|k| = {q}: {a:.5f} to {b:.5f}' for q, (a, b, _) in zip((0.1, 0.05, 0.025), r)))
    report('Theorem 5.5 (ii)', ok4, 'the six metric eigenvalues of D(k) on the complement of the vertex moves, divided by |k|^4, '
           'settle as |k| -> 0; ' + '; '.join(lines))
    gap = d0[10:]
    report('Non-metric modes', abs(gap[1] - 65 / 9) < 1e-12 and abs(gap[3] - 65 / 9) < 1e-12 and min(low) > 4.9,
           f'at k -> 0 the five non-metric eigenvalues of D are {", ".join(f"{x:.6f}" for x in gap)} (three equal to 65/9); '
           f'the lowest of them lies between {min(low):.4f} and {gap[0]:.4f} at |k| = 0.1, 0.05 and 0.025 on the three lines')
    damping, inverse = [], []
    for c in range(4):
        k = np.pi * np.eye(4)[c]
        v = doubler(c)
        damping.append(float((v @ D(k) @ v).real / (v @ v)))
        inverse.append((np.linalg.svd(D(k) - 1j * A(k), compute_uv=False).min(), np.linalg.svd(A(k), compute_uv=False).min()))
    exact = (10 + 3 * SQRT3) / 124
    report('Theorem 5.5 (iii)', min(damping) > 0 and max(abs(x - exact) for x in damping) < 1e-13 and min(x[0] for x in inverse) > 1e-3,
           f'v_c^dagger D v_c / |v_c|^2 = {", ".join(f"{x:.13f}" for x in damping)} at the four axis corners; (10 + 3 sqrt3)/124 = '
           f'{exact:.13f}; smallest singular value '
           f'of D - i A there {min(x[0] for x in inverse):.4f}, of A alone {max(x[1] for x in inverse):.1e} (Theorem 5.2, a = 1)')
    ranks = [int(np.linalg.matrix_rank(deficit_rows(np.pi * np.eye(4)[c]), tol=1e-9)) for c in range(4)]
    ranks += [int(np.linalg.matrix_rank(deficit_rows(k), tol=1e-9)) for k in WAVES]
    report('After Theorem 5.6', ranks == [11] * len(ranks),
           f'rank of the linearized deficits {ranks} at the four axis corners and {len(WAVES)} sampled waves, so their kernel is the '
           f'four vertex moves')


def figure_data(path, n=200):
    out = {'program': 'curvature_weight.py', 'paper': 'https://doi.org/10.5281/zenodo.23813297',
           'line': 'k = s pi e_0, s from 0 to 1 (the axis corner)', 'counts': 'dual cycle, n_h = valence',
           's': [], 'weight_eigenvalues': [], 'action_softest': [], 'weight_on_softest': []}
    for j in range(1, n + 1):
        s = j / n
        ev, a, w = on_complement(np.array([s * math.pi, 0, 0, 0]))
        out['s'].append(s)
        out['weight_eigenvalues'].append([float(f'{x:.6e}') for x in ev])
        out['action_softest'].append(float(f'{a:.6e}'))
        out['weight_on_softest'].append(float(f'{w:.6e}'))
    v = doubler(0)
    out['doubler_damping'] = float(f'{(v @ D(np.pi * np.eye(4)[0]) @ v).real / (v @ v):.7f}')
    with open(path, 'w') as f:
        json.dump(out, f, separators=(',', ':'))
    print(f'wrote {path}: {n} points on the axis line')


def main():
    parser = argparse.ArgumentParser(description=__doc__.splitlines()[0])
    parser.add_argument('--json', metavar='FILE', help='write the figure data')
    parser.add_argument('--compare', action='store_true', help='compare the Hessian with kuhn_bloch.py (same folder)')
    args = parser.parse_args()
    run_checks(args.compare)
    if args.json:
        figure_data(args.json)
    print(f'\n{sum(results)} of {len(results)} checks passed. The values are numerical evaluations of the exact formulas; the weight '
          'rests on the identifications (I1) and (I2), the dual-cycle channel count and the scaling of the phase.')
    if not all(results):
        sys.exit(1)


if __name__ == '__main__':
    main()
