"""Example 43 — Prime-ladder Hamiltonian (TNFR-Riemann P14).

Demonstrates the operational closure of gap G1 in the TNFR-Riemann
research programme:

* construct the canonical TNFR internal Hamiltonian on the prime-ladder
  graph,
* certify that its spectrum reproduces the prime-ladder spectrum
  ``{k * log(p)}`` exactly,
* certify that its weighted spectral trace
  ``Tr(W exp(-s H))`` reproduces the von Mangoldt weighted Dirichlet
  trace ``-zeta'(s)/zeta(s)`` (in the convergent regime),
* verify the canonical TNFR property that the time-evolution operator
  ``U(t) = exp(-i t H / hbar_str)`` is unitary (self-adjointness check
  at the dynamics level),
* probe how a small inter-node ladder coupling perturbs the spectrum.

Run with ``PYTHONPATH=./src`` so the local tnfr package is used
instead of any stale installed copy.
"""

from __future__ import annotations

import math
import sys

import numpy as np

from tnfr.riemann import (
    build_prime_ladder_hamiltonian,
    tnfr_log_zeta_derivative,
    verify_hamiltonian_reproduces_prime_ladder,
    weighted_spectral_trace,
)


def _safe_utf8_stdout() -> None:
    try:
        sys.stdout.reconfigure(encoding="utf-8")  # type: ignore[attr-defined]
    except Exception:
        pass


def section(title: str) -> None:
    bar = "=" * 78
    print(f"\n{bar}\n {title}\n{bar}")


def main() -> int:
    _safe_utf8_stdout()

    n_primes = 12
    max_power = 6

    # ------------------------------------------------------------------
    section("1. Construct the prime-ladder Hamiltonian (decoupled, J0=0)")
    # ------------------------------------------------------------------
    bundle = build_prime_ladder_hamiltonian(
        n_primes,
        max_power=max_power,
        coupling=0.0,
    )
    H = bundle.hamiltonian
    W = bundle.weight_operator
    print(f"  number of primes      : {n_primes}")
    print(f"  REMESH echo cap K     : {max_power}")
    print(f"  Hilbert dimension N   : {H.N}")
    print(f"  H is Hermitian        : True (verified at construction)")
    print(f"  ladder coupling J0    : {bundle.coupling}")
    print(f"  primes used           : {list(bundle.spectrum.primes)}")

    # ------------------------------------------------------------------
    section("2. Spectrum reproduction certificate (gap G1, exact at J0=0)")
    # ------------------------------------------------------------------
    cert = verify_hamiltonian_reproduces_prime_ladder(
        bundle,
        s_values=(2.0, 3.0, 5.0, 10.0),
    )
    print(f"  spectrum max abs error: {cert.spectrum_max_abs_error:.3e}")
    print(f"  spectrum reproduced   : {cert.spectrum_reproduced}")
    print()
    print("  Sample sorted eigenvalues (Hamiltonian vs ladder reference):")
    eigvals_ham, _ = H.get_spectrum()
    eigvals_ref = np.sort(bundle.spectrum.eigenvalues)
    eigvals_ham_sorted = np.sort(np.real(eigvals_ham))
    for i in range(0, len(eigvals_ref), max(1, len(eigvals_ref) // 6)):
        print(
            f"    n={i:3d}  E_Ham={eigvals_ham_sorted[i]:.10f}  "
            f"E_ref={eigvals_ref[i]:.10f}  "
            f"|diff|={abs(eigvals_ham_sorted[i] - eigvals_ref[i]):.2e}"
        )

    # ------------------------------------------------------------------
    section("3. Weighted spectral trace vs classical -zeta'(s)/zeta(s)")
    # ------------------------------------------------------------------
    print(f"  s values              : {list(cert.s_values)}")
    print(f"  max relative error    : {cert.trace_max_rel_error:.3e}")
    print(f"  trace reproduced      : {cert.trace_reproduced}")
    print(f"  overall certificate   : {cert.overall_ok}")
    print()
    print("  Per-s comparison Z_H(s) vs Z_vM(s):")
    for s in cert.s_values:
        z_ham = weighted_spectral_trace(H.H_int, W, float(s))
        z_ref = float(tnfr_log_zeta_derivative(bundle.spectrum, float(s)))
        print(
            f"    s={s:5.2f}  Z_H={z_ham:.12f}  Z_vM={z_ref:.12f}  "
            f"|rel|={abs(z_ham - z_ref) / abs(z_ref):.2e}"
        )

    # ------------------------------------------------------------------
    section("4. Unitarity of time evolution U(t) = exp(-i t H / hbar_str)")
    # ------------------------------------------------------------------
    for t in (0.1, 1.0, 5.0):
        U = H.time_evolution_operator(t)
        U_dag = U.conj().T
        I = np.eye(H.N, dtype=complex)
        unit_err = float(np.max(np.abs(U_dag @ U - I)))
        print(f"    t={t:5.2f}   ||U^H U - I||_inf = {unit_err:.3e}")

    # ------------------------------------------------------------------
    section("5. Perturbative probe: small ladder coupling J0 > 0")
    # ------------------------------------------------------------------
    print("  Spectrum deviation from prime ladder as J0 increases:")
    print("  (departure is EXPECTED — Euler product only holds at J0=0)")
    print()
    print("    J0        |E_max_dev|   |trace_rel_dev| @ s=2")
    for j0 in (0.0, 1e-4, 1e-3, 1e-2, 5e-2):
        b = build_prime_ladder_hamiltonian(
            n_primes,
            max_power=max_power,
            coupling=j0,
        )
        ev, _ = b.hamiltonian.get_spectrum()
        spec_dev = float(
            np.max(np.abs(np.sort(np.real(ev)) - np.sort(b.spectrum.eigenvalues)))
        )
        z_ham = weighted_spectral_trace(b.hamiltonian.H_int, b.weight_operator, 2.0)
        z_ref = float(tnfr_log_zeta_derivative(b.spectrum, 2.0))
        rel_dev = abs(z_ham - z_ref) / abs(z_ref)
        print(f"    {j0:9.1e}  {spec_dev:.3e}    {rel_dev:.3e}")

    section("Conclusion")
    print(
        "  The canonical TNFR internal Hamiltonian, instantiated on the\n"
        "  prime-ladder graph with J0 = 0, is a finite-dimensional\n"
        "  self-adjoint operator whose spectrum exactly equals the\n"
        "  prime-ladder spectrum {k log p} and whose weighted spectral\n"
        "  trace equals -zeta'(s)/zeta(s) in the convergent regime.\n"
        "  This is the operational closure of gap G1 in the TNFR-Riemann\n"
        "  programme (the existence of an explicit self-adjoint operator\n"
        "  carrying the von Mangoldt spectral data).  Closure of the\n"
        "  remaining gaps G3-G4 (localising zeros on Re(s) = 1/2) is the\n"
        "  substance of RH itself and is not addressed by this construction."
    )
    return 0


if __name__ == "__main__":
    raise SystemExit(main())