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

A306 forced additive composition GIVEN the independent-decoupling axiom.
A308 tests whether that axiom follows from P18's R1-R8 coverage
requirements, so it is REPLACED by a coverage requirement rather than
assumed. P18 S1.2: (R1)-(R8) are categorical coverage requirements, each
an INDEPENDENT category. The operator reading is "no category bleed"
(IC2): sector i's coupling does not enter sector j's operator
contribution, measured by the operator mixed partial

  B_ij = || d^2 O / dc_i dc_j ||_F .

Recomputed from scratch on two declared sectors of P18's radial operator:
density coupling alpha*rho(r) [R5] (S1) and renormalization generator
zeta*(r d/dr) [R7] (S2), over the fixed (kinetic + V_self) operator K, at
canonical alpha = 1/137.036, zeta = alpha^{1.25}.

A declared, finite, linear-coupling family is tested:
  ADD   K + c1 S1 + c2 S2
  MULT  K (I + c1 S1)(I + c2 S2)
  AFFP  K + c1 S1 + c2 S2 + c1 c2 S1 S2
  OPER  K + c1 S1 + c2 S2 + c1 c2 (S1 S2 + S2 S1)/2

Asserts: additive B_ij = 0 (no bleed) and IC1-clean; every non-additive
composition has B_ij clearly nonzero (category bleed); within the family
no-bleed <=> additive; the c2 -> 0 coverage check. Hence IC2 (a coverage
requirement) forces additivity with NO separate decoupling axiom.

Tolerances: tol_zero = 1e-4 (vanishing), tol_nonzero = 1e-2 (forced
bleed). Honest boundary: still assumes linear-strength couplings (P18's
own treatment), tests a finite declared family not every operad;
sub-results (i),(iii) unchanged; A294's re-typing stands. Strictly
stronger than A306 -- the assumed axiom is now grounded in R1-R8.

  S1  Sector pieces and the declared family        - checks 1-2
  S2  Category bleed B_ij at the operator level     - checks 3-5
  S3  Coverage on removal and the equivalence       - checks 6-8
"""
import sys
import math

import numpy as np

PASS = FAIL = 0


def check(n, desc, cond):
    global PASS, FAIL
    ok = bool(cond)
    PASS += ok
    FAIL += (not ok)
    print(f"  [{'PASS' if ok else 'FAIL'}] {n:>2}. {desc}")


PI = math.pi
N = 500
NMODES = 6
A0 = 1.0 / 137.036
Z0 = A0 ** 1.25
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):
    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))
    return K, np.diag(rho(r)), np.diag(r) @ d1


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


K, Rho, Ren = build_pieces(l=0)
EYE = np.eye(K.shape[0])
S1 = Rho / float(np.linalg.norm(Rho, 2))      # density [R5]
S2 = Ren / float(np.linalg.norm(Ren, 2))      # renorm [R7]


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):
    return K + c1 * S1 + c2 * S2 + c1 * c2 * (S1 @ S2)


def O_oper(c1, c2):
    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):
    return (builder(+h, +h) - builder(+h, -h)
            - builder(-h, +h) + builder(-h, -h)) / (4 * h * h)


bleed = {name: float(np.linalg.norm(mixed_d2_operator(b), "fro"))
         for name, b in family.items()}

ref_mult = float(np.linalg.norm(K @ S1 @ S2, "fro"))
ref_affp = float(np.linalg.norm(S1 @ S2, "fro"))

a = A0
spec_reduced = lowest_eigs(K + a * S1)
ic1 = {name: float(np.max(np.abs(lowest_eigs(b(a, 0.0)) - spec_reduced)))
       for name, b in family.items()}

nobleed = {name: (bleed[name] < TOL_ZERO) for name in family}
equivalence = all(nobleed[name] == (name == "additive") for name in family)
nonadd = [n for n in family if n != "additive"]

print("S1  Sector pieces and the declared family")
check(1, "two distinct sectors built: density S1 [R5], renorm S2 [R7] "
      "(norms > 0, S1 != S2)",
      np.linalg.norm(S1) > 0 and np.linalg.norm(S2) > 0
      and not np.allclose(S1, S2))
check(2, "declared linear-coupling family of 4 compositions built "
      "(additive, multiplicative, affine-product, operadic); each "
      "non-additive carries a nonzero c1*c2 cross-structure (differs from "
      "additive in its operator mixed partial)",
      len(family) == 4
      and not np.allclose(mixed_d2_operator(O_mult), 0.0)
      and not np.allclose(mixed_d2_operator(O_affp), 0.0)
      and not np.allclose(mixed_d2_operator(O_oper), 0.0)
      and np.allclose(mixed_d2_operator(O_add), 0.0))

print("S2  Category bleed B_ij = ||d2O/dc1dc2||_F at the operator level")
check(3, "additive B_ij = %.3e is the zero operator (no category bleed, "
      "< %.0e: additive O is affine in the couplings, exact identity)"
      % (bleed["additive"], TOL_ZERO),
      bleed["additive"] < TOL_ZERO)
check(4, "multiplicative B_ij = %.3e is forced bleed (> %.0e), matching "
      "analytic ||K S1 S2||_F = %.3e to 1e-6 rel"
      % (bleed["multiplicative"], TOL_NONZERO, ref_mult),
      bleed["multiplicative"] > TOL_NONZERO
      and abs(bleed["multiplicative"] - ref_mult) / ref_mult < 1e-6)
check(5, "affine-product B_ij = %.3e and operadic B_ij = %.3e are forced "
      "bleed (both > %.0e); affp matches analytic ||S1 S2||_F = %.3e to "
      "1e-5 rel (finite-difference O(h^2) residual)"
      % (bleed["affine_product"], bleed["operadic"], TOL_NONZERO, ref_affp),
      bleed["affine_product"] > TOL_NONZERO
      and bleed["operadic"] > TOL_NONZERO
      and abs(bleed["affine_product"] - ref_affp) / ref_affp < 1e-5)

print("S3  Coverage on removal (c2->0) and the equivalence")
check(6, "additive c2->0 leaves category 1 intact: spectrum equals bare "
      "reduced K + alpha S1 (IC1 residual %.3e < %.0e)"
      % (ic1["additive"], TOL_ZERO),
      ic1["additive"] < TOL_ZERO)
check(7, "every non-additive composition is disqualified -- bleed > %.0e "
      "OR fails coverage IC1 > %.0e (mult: bleed %.3e / IC1 %.3e; affp: "
      "bleed %.3e / IC1 %.3e; oper: bleed %.3e / IC1 %.3e)"
      % (TOL_NONZERO, TOL_NONZERO,
         bleed["multiplicative"], ic1["multiplicative"],
         bleed["affine_product"], ic1["affine_product"],
         bleed["operadic"], ic1["operadic"]),
      all((bleed[n] > TOL_NONZERO) or (ic1[n] > TOL_NONZERO)
          for n in nonadd))
check(8, "within the linear-coupling family no-bleed <=> additive, so IC2 "
      "(a P18 coverage requirement) forces additivity with NO separate "
      "decoupling axiom; honest boundary: linear-strength couplings "
      "assumed, finite declared family",
      equivalence and bleed["additive"] < TOL_ZERO
      and all(bleed[n] > TOL_NONZERO for n in nonadd))

print(f"\n{'='*60}\nRESULT: {PASS} PASS / {FAIL} FAIL")
sys.exit(0 if FAIL == 0 else 1)