#!/usr/bin/env python3
"""
N3 — THE FEE SCHEDULE, MEASURED (the thermometer that worked)
===============================================================
Monogamy Gravity, Folio IV numerical experiment.

CLAIM UNDER TEST (Postulate 4, "The Fee"):
    A causal region of size L bills at the diamond temperature:
    T(L) ~ 1/L, with centre constant T*L = 2/pi (CHM / Hislop-Longo).

GATE (written before the run):
    If the temperature extracted from the vacuum's own entanglement
    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.

LAB NOTEBOOK (kept honest):
    Attempt 1 — harmonic chain (n3_fee.py): FAILED. The bosonic zero
    mode obstructs the lattice EH -> Bisognano-Wichmann correspondence
    (wedge calibration gives slope 0.05-0.2 x theorem). Same culprit as
    the Folio V critical anomaly. Filed.
    Attempt 2 — free fermions in float64: FAILED ABOVE L~25. The EH is
    h = ln((1-C)/C); once eigenvalues of C reach machine epsilon the log
    saturates (~32) and fabricates a temperature plateau. A numerical
    lie, caught by the L=21 row disagreeing with the trend.
    Attempt 3 — free fermions in arbitrary precision (this file): PASSED.

METHOD:
    Infinite hopping chain at half filling (v_F = 1), exact sine-kernel
    correlations C_ij = sin(pi(i-j)/2)/(pi(i-j)). For an interval of L
    sites, the entanglement Hamiltonian is h = ln((1-C)/C), computed
    with mpmath at 30-60 significant digits. CHM predicts the
    nearest-neighbour EH amplitude profile |t_i| ~ pi * f(i),
    f(x) = (x-a)(b-x)/(b-a); the centre inverse temperature is
    beta_c = 2*pi*f_max * coef. Fit coef centre-weighted, sweep L.
"""

import numpy as np
import json
import time
from mpmath import mp, mpf, sin, pi as mpi, log, matrix, eigsy


def eh_offdiagonal(L, dps):
    """Nearest-neighbour amplitudes of the EH of L sites (sine kernel)."""
    mp.dps = dps
    C = matrix(L, L)
    for i in range(L):
        for j in range(L):
            d = i - j
            C[i, j] = mpf('0.5') if d == 0 else sin(mpi * d / 2) / (mpi * d)
    E, V = eigsy(C)
    g = [log((1 - E[k]) / E[k]) for k in range(L)]
    out = []
    for i in range(L - 1):
        s = mpf(0)
        for k in range(L):
            s += V[i, k] * g[k] * V[i + 1, k]
        out.append(float(s))
    return np.array(out)


if __name__ == '__main__':
    print("N3 — the fee schedule, measured (fermionic thermometer, mpmath)")
    print(f"{'L':>4} {'dps':>5} {'beta_c':>10} {'T*L':>9} {'corr':>8} {'time':>6}")
    rows = []
    for L, dps in ((11, 30), (15, 35), (21, 40), (27, 50), (35, 60)):
        t0 = time.time()
        t = np.abs(eh_offdiagonal(L, dps))
        xs = np.arange(L - 1) + 0.5
        f = xs * (L - 1 - xs) / (L - 1)
        wgt = np.exp(-((xs - (L - 1) / 2) / (L / 4)) ** 2)
        coef = float(np.sum(wgt * t * np.pi * f) / np.sum(wgt * (np.pi * f) ** 2))
        corr = float(np.corrcoef(t, f)[0, 1])
        beta = coef * 2 * np.pi * (L - 1) / 4
        rows.append(dict(L=L, beta=beta, TL=L / beta, corr=corr))
        print(f"{L:4d} {dps:5d} {beta:10.3f} {L / beta:9.4f} {corr:8.4f} "
              f"{time.time() - t0:5.0f}s", flush=True)

    LL = np.array([r['L'] for r in rows], float)
    TT = 1 / np.array([r['beta'] for r 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: (2/pi)/L = 0.6366/L.  "
          f"Largest-L constant T*L = {rows[-1]['TL']:.4f} vs 2/pi = {2 / np.pi:.4f}")
    gate = abs(pw + 1.0) < 0.15
    print("GATE: " + ("PASSED — the vacuum bills at 1/L." if gate
                      else "TRIPPED — the invoice is mispriced."))
    json.dump(dict(rows=rows, power=float(pw), A=float(np.exp(lnA))),
              open('/Users/antoine/agi/ledger/n3_fermion_hp.json', 'w'), indent=1)
    print("-> n3_fermion_hp.json")