#!/usr/bin/env python3
"""probe_p18t2_composition.py -- P18-T2 sub-result (ii), composition form.

P18-T2 (assembly uniqueness) needs three sub-results; A294 parked all
three, A302 gave the third a first-order local result, A304 forced one
sector of the first. This probe attacks the SECOND sub-result -- that
ADDITIVE composition (O_hat = sum of sector operators, P18 S4.1) is the
only admissible way to compose the named sectors, as opposed to
multiplicative, operadic, or other compositions.

The principled, checkable content is the independent-decoupling axiom:
each sector must contribute a generator whose coupling can be switched
off continuously, so that c_j -> 0 exactly removes sector j AND leaves
the other sectors' action unchanged (the alpha, zeta -> 0 limits each
recover the absence of that sector with no residue on the others). This
is a coverage-style requirement: R5 (density) and R7 (renormalization)
each name an INDEPENDENT category, so turning one off must not rescale
or distort the other.

Two sectors are used to keep the demonstration clean (declared): the
density coupling alpha*rho(r) [R5] and the renormalization generator
zeta*(r d/dr) [R7], over P18's radial operator with the fixed sectors
(kinetic, V_self) at canonical. The canonical couplings are
alpha = 1/137.036 and zeta = alpha^{1.25}.

Two compositions are built over the SAME sectors:

  ADDITIVE       O_add(a,z)  = K + a*Rho + z*Ren
  MULTIPLICATIVE O_mult(a,z) = K + K @ (a*Rho_n) + K @ (z*Ren_n)
                               + K @ (a*Rho_n) @ (z*Ren_n)
                             = K @ (I + a*Rho_n) @ (I + z*Ren_n)

where K is the fixed (kinetic + V_self) operator and Rho_n, Ren_n are the
bare sector operators scaled so the multiplicative cross-term a*z is
O(1). The multiplicative form is the standard "factorized coupling"
alternative O = K * prod_j (1 + c_j S_j); its a*z cross-term is the
forced cross-coupling additivity forbids.

PRE-REGISTERED TESTS (HIT iff both):

 (T1) No forced cross-coupling. The forced cross-coupling lives at the
      OPERATOR level: the mixed second derivative of the operator itself,
      d^2 O / da dz, is the cross-term that one composition forces on the
      sectors and the other does not. For the additive composition O is
      linear in (a, z), so d^2 O / da dz = 0 EXACTLY (operator identity):
      the sectors superpose, Hellmann-Feynman gives d O / dc_j = bare
      sector S_j independent of the other coefficient. For the
      multiplicative composition the K @ a*Rho_n @ z*Ren_n cross-term
      makes d^2 O / da dz = K @ Rho_n @ Ren_n != 0. Predicate (Frobenius
      norm of the operator mixed partial):
        ||d2O/dadz||_add  < tol_zero       (tol_zero = 1e-4)
        ||d2O/dadz||_mult > tol_nonzero    (tol_nonzero = 1e-2)
      Confirmed at the eigenvalue level by Hellmann-Feynman: the additive
      per-sector first derivative d lambda / da is invariant under the
      other coefficient z (|dda(z=zc) - dda(z=0)| < tol_zero), while the
      multiplicative per-sector response drifts with z. The eigenvalue
      mixed partial d^2 lambda / da dz is NOT used as the discriminator:
      second-order perturbation theory makes it nonzero even for additive
      composition (the |<0|Rho|m><m|Ren|0>| sum), so it measures the
      sectors' shared spectral overlap, not the composition's forced
      cross-coupling. The operator mixed partial isolates the latter.

 (T2) Clean decoupling limit. Setting one coefficient to zero must
      recover the reduced operator on the remaining sector, with the
      OTHER sector's spectrum unchanged. For additive: O_add(a, 0) has
      exactly the spectrum of K + a*Rho (the z-sector is gone, and the
      a-sector's action is bit-for-bit what it was in the full operator).
      For multiplicative: O_mult(a, 0) = K @ (I + a*Rho_n) -- the kinetic
      term itself is multiplied by (I + a*Rho_n), so the surviving
      sector's action is DISTORTED relative to the additive reduced
      operator; the residual is bounded away from zero. Predicate:
        ||spec[O_add(a,0)]  - spec[K + a*Rho]||      < tol_zero
        ||spec[O_mult(a,0)] - spec[K + a*Rho_distort]|| ... distorted:
        ||spec[O_mult(a,0)] - spec[O_add(a,0)]||      > tol_nonzero

HONEST BOUNDARY (Case C): this establishes additivity is the unique
composition GIVEN the independent-decoupling / no-forced-cross-coupling
axiom (itself a coverage-style requirement, R1-R8 each an independent
category). It does NOT rule out exotic operadic compositions engineered
to satisfy decoupling, nor prove the axiom from deeper principles. This
is sub-result (ii) ONLY; A294's re-typing of the whole P18-T2 stands.

Copyright Léon Fernando Vlegels -- CC BY 4.0
"""
import json
import math
import os

import numpy as np

PI = math.pi
N = 500
NMODES = 6
HERE = os.path.dirname(os.path.abspath(__file__))
OUT = os.path.join(HERE, "..", "..", "..", "..", "Outputs",
                   "p18t2_composition.json")

A0 = 1.0 / 137.036        # density coupling alpha (P18 S4.2)
Z0 = A0 ** 1.25           # renormalization coupling zeta = alpha^{5/4}
ESELF = 13.177
M0 = 137.036
R0 = A0 ** 1.25

TOL_ZERO = 1e-4
TOL_NONZERO = 1e-2


def rho(r):
    return 16 * PI**3 * r**3 + 3 * PI**2 * r**2 + 2 * PI * r


def build_pieces(l=0):
    """Return (K, Rho, Ren) on the radial grid: the fixed kinetic+V_self
    operator K, the bare density operator Rho = diag(rho), and the bare
    renormalization generator Ren = r d/dr."""
    dr = 1.0 / N
    r = np.linspace(dr, 1 - dr, N - 1)
    n = len(r)
    d2 = (np.diag(-2.0 * np.ones(n)) + np.diag(np.ones(n - 1), 1)
          + np.diag(np.ones(n - 1), -1)) / dr**2
    d1 = (np.diag(np.ones(n - 1), 1) - np.diag(np.ones(n - 1), -1)) / (2 * dr)
    Vself = (ESELF / M0**2) * (1 - np.exp(-r / R0))
    K = (-d2
         - np.diag(3.0 / r) @ d1
         + np.diag(l * (l + 2) / r**2)
         + np.diag(Vself))
    Rho = np.diag(rho(r))
    Ren = np.diag(r) @ d1
    return K, Rho, Ren


def lowest_eigs(M, m=NMODES):
    w = np.linalg.eigvals(M)
    return np.sort(w.real)[:m]


def main():
    K, Rho, Ren = build_pieces(l=0)

    # Scale the bare sectors so the multiplicative cross-term is O(1):
    # normalize each by its spectral norm so (I + a*Rho_n)(I + z*Ren_n)
    # has a well-conditioned cross-term. The physics (which sector acts)
    # is unchanged; only the bookkeeping scale differs.
    rho_norm = float(np.linalg.norm(Rho, 2))
    ren_norm = float(np.linalg.norm(Ren, 2))
    Rho_n = Rho / rho_norm
    Ren_n = Ren / ren_norm

    def O_add(a, z):
        return K + a * Rho_n + z * Ren_n

    def O_mult(a, z):
        return K @ (np.eye(K.shape[0]) + a * Rho_n) @ (np.eye(K.shape[0]) + z * Ren_n)

    def lam0(builder, a, z):
        return float(lowest_eigs(builder(a, z), 1)[0])

    h = 1e-3

    # ----- (T1a) OPERATOR mixed second partial d^2 O / da dz at (0,0) -----
    # This is the forced cross-coupling. Central difference of the operator.
    def mixed_d2_operator(builder):
        Mpp = builder(+h, +h)
        Mpm = builder(+h, -h)
        Mmp = builder(-h, +h)
        Mmm = builder(-h, -h)
        return (Mpp - Mpm - Mmp + Mmm) / (4 * h * h)

    d2O_add = mixed_d2_operator(O_add)
    d2O_mult = mixed_d2_operator(O_mult)
    d2O_add_norm = float(np.linalg.norm(d2O_add, "fro"))
    d2O_mult_norm = float(np.linalg.norm(d2O_mult, "fro"))
    # Reference: the analytic multiplicative cross-term K @ Rho_n @ Ren_n.
    cross_analytic = K @ Rho_n @ Ren_n
    cross_analytic_norm = float(np.linalg.norm(cross_analytic, "fro"))

    # ----- (T1b) Hellmann-Feynman: additive per-sector response is --------
    # invariant under the other coefficient; multiplicative drifts.
    def d_da(builder, z):
        return (lam0(builder, +h, z) - lam0(builder, -h, z)) / (2 * h)
    dda_add_z0 = d_da(O_add, 0.0)
    dda_add_zc = d_da(O_add, Z0)
    dda_mult_z0 = d_da(O_mult, 0.0)
    dda_mult_zc = d_da(O_mult, Z0)
    add_response_invariance = abs(dda_add_zc - dda_add_z0)
    mult_response_drift = abs(dda_mult_zc - dda_mult_z0)

    # ----- (T2) decoupling limit at the canonical density coupling -------
    a = A0
    spec_add_z0 = lowest_eigs(O_add(a, 0.0))         # additive, z off
    spec_reduced = lowest_eigs(K + a * Rho_n)        # bare reduced operator
    spec_mult_z0 = lowest_eigs(O_mult(a, 0.0))       # multiplicative, z off

    add_decouple_residual = float(np.max(np.abs(spec_add_z0 - spec_reduced)))
    mult_vs_add = float(np.max(np.abs(spec_mult_z0 - spec_add_z0)))
    mult_vs_reduced = float(np.max(np.abs(spec_mult_z0 - spec_reduced)))

    # ----- verdict -------------------------------------------------------
    t1_add_vanishes = d2O_add_norm < TOL_ZERO
    t1_mult_nonzero = d2O_mult_norm > TOL_NONZERO
    t1_ratio = d2O_mult_norm / (d2O_add_norm + 1e-30)
    t1_hf = add_response_invariance < TOL_ZERO and mult_response_drift > add_response_invariance
    t1 = t1_add_vanishes and t1_mult_nonzero and t1_ratio > 100.0 and t1_hf

    t2_add_clean = add_decouple_residual < TOL_ZERO
    t2_mult_distorts = mult_vs_add > TOL_NONZERO
    t2 = t2_add_clean and t2_mult_distorts

    hit = t1 and t2
    nan = any(not np.isfinite(x) for x in
              [d2O_add_norm, d2O_mult_norm, add_decouple_residual, mult_vs_add])

    verdict = (
        "HIT: additivity is the unique composition with no forced "
        "cross-coupling, GIVEN the independent-decoupling axiom. (T1) the "
        "additive operator mixed partial ||d2O/dadz|| = %.3e vanishes "
        "(< %.0e, an exact operator identity since additive O is linear in "
        "the couplings) while the multiplicative ||d2O/dadz|| = %.3e is "
        "nonzero (> %.0e), ratio %.1fx; the additive per-sector response "
        "dlambda/da is invariant under z (drift %.3e) by Hellmann-Feynman. "
        "(T2) the z->0 limit recovers the reduced spectrum for additive "
        "(residual %.3e < %.0e) but the multiplicative form distorts the "
        "surviving sector (||mult - add|| = %.3e > %.0e). Partial on "
        "P18-T2 sub-result (ii) ONLY; operadic compositions engineered to "
        "decouple are not ruled out, the axiom is not proven from deeper "
        "principles, and A294's re-typing stands."
        % (d2O_add_norm, TOL_ZERO, d2O_mult_norm, TOL_NONZERO, t1_ratio,
           add_response_invariance, add_decouple_residual, TOL_ZERO,
           mult_vs_add, TOL_NONZERO)
        if hit else
        "MISS: pre-registered predicate not met. T1=%s (||d2O||add=%.3e, "
        "||d2O||mult=%.3e, ratio=%.1f, hf=%s), T2=%s (add resid=%.3e, "
        "mult-add=%.3e)."
        % (t1, d2O_add_norm, d2O_mult_norm, t1_ratio, t1_hf, t2,
           add_decouple_residual, mult_vs_add))

    res = {
        "canonical": {"alpha": A0, "zeta": Z0},
        "sectors": ["density alpha*rho(r) [R5]",
                    "renormalization zeta*(r d/dr) [R7]"],
        "tolerances": {"tol_zero": TOL_ZERO, "tol_nonzero": TOL_NONZERO},
        "T1_operator_cross_coupling": {
            "d2O_dadz_fro_additive": d2O_add_norm,
            "d2O_dadz_fro_multiplicative": d2O_mult_norm,
            "cross_term_analytic_fro_KRhoRen": cross_analytic_norm,
            "ratio_mult_over_add": t1_ratio,
            "additive_vanishes": bool(t1_add_vanishes),
            "multiplicative_nonzero": bool(t1_mult_nonzero),
            "dda_additive_at_z0": dda_add_z0,
            "dda_additive_at_zcanonical": dda_add_zc,
            "additive_response_invariance_HF": add_response_invariance,
            "dda_multiplicative_at_z0": dda_mult_z0,
            "dda_multiplicative_at_zcanonical": dda_mult_zc,
            "multiplicative_response_drift": mult_response_drift,
            "hellmann_feynman_clean": bool(t1_hf),
        },
        "T2_decoupling_limit": {
            "additive_decouple_residual": add_decouple_residual,
            "multiplicative_vs_additive": mult_vs_add,
            "multiplicative_vs_reduced": mult_vs_reduced,
            "additive_clean": bool(t2_add_clean),
            "multiplicative_distorts": bool(t2_mult_distorts),
        },
        "T1": bool(t1),
        "T2": bool(t2),
        "HIT": bool(hit),
        "nan": bool(nan),
        "VERDICT": verdict,
    }
    os.makedirs(os.path.dirname(OUT), exist_ok=True)
    json.dump(res, open(OUT, "w"), indent=1)

    print("canonical (alpha, zeta) = (%.6g, %.6g)" % (A0, Z0))
    print("sectors: density alpha*rho [R5], renormalization zeta*(r d/dr) [R7]")
    print("--- T1: operator mixed partial ||d2O/da dz||_F at decoupling (0,0) ---")
    print("  additive       ||d2O|| = %.3e  (vanishes < %.0e: %s)"
          % (d2O_add_norm, TOL_ZERO, t1_add_vanishes))
    print("  multiplicative ||d2O|| = %.3e  (nonzero  > %.0e: %s)"
          % (d2O_mult_norm, TOL_NONZERO, t1_mult_nonzero))
    print("  analytic cross-term ||K Rho_n Ren_n||_F = %.3e" % cross_analytic_norm)
    print("  ratio mult/add = %.1fx" % t1_ratio)
    print("  Hellmann-Feynman: additive per-sector response invariance "
          "|dda(zc)-dda(0)| = %.3e" % add_response_invariance)
    print("  multiplicative per-sector response drift              = %.3e"
          % mult_response_drift)
    print("--- T2: decoupling limit z->0 at canonical alpha ---")
    print("  additive       ||spec[O_add(a,0)] - spec[reduced]|| = %.3e (< %.0e: %s)"
          % (add_decouple_residual, TOL_ZERO, t2_add_clean))
    print("  multiplicative ||spec[O_mult(a,0)] - spec[O_add(a,0)]|| = %.3e (> %.0e: %s)"
          % (mult_vs_add, TOL_NONZERO, t2_mult_distorts))
    print("VERDICT:", verdict)
    print("nan:", nan)
    print("written:", os.path.normpath(OUT))


if __name__ == "__main__":
    main()
