r"""Richardson-style convergence study for the cubic Taylor coefficient
of the nonlinear Regge action on a 3x3x3 periodic Freudenthal torus.

Outputs a table of d3 at multiple step sizes and a Richardson extrapolation,
plus sanity checks at directions where the cubic must vanish.

Seed is fixed for reproducibility. The probe direction xi is drawn once with
np.random.seed(SEED), scaled by SCALE, and reused across all step sizes.

Conventions:
- Vertex potentials live in R^27 with index v_idx(x, y, z) = x*9 + y*3 + z, x,y,z in {0,1,2}.
- Conformal ansatz: L_ij(xi) = L0_ij * exp((xi_i + xi_j) / 2).
- Base squared edge length L0_ij^2 = shortest periodic distance squared in {0,1,2}^3.
- Deficit angles: positive ε_e means concentrated positive curvature at hinge e.
- Regge action S(xi) = sum_e L_e(xi) * (2*pi - sum_{tau \supset e} theta_{e,tau}(xi)).
"""

import numpy as np

SEED = 42
SCALE = 0.1

N = 3
V = N ** 3


def v_idx(x, y, z):
    return (x % N) * N * N + (y % N) * N + (z % N)


# Six-tet decomposition of each unit cube (Freudenthal subdivision)
perms = [
    ((0, 0, 0), (1, 0, 0), (1, 1, 0), (1, 1, 1)),
    ((0, 0, 0), (1, 0, 0), (1, 0, 1), (1, 1, 1)),
    ((0, 0, 0), (0, 1, 0), (1, 1, 0), (1, 1, 1)),
    ((0, 0, 0), (0, 1, 0), (0, 1, 1), (1, 1, 1)),
    ((0, 0, 0), (0, 0, 1), (1, 0, 1), (1, 1, 1)),
    ((0, 0, 0), (0, 0, 1), (0, 1, 1), (1, 1, 1)),
]

tets = []
for x in range(N):
    for y in range(N):
        for z in range(N):
            for p in perms:
                tet = tuple(v_idx(x + dx, y + dy, z + dz) for dx, dy, dz in p)
                tets.append(tet)


def decode(idx):
    x = idx // (N * N)
    y = (idx // N) % N
    z = idx % N
    return x, y, z


def base_sq_len(u, v):
    ux, uy, uz = decode(u)
    vx, vy, vz = decode(v)
    dx = min((vx - ux) % N, (ux - vx) % N)
    dy = min((vy - uy) % N, (uy - vy) % N)
    dz = min((vz - uz) % N, (uz - vz) % N)
    return dx ** 2 + dy ** 2 + dz ** 2


def cayley_menger_matrix(sq_lengths):
    M = np.ones((5, 5))
    M[0, 0] = 0
    M[1:, 1:] = sq_lengths
    return M


def dihedral_angle(sq_lengths, i, j):
    others = [x for x in range(4) if x != i and x != j]
    k, l = others
    M = cayley_menger_matrix(sq_lengths)

    def cofactor(r, c):
        subM = np.delete(np.delete(M, r, axis=0), c, axis=1)
        return ((-1) ** (r + c)) * np.linalg.det(subM)

    C_kk = cofactor(k + 1, k + 1)
    C_ll = cofactor(l + 1, l + 1)
    C_kl = cofactor(k + 1, l + 1)

    cos_theta = C_kl / np.sqrt(C_kk * C_ll)
    cos_theta = np.clip(cos_theta, -1.0, 1.0)
    return np.arccos(cos_theta)


def regge_action(xi):
    edge_deficits = {}
    edge_lengths = {}
    for tet in tets:
        sq_lengths = np.zeros((4, 4))
        for i in range(4):
            for j in range(4):
                if i != j:
                    u, v = tet[i], tet[j]
                    L0_sq = base_sq_len(u, v)
                    sq_lengths[i, j] = L0_sq * np.exp(xi[u] + xi[v])
        for i in range(4):
            for j in range(i + 1, 4):
                u, v = tet[i], tet[j]
                edge = tuple(sorted((u, v)))
                theta = dihedral_angle(sq_lengths, i, j)
                if edge not in edge_deficits:
                    edge_deficits[edge] = 2 * np.pi
                    edge_lengths[edge] = np.sqrt(sq_lengths[i, j])
                edge_deficits[edge] -= theta
    action = 0.0
    for edge, deficit in edge_deficits.items():
        action += edge_lengths[edge] * deficit
    return action


def j_cost_action(xi):
    edge_lengths = {}
    for tet in tets:
        for i in range(4):
            for j in range(i + 1, 4):
                u, v = tet[i], tet[j]
                edge = tuple(sorted((u, v)))
                if edge not in edge_lengths:
                    edge_lengths[edge] = np.sqrt(base_sq_len(u, v))
    J = 0.0
    for edge, L0 in edge_lengths.items():
        u, v = edge
        J += L0 * (np.cosh(xi[u] - xi[v]) - 1)
    return J


def central_fd_cubic(f, eta, h):
    """5-point central FD estimate of f'''(0; eta) with O(h^2) truncation error."""
    Sp1 = f(h * eta)
    Sm1 = f(-h * eta)
    Sp2 = f(2 * h * eta)
    Sm2 = f(-2 * h * eta)
    return (Sp2 - 2 * Sp1 + 2 * Sm1 - Sm2) / (2 * h ** 3)


def richardson_extrap(d_hi, d_lo):
    """Order-2 Richardson: hi at step h, lo at step h/10.
    If truncation error is O(h^2), the extrapolation is
    (100 * d_lo - d_hi) / 99."""
    return (100.0 * d_lo - d_hi) / 99.0


def run_table(f, eta, label):
    print(f"\n==== {label} ====")
    print(f"{'h':>12} {'cubic estimate':>22}")
    estimates = {}
    for h in [1e-2, 1e-3, 1e-4, 1e-5]:
        d3 = central_fd_cubic(f, eta, h)
        estimates[h] = d3
        print(f"{h:12.0e} {d3:22.10f}")
    print()
    print("Richardson extrapolations (assuming O(h^2) truncation):")
    pairs = [(1e-2, 1e-3), (1e-3, 1e-4), (1e-4, 1e-5)]
    for h_hi, h_lo in pairs:
        rex = richardson_extrap(estimates[h_hi], estimates[h_lo])
        print(f"  h={h_hi:.0e} -> h={h_lo:.0e}: {rex:22.10f}")
    return estimates


if __name__ == "__main__":
    # Probe direction: legacy NumPy global RNG, seed 42, np.random.randn(V) * 0.1.
    # We use the legacy global-RNG form so the probe is bitwise reproducible
    # from the printed two-line "np.random.seed(42); eta = np.random.randn(27) * 0.1"
    # in the manuscript.
    np.random.seed(SEED)
    eta_random = np.random.randn(V) * SCALE
    print("Probe direction eta (np.random.seed=42, eta = np.random.randn(27) * 0.1):")
    for i, x in enumerate(eta_random):
        print(f"  eta[{i:2d}] = {x:+.10f}")

    run_table(regge_action, eta_random, "Regge cubic along seed-42 direction")
    run_table(j_cost_action, eta_random, "J-cost cubic along seed-42 direction (must be ~0)")

    # Sanity check: along the all-ones direction, conformal scale invariance
    # of dihedral angles together with flat background gives Regge(t * 1) = 0
    # identically, so every derivative is zero.
    eta_constant = np.ones(V) * SCALE
    run_table(regge_action, eta_constant, "Regge cubic along constant direction (must be ~0)")
    run_table(j_cost_action, eta_constant, "J-cost cubic along constant direction (must be 0 identically)")