#!/usr/bin/env python3
"""
N3 — THE FEE SCHEDULE, MEASURED
================================
Monogamy Gravity, Folio IV numerical experiment.

CLAIM UNDER TEST (Postulate 4, "The Fee"):
    Every causal region of size L bills transactions at the diamond
    temperature T_L = hbar*c / (2*pi*k_B*L) -- up to the geometric O(1)
    profile factor. CHM (Casini-Huerta-Myers / Hislop-Longo) make this
    exact for an interval: the entanglement Hamiltonian is
        K = 2*pi * Integral f(x) T_00(x) dx,  f(x) = (x-a)(b-x)/(b-a),
    i.e. a LOCAL inverse temperature beta(x) = 2*pi*f(x). At the centre
    f = L/4, so the centre temperature is
        T_centre(L) = 2/(pi*L)        [c = k_B = hbar = 1]
    => the fee schedule: T * L = 2/pi ~ 0.6366, distance-law 1/L.

GATE (written before the run):
    If the temperature extracted from the lattice's own modular
    Hamiltonians does not fall as 1/L (power within ~15% of -1), the fee
    schedule of Postulate 4 has the wrong distance law and the invoice
    of Folio I is mispriced.

METHOD:
    Exact single-particle entanglement Hamiltonian of an interval
    (Peschel construction for Gaussian states with <qp> = 0):
        M = X^{1/2} P X^{1/2},  eigendecomposition M = V nu^2 V^T,
        h_pp = X^{1/2} V [ ln((nu+1/2)/(nu-1/2)) / (2 nu) ] V^T X^{1/2} * ...
    We use the standard form: the p-quadratic part of K has weight
    profile ~ 2*pi*f(x) on the diagonal. We extract the profile,
    fit the parabola, and read the centre temperature -- for a family
    of interval sizes L. Then fit T(L) ~ A / L^p.

Plain numpy.
"""

import numpy as np
import json
import warnings
warnings.filterwarnings("ignore", category=RuntimeWarning)

N, m = 400, 1e-3
Ls = (21, 41, 61, 81, 101, 121)

# ---------------- vacuum of the chain
K0 = np.zeros((N, N))
for i in range(N):
    K0[i, i] = m * m + 2.0
    K0[i, (i + 1) % N] = -1.0
    K0[(i + 1) % N, i] = -1.0
w2, U = np.linalg.eigh(K0)
w = np.sqrt(np.clip(w2, 1e-30, None))
X = (U * (0.5 / w)) @ U.T
P = (U * (0.5 * w)) @ U.T

def sqrtm_sym(A):
    e, V = np.linalg.eigh(A)
    return (V * np.sqrt(np.clip(e, 1e-30, None))) @ V.T

def modular_hpp(region):
    """Single-particle entanglement Hamiltonian, p-quadratic block.
    K = 1/2 sum p h_pp p + 1/2 sum q h_qq q  (for <qp>=0 states).
    h_pp = P^{-1/2} g(sqrt(P^{1/2} X P^{1/2})) P^{-1/2} with
    g(nu) = nu * ln((nu+1/2)/(nu-1/2))  -- standard Peschel form,
    roles of X and P swapped relative to h_qq."""
    idx = np.ix_(region, region)
    Xr, Pr = X[idx], P[idx]
    Ph = sqrtm_sym(Pr)
    e, V = np.linalg.eigh(Ph @ Xr @ Ph)
    nu = np.sqrt(np.clip(e, 0.25 + 1e-12, None))
    g = nu * np.log((nu + 0.5) / (nu - 0.5))
    Phi = np.linalg.inv(Ph)
    return Phi @ ((V * g) @ V.T) @ Phi

rows = []
print(f"N3 — the fee schedule, measured. chain N={N}, m={m}")
print(f"{'L':>5} {'beta_centre':>12} {'T_centre':>11} {'T*L':>9} {'profile corr':>13}")
for L in Ls:
    a = (N - L) // 2
    region = list(range(a, a + L))
    hpp = modular_hpp(region)
    prof = np.diag(hpp)
    xs = np.arange(L, dtype=float)
    f = (xs) * (L - 1 - xs) / (L - 1)            # CHM parabola on the lattice
    # fit prof = coef * 2*pi*f  (centre-weighted to dodge boundary artifacts)
    wgt = np.exp(-((xs - (L - 1) / 2) / (L / 4)) ** 2)
    coef = float(np.sum(wgt * prof * 2 * np.pi * f) /
                 np.sum(wgt * (2 * np.pi * f) ** 2))
    corr = float(np.corrcoef(prof, f)[0, 1])
    beta_c = coef * 2 * np.pi * (L - 1) / 4       # beta at centre
    T_c = 1.0 / beta_c
    rows.append(dict(L=L, beta=beta_c, T=T_c, TL=T_c * L, corr=corr))
    print(f"{L:5d} {beta_c:12.3f} {T_c:11.5f} {T_c * L:9.4f} {corr:13.4f}")

# ---------------- the fee schedule fit
LL = np.array([q['L'] for q in rows], float)
TT = np.array([q['T'] for q in rows])
pw, lnA = np.polyfit(np.log(LL), np.log(TT), 1)
print(f"\nfee schedule: T(L) = {np.exp(lnA):.4f} * L^{pw:+.3f}")
print(f"CHM prediction: T(L) = (2/pi)/L = 0.6366 * L^-1")
print(f"measured constant T*L (largest L): {rows[-1]['TL']:.4f}  vs 2/pi = {2/np.pi:.4f}")
gate = abs(pw + 1.0) < 0.15
print("\nGATE: " + ("PASSED — the vacuum bills at 1/L." if gate else
                    "TRIPPED — wrong distance law; the invoice is mispriced."))
json.dump(dict(N=N, m=m, rows=rows, power=float(pw), A=float(np.exp(lnA))),
          open('/Users/antoine/agi/ledger/n3_results.json', 'w'), indent=1)
print("-> n3_results.json")