#!/usr/bin/env python3
"""probe_p18t2_decoupling.py -- P18-T2 sub-result (ii), discharging the
decoupling axiom of A306.

A306 forced additive composition GIVEN one assumption: the
independent-decoupling axiom (each c_j -> 0 removes sector j and leaves
the other sectors unchanged). That axiom was assumed, not derived. This
probe tests whether it FOLLOWS from P18's R1-R8 coverage requirements
themselves, so that the axiom is replaced by a P18 coverage requirement
rather than assumed afresh.

P18 S1.2 states (R1)-(R8) are categorical coverage requirements, not
algebraic constraints: each names an INDEPENDENT category of structure
(bulk, boundary, fiber, self-lensing, density, Z_3 family,
renormalization, moment hierarchy). The operator reading of "independent
category" is "no category bleed":

  (IC1) coverage survives removal: C with c_j = 0 still covers the
        categories i != j (removing one sector does not remove another's
        contribution);
  (IC2) no category bleed: the presence/coupling of sector i does not
        alter WHICH category sector j covers -- the sector-j contribution
        to the operator is a function of c_j and S_j alone, not of c_i.

IC2 is the direct operator reading of "the categories are independent."
The operator-level measure of category bleed between two sectors is the
forced cross term

  B_ij = || d^2 O / dc_i dc_j ||_F        (Frobenius norm of the OPERATOR
                                           mixed partial, A306's correction:
                                           the operator, not eigenvalues).

B_ij = 0 means sector i's coupling does not enter sector j's operator
contribution: no bleed. B_ij != 0 means turning on sector i changes
sector j's contribution: bleed.

PRE-REGISTERED TEST (HIT iff all three):

 (1) additive composition has B_ij = 0 (no bleed) AND recovers each
     category in the c_j -> 0 limit (IC1 holds, IC2 holds);

 (2) every NON-additive composition in the declared linear-coupling
     family has B_ij clearly nonzero (category bleed) OR fails
     independent coverage IC1;

 (3) therefore within the linear-coupling class no-bleed <=> additive,
     so IC2 (a P18 coverage requirement) forces additivity WITHOUT a
     separate decoupling axiom.

DECLARED FAMILY (finite, fixed; each c_j a linear coupling STRENGTH
scaling its sector, per P18's treatment of alpha,gamma,zeta,beta as
multiplicative strengths):

  ADD   O = K + c1 S1 + c2 S2                  (additive, P18 S4.1)
  MULT  O = K (I + c1 S1)(I + c2 S2)           (factorized coupling)
  AFFP  O = K + c1 S1 + c2 S2 + c1 c2 S1 S2    (affine-with-product)
  OPER  O = K + c1 S1 + c2 S2 + c1 c2 (S1 S2 + S2 S1)/2  (operadic,
                                       symmetrized binary composition)

Two sectors (declared) of P18's radial operator: density coupling
c1 = alpha, S1 = rho(r) [R5]; renormalization c2 = zeta, S2 = r d/dr
[R7]; over the fixed (kinetic + V_self) operator K. Canonical couplings
alpha = 1/137.036, zeta = alpha^{1.25}.

HONEST BOUNDARY (Case C, strictly stronger than A306): this discharges
the decoupling axiom by REPLACING it with IC2 (no category bleed), the
direct operator reading of "R1-R8 are independent categories" -- a P18
coverage requirement, not a new assumption. It does NOT eliminate all
assumptions: (a) we still assume couplings act as linear strengths,
which is P18's own treatment; (b) we test a finite declared family of
compositions, not literally every conceivable operad; (c) sub-results
(i) and (iii) are unchanged and A294's re-typing stands. Strictly
stronger than A306 because the assumed decoupling axiom is now grounded
in R1-R8 rather than asserted.

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_decoupling.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           # operator/spectrum residual below this: vanishing
TOL_NONZERO = 1e-2        # cross term above this: forced bleed


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): fixed kinetic+V_self operator K, bare density
    operator Rho = diag(rho) [R5], bare renormalization generator
    Ren = r d/dr [R7]."""
    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)
    EYE = np.eye(K.shape[0])

    # Scale bare sectors so cross-terms are O(1); the physics (which sector
    # acts) is unchanged, only the bookkeeping scale differs (A306 convention).
    S1 = Rho / float(np.linalg.norm(Rho, 2))     # density sector [R5]
    S2 = Ren / float(np.linalg.norm(Ren, 2))     # renorm sector [R7]

    # ---- the declared family of compositions (linear-coupling class) ----
    def O_add(c1, c2):
        return K + c1 * S1 + c2 * S2

    def O_mult(c1, c2):
        return K @ (EYE + c1 * S1) @ (EYE + c2 * S2)

    def O_affp(c1, c2):                          # affine-with-product
        return K + c1 * S1 + c2 * S2 + c1 * c2 * (S1 @ S2)

    def O_oper(c1, c2):                           # operadic (symmetrized)
        return (K + c1 * S1 + c2 * S2
                + c1 * c2 * 0.5 * (S1 @ S2 + S2 @ S1))

    family = {"additive": O_add, "multiplicative": O_mult,
              "affine_product": O_affp, "operadic": O_oper}

    h = 1e-3

    def mixed_d2_operator(builder):
        """B_ij = d^2 O / dc1 dc2 via central differences of the OPERATOR."""
        return (builder(+h, +h) - builder(+h, -h)
                - builder(-h, +h) + builder(-h, -h)) / (4 * h * h)

    # ---- (1) operator-level category bleed B_ij for each composition ----
    bleed = {}
    for name, builder in family.items():
        B = mixed_d2_operator(builder)
        bleed[name] = float(np.linalg.norm(B, "fro"))

    # Reference cross-terms (analytic mixed partials):
    #   mult: d^2/dc1dc2 [K(I+c1 S1)(I+c2 S2)] = K S1 S2
    #   affp: d^2/dc1dc2 [... + c1 c2 S1 S2]   = S1 S2
    #   oper:                                  = (S1 S2 + S2 S1)/2
    ref_mult = float(np.linalg.norm(K @ S1 @ S2, "fro"))
    ref_affp = float(np.linalg.norm(S1 @ S2, "fro"))
    ref_oper = float(np.linalg.norm(0.5 * (S1 @ S2 + S2 @ S1), "fro"))

    # ---- IC1: coverage survives removal (c2 -> 0 still covers cat 1) ----
    # additive: O(c1, 0) = K + c1 S1 = bare reduced operator (sector 1 intact)
    a = A0
    spec_reduced = lowest_eigs(K + a * S1)
    ic1 = {}
    for name, builder in family.items():
        spec_c2off = lowest_eigs(builder(a, 0.0))
        # IC1 residual: does removing sector 2 leave sector 1's contribution
        # equal to the bare reduced operator (no distortion of category 1)?
        ic1[name] = float(np.max(np.abs(spec_c2off - spec_reduced)))

    # ---- the equivalence within the family ------------------------------
    # no-bleed (B_ij < tol_zero) holds for which compositions?
    nobleed = {name: (bleed[name] < TOL_ZERO) for name in family}
    is_additive = {name: (name == "additive") for name in family}
    # equivalence: no-bleed <=> additive across the declared family
    equivalence = all(nobleed[name] == is_additive[name] for name in family)

    # additive: no bleed AND clean coverage on removal (IC1, IC2)
    add_nobleed = bleed["additive"] < TOL_ZERO
    add_ic1_clean = ic1["additive"] < TOL_ZERO
    add_matches_ic2_and_ic1 = add_nobleed and add_ic1_clean

    # every non-additive: bleed clearly nonzero OR fails IC1
    nonadd = [n for n in family if n != "additive"]
    nonadd_disqualified = {
        n: (bleed[n] > TOL_NONZERO) or (ic1[n] > TOL_NONZERO)
        for n in nonadd
    }
    all_nonadd_disqualified = all(nonadd_disqualified.values())

    hit = (add_matches_ic2_and_ic1 and all_nonadd_disqualified
           and equivalence)

    vals = list(bleed.values()) + list(ic1.values())
    nan = any(not np.isfinite(x) for x in vals)

    if hit:
        verdict = (
            "HIT: within the linear-coupling family no-bleed <=> additive, "
            "so IC2 (no category bleed -- the operator reading of R1-R8 as "
            "independent categories) forces additivity WITHOUT a separate "
            "decoupling axiom. (1) additive B_ij = %.3e (no bleed, < %.0e) "
            "and IC1-clean (c2->0 residual %.3e < %.0e: category 1 intact). "
            "(2) every non-additive composition has category bleed: "
            "multiplicative B_ij = %.3e, affine-product B_ij = %.3e, "
            "operadic B_ij = %.3e (all > %.0e). (3) additivity is the unique "
            "no-bleed composition in the family, so the prior decoupling "
            "axiom is REPLACED by IC2, a P18 coverage requirement. Case C, "
            "strictly stronger than A306: still assumes linear-strength "
            "couplings (P18's own treatment), tests a finite declared "
            "family not every operad; sub-results (i),(iii) unchanged; "
            "A294 stands."
            % (bleed["additive"], TOL_ZERO, ic1["additive"], TOL_ZERO,
               bleed["multiplicative"], bleed["affine_product"],
               bleed["operadic"], TOL_NONZERO))
    else:
        verdict = (
            "MISS: pre-registered predicate not met. add_nobleed=%s "
            "(B=%.3e), add_ic1=%s (resid=%.3e), nonadd_disqualified=%s, "
            "equivalence=%s. bleed=%s ic1=%s"
            % (add_nobleed, bleed["additive"], add_ic1_clean,
               ic1["additive"], nonadd_disqualified, equivalence,
               bleed, ic1))

    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},
        "declared_family": list(family.keys()),
        "category_bleed_Bij_fro": bleed,
        "analytic_cross_term_fro": {
            "multiplicative_K_S1_S2": ref_mult,
            "affine_product_S1_S2": ref_affp,
            "operadic_sym_S1S2": ref_oper,
        },
        "IC1_coverage_residual_c2_to_0": ic1,
        "no_bleed_flags": nobleed,
        "no_bleed_iff_additive": bool(equivalence),
        "additive_IC2_and_IC1_clean": bool(add_matches_ic2_and_ic1),
        "nonadditive_all_disqualified": bool(all_nonadd_disqualified),
        "nonadditive_disqualified": {k: bool(v)
                                     for k, v in nonadd_disqualified.items()},
        "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("--- (1) operator-level category bleed B_ij = ||d2O/dc1dc2||_F ---")
    for name in family:
        flag = "no bleed" if nobleed[name] else "BLEED"
        print("  %-16s B_ij = %.3e   (%s)" % (name, bleed[name], flag))
    print("  analytic refs: mult ||K S1 S2|| = %.3e, affp ||S1 S2|| = %.3e, "
          "oper = %.3e" % (ref_mult, ref_affp, ref_oper))
    print("--- IC1 coverage on removal (c2->0, residual vs bare reduced) ---")
    for name in family:
        print("  %-16s IC1 residual = %.3e" % (name, ic1[name]))
    print("--- equivalence within the linear-coupling family ---")
    print("  no-bleed <=> additive : %s" % equivalence)
    print("VERDICT:", verdict)
    print("nan:", nan)
    print("written:", os.path.normpath(OUT))


if __name__ == "__main__":
    main()
