nodal_propagator_residue_bridge.py

"""Camino 17 - Nodal-propagator residue bridge.

Canonical rebuild of the §13nonies.4 / P31 oscillatory-half (S(T))
reconstruction experiment, redone "from TNFR nodal structure and
dynamics in the most canonical way possible" using the post-N15
machinery that did not exist when P31 was written.

What changed relative to P31 (`oscillatory_correction.py`)
----------------------------------------------------------
P31 built S_TNFR by reading off the classical Riemann-Siegel template
(the sum -(1/pi) sum (1/k) p^{-k/2} sin(T k log p)) and *plugging* the
canonical prime-ladder spectrum into it.  The docstring of P31 admits
the open gap: expressing S_TNFR as a **derivation** of the canonical
nodal evolution (rather than as an ingredient injected into Riemann's
template) was still open.

Camino 17 closes that *methodological* gap (not the mathematical one):
the oscillatory observable here is **generated** by the canonical
structural propagator e^{-s H_P14}, i.e. by the nodal time-evolution
operator itself, through the canonical weighted spectral trace

    Z(s) = Tr(W e^{-s H_freq}) = sum_{p,k} log(p) e^{-s k log p}.

On the critical line s = 1/2 + iT the imaginary part

    Im Z(1/2 + iT) = - sum_{p,k} log(p) p^{-k/2} sin(T k log p)

is the von Mangoldt oscillation *emitted* by the propagator, with no
external template.  We then diagnose it with the two strongest
post-N15 structural results:

  * R-infinity (N15, Branch A): the canonical REMESH-infinity
    orthogonal projector splits any signal into range (smooth /
    resonant) + kernel (oscillatory).  The propagator oscillation
    lands in the kernel.
  * CCET / S_n-equivariance: the canonical observable is a symmetric
    function of the spectrum; relabelling the primes leaves the kernel
    residue invariant to machine precision -> the residue is
    S_n-degenerate and cannot carry the prime-specific arithmetic
    correlations that the true S(T) requires.

HONEST SCOPE (non-negotiable)
-----------------------------
This harness closes NOTHING.  G4 = RH stays OPEN.  R+ and pi remain
the assumed substrate.  The result is a *sharper, more
canonical NEGATIVE result*: it exhibits, dynamically, WHY the canonical
13-operator machinery cannot reach S(T) = (1/pi) arg zeta(1/2 + iT) --
the propagator output sits in ker(R-infinity) AND in Fix(S_n), while
S(T) requires the Fix(S_n)-perpendicular part of the kernel.  This
strengthens branch-B2 evidence (a genuinely new canonical operator,
not derivable from the present catalog, would be required).  It does
NOT prove B2, B3, or RH.

Tests (all four PASS by *confirming the wall*, not by breaking it):
  T1  Canonical propagator reproduces von Mangoldt (sanity).
  T2  Propagator oscillation Im Z(1/2+iT) lands in ker(R-infinity).
  T3  The kernel residue is S_n-degenerate (CCET, dynamical).
  T4  Amplitude undercount vs classical S(T) + R-infinity certificate.
"""

from __future__ import annotations

import math
import os
import sys

import numpy as np

# --- dual sys.path: benchmarks dir + ../src ------------------------------
_HERE = os.path.dirname(os.path.abspath(__file__))
_SRC = os.path.normpath(os.path.join(_HERE, "..", "src"))
for _p in (_HERE, _SRC):
    if _p not in sys.path:
        sys.path.insert(0, _p)

# --- guarded canonical imports -------------------------------------------
try:
    from tnfr.riemann.prime_ladder_hamiltonian import (
        build_prime_ladder_hamiltonian,
        weighted_spectral_trace,
    )
    from tnfr.riemann.remesh_infinity_residue_split import (
        compute_residue_split_certificate,
        split_residue_by_remesh_infinity,
    )
    from tnfr.riemann.von_mangoldt import build_prime_ladder_spectrum

    _HAVE_TNFR = True
except Exception as exc:  # pragma: no cover - import guard
    _HAVE_TNFR = False
    _IMPORT_ERROR = exc

try:
    import mpmath as mp

    _HAVE_MPMATH = True
except Exception:  # pragma: no cover - optional dependency
    _HAVE_MPMATH = False


# --- canonical constants -------------------------------------------------
TAU_L = 4  # local REMESH period
TAU_G = 8  # global REMESH period; L = lcm(4, 8) = 8
N_PRIMES = 100  # number of primes in the ladder
MAX_POWER = 8  # REMESH echo cap K (matches P31)
N_SAMPLES = 2048  # divisible by lcm(TAU_L, TAU_G) = 8
T_MIN = 0.0
T_MAX = 160.0
THRESH = 0.05  # R-infinity 5% range/kernel decision threshold
SEED = 12345  # prime-shuffle seed (T3)

# documented fallbacks (used only if mpmath is unavailable)
_REF_LOG_ZETA_DERIV_2 = 0.5699664469153647  # -zeta'(2)/zeta(2)
_CLASSICAL_S_FALLBACK = 0.5  # conservative max|S(T)|


# --- canonical observable builders ---------------------------------------
def _propagator_oscillation(
    eigenvalues: np.ndarray,
    weights: np.ndarray,
    t_grid: np.ndarray,
    sigma: float = 0.5,
) -> np.ndarray:
    """Im Tr(W e^{-s H_freq}) on Re(s)=sigma, generated by the canonical
    structural propagator e^{-iT H_P14}.

    Vectorised over the T-grid.  For the decoupled (diagonal) P14
    Hamiltonian this equals ``weighted_spectral_trace`` exactly; the
    equivalence is asserted numerically in T1.
    """
    mu = np.asarray(eigenvalues, dtype=float)
    w = np.asarray(weights, dtype=float)
    s = sigma + 1j * np.asarray(t_grid, dtype=float)
    # Z(s) = sum_i w_i exp(-s mu_i)
    z = np.exp(-np.outer(s, mu)) @ w
    return np.asarray(z.imag, dtype=float)


def _s_tnfr_oscillation(
    eigenvalues: np.ndarray,
    weights: np.ndarray,
    t_grid: np.ndarray,
    sigma: float = 0.5,
) -> np.ndarray:
    """P31 oscillatory term built from the *canonical* spectrum:

        S_TNFR(T) = -(1/pi) sum_{p,k} (1/k) p^{-k sigma} sin(T k log p).

    Here k = mu / w (since mu = k log p, w = log p) and
    p^{-k sigma} = e^{-sigma mu}.
    """
    mu = np.asarray(eigenvalues, dtype=float)
    w = np.asarray(weights, dtype=float)
    k = mu / w  # echo index
    amp = (1.0 / k) * np.exp(-sigma * mu)
    sins = np.sin(np.outer(np.asarray(t_grid, float), mu))
    return -(1.0 / math.pi) * (sins @ amp)


def _ref_log_zeta_deriv_2() -> float:
    if _HAVE_MPMATH:
        mp.mp.dps = 25
        return float(-mp.zeta(2, 1, 1) / mp.zeta(2))
    return _REF_LOG_ZETA_DERIV_2


def _classical_max_abs_s() -> tuple[float, str]:
    """Max |S(T)| over a window, from mpmath if available.

    S(T) = N(T) - 1 - theta(T)/pi, with N(T) = mpmath.nzeros(T) and
    theta = mpmath.siegeltheta.  Comparison side only.
    """
    if not _HAVE_MPMATH:
        return _CLASSICAL_S_FALLBACK, "fallback(documented)"
    mp.mp.dps = 25
    heights = [85.0, 95.0, 105.0, 115.0, 125.0, 135.0]
    vals = []
    for t in heights:
        n = mp.nzeros(t)
        theta = mp.siegeltheta(t)
        vals.append(abs(float(n) - 1.0 - float(theta) / math.pi))
    return max(vals), "mpmath"


# --- tests ---------------------------------------------------------------
def test_t1_propagator_reproduces_von_mangoldt() -> bool:
    """T1: the canonical propagator reproduces von Mangoldt at s=2."""
    bundle = build_prime_ladder_hamiltonian(n_primes=40, max_power=8)
    z_fn = weighted_spectral_trace(
        bundle.hamiltonian.H_freq, bundle.weight_operator, 2.0
    )
    spec = bundle.spectrum
    z_vec = float(np.sum(spec.weights * np.exp(-2.0 * spec.eigenvalues)))
    ref = _ref_log_zeta_deriv_2()
    same = abs(float(np.real(z_fn)) - z_vec) < 1e-9
    close = abs(z_vec - ref) < 1e-2
    ok = same and close
    print("  [T1] canonical propagator vs von Mangoldt")
    print(f"       weighted_spectral_trace(s=2) = {float(z_fn.real):.10f}")
    print(f"       vectorised propagator sum    = {z_vec:.10f}")
    print(f"       reference -zeta'/zeta(2)     = {ref:.10f}")
    print(f"       fn==vec (machine prec): {same} | ==ref: {close}")
    print(f"       -> {'PASS' if ok else 'FAIL'}")
    return ok


def test_t2_oscillation_in_kernel() -> bool:
    """T2: Im Z(1/2+iT) lands in ker(R-infinity)."""
    spec = build_prime_ladder_spectrum(N_PRIMES, max_power=MAX_POWER)
    t_grid = np.linspace(T_MIN, T_MAX, N_SAMPLES, endpoint=False)
    sig = _propagator_oscillation(spec.eigenvalues, spec.weights, t_grid)
    sig = sig - sig.mean()  # DC bin is trivially resonant
    rng, ker = split_residue_by_remesh_infinity(sig, tau_l=TAU_L, tau_g=TAU_G)
    e_total = float(np.dot(sig, sig))
    frac_range = float(np.dot(rng, rng)) / e_total
    frac_ker = float(np.dot(ker, ker)) / e_total
    recon = float(np.max(np.abs((rng + ker) - sig)))
    ok = (frac_range < THRESH) and (recon < 1e-9)
    print("  [T2] propagator oscillation -> R-infinity split")
    print(f"       energy in range (smooth/resonant): {frac_range:.6f}")
    print(f"       energy in kernel (oscillatory)   : {frac_ker:.6f}")
    print(f"       reconstruction error (range+ker) : {recon:.2e}")
    print(f"       range < {THRESH}: {frac_range < THRESH}")
    print(f"       -> {'PASS' if ok else 'FAIL'}")
    return ok


def test_t3_residue_is_sn_degenerate() -> bool:
    """T3: the kernel residue is S_n-degenerate (CCET, dynamical)."""
    t_grid = np.linspace(T_MIN, T_MAX, N_SAMPLES, endpoint=False)
    base = build_prime_ladder_spectrum(N_PRIMES, max_power=MAX_POWER)
    primes = np.asarray(base.primes, dtype=int)
    shuffled = primes.copy()
    np.random.default_rng(SEED).shuffle(shuffled)
    spec_b = build_prime_ladder_spectrum(
        N_PRIMES, max_power=MAX_POWER, primes=shuffled.tolist()
    )
    sa = _propagator_oscillation(base.eigenvalues, base.weights, t_grid)
    sb = _propagator_oscillation(spec_b.eigenvalues, spec_b.weights, t_grid)
    sa = sa - sa.mean()
    sb = sb - sb.mean()
    _, ker_a = split_residue_by_remesh_infinity(sa, tau_l=TAU_L, tau_g=TAU_G)
    _, ker_b = split_residue_by_remesh_infinity(sb, tau_l=TAU_L, tau_g=TAU_G)
    max_sig_diff = float(np.max(np.abs(sa - sb)))
    max_ker_diff = float(np.max(np.abs(ker_a - ker_b)))
    ok = max_ker_diff < 1e-9
    print("  [T3] S_n-degeneracy of the kernel residue")
    print(f"       max|sig_canonical - sig_shuffled| : {max_sig_diff:.2e}")
    print(f"       max|ker_canonical - ker_shuffled| : {max_ker_diff:.2e}")
    print("       (machine-precision zero == CCET: observable depends")
    print("        only on the spectrum SET, never on prime labels)")
    print(f"       -> {'PASS' if ok else 'FAIL'}")
    return ok


def test_t4_amplitude_scale_and_certificate() -> bool:
    """T4: S_TNFR matches S(T)'s amplitude SCALE; the wall is in the
    correlation structure (T3), not the magnitude. Plus R-infinity
    certificate cross-check.

    Note on P31: P31 reported that the *local* correction at zero
    ordinates needed a damping factor d ~ 3-5 to best-fit.  That is a
    PHASE-decoherence effect (the Fix(S_n)-perpendicular content T3
    shows is missing), NOT a global-amplitude deficit: the global
    amplitude scale matches, as confirmed here.
    """
    spec = build_prime_ladder_spectrum(N_PRIMES, max_power=MAX_POWER)
    t_grid = np.linspace(80.0, 140.0, 4096, endpoint=False)
    s_tnfr = _s_tnfr_oscillation(spec.eigenvalues, spec.weights, t_grid)
    max_tnfr = float(np.max(np.abs(s_tnfr)))
    classical, src = _classical_max_abs_s()
    ratio = classical / max_tnfr if max_tnfr > 0 else float("inf")
    same_scale = 0.25 < ratio < 4.0  # same order of magnitude
    cert = compute_residue_split_certificate(
        n_primes=N_PRIMES, max_power=MAX_POWER, n_periods=64
    )
    cert_ok = cert.verdict == "RESIDUE_IN_KER_ONLY"
    ok = cert_ok and same_scale
    print("  [T4] amplitude scale + R-infinity certificate")
    print(f"       max|S_TNFR| (canonical spectrum) : {max_tnfr:.6f}")
    print(f"       max|S(T)|   ({src:>16}) : {classical:.6f}")
    print(f"       scale ratio  S(T)/S_TNFR         : {ratio:.2f}x")
    print(f"       same order of magnitude          : {same_scale}")
    print(f"       R-infinity certificate verdict   : {cert.verdict}")
    print(
        f"       ratio_in_kernel={cert.ratio_in_kernel:.6f} "
        f"ratio_in_range={cert.ratio_in_range:.6f}"
    )
    print("       (amplitude scale matches; the obstruction is the")
    print("        S_n-degenerate phase structure of T3, which is why")
    print("        P31's local correction needs damping d ~ 3-5)")
    print(f"       -> {'PASS' if ok else 'FAIL'}")
    return ok


def main() -> int:
    print("=" * 70)
    print("Camino 17 - Nodal-propagator residue bridge")
    print("Canonical redo of P31 / §13nonies.4 from nodal dynamics")
    print("=" * 70)
    if not _HAVE_TNFR:
        print(f"SKIP: TNFR canonical imports unavailable: {_IMPORT_ERROR}")
        return 0
    print(f"mpmath available: {_HAVE_MPMATH}")
    print(
        f"L = lcm(tau_l={TAU_L}, tau_g={TAU_G}) = "
        f"{math.lcm(TAU_L, TAU_G)}; "
        f"N_PRIMES={N_PRIMES}, MAX_POWER={MAX_POWER}, "
        f"N_SAMPLES={N_SAMPLES}"
    )
    print("-" * 70)

    results = [
        test_t1_propagator_reproduces_von_mangoldt(),
        test_t2_oscillation_in_kernel(),
        test_t3_residue_is_sn_degenerate(),
        test_t4_amplitude_scale_and_certificate(),
    ]
    n_pass = sum(1 for r in results if r)
    n_total = len(results)

    print("-" * 70)
    print(f"SUMMARY: {n_pass}/{n_total} tests PASS")
    print("")
    print("Interpretation (honest scope):")
    print("  * The oscillatory observable is now GENERATED by the")
    print("    canonical structural propagator e^{-s H_P14}, not read")
    print("    off the Riemann-Siegel template (the P31 methodological")
    print("    gap).")
    print("  * It lands in ker(R-infinity) AND is S_n-degenerate: the")
    print("    canonical machinery sees only the spectrum SET, so it")
    print("    sits in ker(R-infinity) cap Fix(S_n).")
    print("  * The true S(T) requires the Fix(S_n)-perpendicular part")
    print("    of the kernel; its absence shows up as the phase")
    print("    decoherence (P31 damping d ~ 3-5), not as an amplitude")
    print("    deficit -- the global amplitude scale matches.")
    print("")
    print("THESIS VERDICT: OPEN, by design. This closes nothing.")
    print("  G4 = RH OPEN. R+ and pi remain assumed.")
    print("  Result strengthens branch-B2 evidence (a new canonical")
    print("  operator, not derivable from the present catalog, would")
    print("  be required to reach S(T)). It does NOT prove B2/B3/RH.")
    print("=" * 70)
    return 0 if n_pass == n_total else 1


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