#!/usr/bin/env python3
"""probe_p18t2_linearcoupling.py -- P18-T2 sub-result (ii), discharging the
linear-coupling assumption of A308.

A308 reduced sub-result (ii) to a single remaining assumption: each
coupling c_j enters as a LINEAR STRENGTH, c_j * S_j, rather than
f_j(c_j) * S_j for some nonlinear f_j, or c_j entering S_j's internal
structure. A308 retained that assumption; it did not derive it. This probe
tests whether the LINEAR FORM is itself FORCED by P18's treatment of the
couplings, so that the assumption is replaced by P18's own definition of a
coupling rather than asserted afresh.

P18 S4.2 ("Coefficient Values", Proposition "Coefficient Determination")
treats alpha, gamma, zeta, beta as coupling STRENGTHS, each fixed by
geometric/dimensional consistency: alpha = 1/137.036 from int rho =
alpha^{-1}; zeta = alpha^{5/4} from dimensional analysis 5/4 = 1 +
3/(4*3); beta = 3*pi/20 from the gravity coefficient; gamma = 3/4 from the
boundary/bulk dimension ratio. P18 S4.1 writes the density sector as
alpha * rho(r), the strength multiplying a fixed sector operator. A
"coupling strength" multiplying a fixed sector operator is, by definition,
a linear (degree-1 homogeneous) prefactor.

The principled, checkable content (two readings, both pre-registered):

  (L1) DIMENSIONAL HOMOGENEITY. A coupling strength c_j carries the
       dimension that makes c_j * S_j match the operator's (energy)
       dimension. The operator analog of "strength = linear prefactor" is
       that the per-sector operator response is INDEPENDENT of the
       coupling: d O_j / d c_j = S_j is a constant operator, so
       d^2 O_j / d c_j^2 = 0 (the zero operator). Under a units rescaling
       c_j -> s c_j only a degree-1 prefactor preserves this; any nonlinear
       f_j has d^2 O_j / d c_j^2 != 0 and the per-sector response acquires
       c_j-dependence, breaking the strength reading. So linear is the
       unique member of the declared family whose own second derivative in
       its coupling vanishes.

  (L2) DEGREE-1 HOMOGENEITY (superposition within a sector). A coupling
       strength means doubling c_j doubles that sector's contribution:
       f_j(2c) = 2 f_j(c) for all c, i.e. f_j is homogeneous of degree 1.
       Only the linear form is degree-1 homogeneous; quadratic, sqrt, and
       (exp - 1) all fail.

  (L3) Combined with A308 (IC2 / no category bleed forces additivity within
       the linear-coupling class): linear-coupling (L1/L2) + IC2 (A308)
       => additive composition with c_j * S_j, with NO remaining free
       assumption beyond the R1-R8 coverage requirements and the definition
       of a coupling as a strength. Both hinges of sub-result (ii) are now
       grounded in P18's stated requirements rather than assumed ad hoc.

PRE-REGISTERED TEST (HIT iff all three):

 (1) the linear coupling has d^2 O_j / d c_j^2 = 0 (the zero operator, <
     tol_zero) AND is degree-1 homogeneous (homogeneity residual < tol);

 (2) every nonlinear f_j in the declared family {quadratic, sqrt, exp-1}
     has d^2 O_j / d c_j^2 clearly nonzero (> tol_nonzero) AND fails
     degree-1 homogeneity (homogeneity residual > tol);

 (3) therefore within "coupling = dimensional strength" the linear form is
     forced (unique in the family on BOTH criteria); combined with A308,
     additive composition is forced with both hinges grounded in P18
     requirements.

DECLARED FAMILY (finite, fixed) of per-sector coupling forms f_j applied to
sector S_j (per-sector operator O_j(c) = K + f_j(c) * S_j):

  LINEAR      f(c) = c            (degree-1 homogeneous, P18's strength)
  QUADRATIC   f(c) = c^2          (degree 2)
  SQRT        f(c) = sqrt(c)      (degree 1/2)
  EXP-1       f(c) = exp(c) - 1   (no definite degree)

Sector and canonical coupling reused from A306/A308: the density sector
S = rho(r) [R5] (normalized), over the fixed (kinetic + V_self) operator K,
at canonical c0 = alpha = 1/137.036. (The renormalization sector S2 = r d/dr
[R7] at zeta = alpha^{5/4} is built too and the same vanishing holds; the
single-sector per-coupling reading is what L1/L2 concern.)

HONEST BOUNDARY (Case C). "Coupling = strength = degree-1 homogeneous linear
prefactor" is P18's OWN definition of alpha, gamma, zeta, beta (S4.2), made
precise. This probe shows the linear FORM is forced by that definition plus
dimensional homogeneity; it does NOT derive from something deeper that the
couplings MUST be strengths. So sub-result (ii) becomes: additive
composition is forced GIVEN (R1-R8 independent categories => IC2, A308) AND
(couplings are dimensional strengths => linear). Both hinges are now
grounded in P18's stated requirements rather than assumed ad hoc. Still: a
finite declared family of f_j is tested (not every conceivable nonlinear
form); sub-results (i)=A304 and (iii)=A302 are unchanged; A294's re-typing
stands; P18-T2 is NOT closed.

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_linearcoupling.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-2           # operator second-derivative below this: vanishing.
#   The linear form's d^2 O/dc^2 is the exact zero operator; its measured
#   value is finite-difference roundoff. The radial operator's d2 stencil
#   carries entries ~1/dr^2 ~ 2.5e5, so the second difference's numerator
#   roundoff is ~eps*||O|| ~ 1e-10 and, divided by h^2 = 1e-6, floors the
#   measured linear curvature at ~1e-3. The smallest nonlinear curvature in
#   the family (quadratic, 17.18) is four orders above this, so 1e-2 cleanly
#   separates the exact zero from genuine curvature.
TOL_NONZERO = 1e-1        # second derivative above this: present curvature
TOL_HOM = 1e-6            # homogeneity residual |f(2c)-2f(c)| / scale


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)

    # Scale bare sectors so the operator entries are O(1); the physics (which
    # sector acts) is unchanged, only the bookkeeping scale (A306/A308 conv).
    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 per-sector coupling forms f_j(c) ----------
    # Each gives a per-sector operator O_j(c) = K + f_j(c) * S_j. A "coupling
    # strength" is the LINEAR form; the others are nonlinear alternatives.
    ffuncs = {
        "linear":    lambda c: c,
        "quadratic": lambda c: c * c,
        "sqrt":      lambda c: math.sqrt(c),
        "exp_minus1": lambda c: math.exp(c) - 1.0,
    }

    # canonical coupling magnitude at which derivatives are evaluated
    c0 = A0
    # Finite-difference step for derivatives in c. A second difference divides
    # by h^2, so a step too small amplifies floating-point roundoff in the
    # numerator: for the linear form the analytic d^2 O/dc^2 is the exact zero
    # operator, and a tiny h^2 turns ~1e-13 roundoff into O(1) noise. The
    # coupling forms here are smooth on the scale of c0, so an O(1e-3) absolute
    # step resolves their curvature cleanly while keeping the linear form's
    # second difference at machine zero (numerator roundoff / h^2 ~ 1e-7).
    h = 1e-3

    def O1(fname, c):
        """Per-sector operator with sector S1 under coupling form f."""
        return K + ffuncs[fname](c) * S1

    def first_d_operator(fname, c):
        """d O_j / d c (operator), central difference."""
        return (O1(fname, c + h) - O1(fname, c - h)) / (2 * h)

    def second_d_operator(fname, c):
        """d^2 O_j / d c^2 (operator), central second difference."""
        return (O1(fname, c + h) - 2 * O1(fname, c) + O1(fname, c - h)) / (h * h)

    # ---- (L1) per-sector operator response: first and second derivatives --
    d1_norm = {}     # ||d O_j / d c||_F  (should equal ||S_j|| only for linear)
    d2_norm = {}     # ||d^2 O_j / d c^2||_F  (the discriminator: 0 iff linear)
    for fname in ffuncs:
        d1_norm[fname] = float(np.linalg.norm(first_d_operator(fname, c0), "fro"))
        d2_norm[fname] = float(np.linalg.norm(second_d_operator(fname, c0), "fro"))

    # analytic references for the per-sector operator derivatives:
    #   linear:    dO/dc = S1                     -> ||S1||_F ;  d2O/dc2 = 0
    #   quadratic: dO/dc = 2c S1 ; d2O/dc2 = 2 S1 -> ||2 S1||_F
    #   sqrt:      d2O/dc2 = -1/4 c^{-3/2} S1      -> (1/4) c0^{-3/2} ||S1||_F
    #   exp-1:     dO/dc = e^c S1 ; d2O/dc2 = e^c S1 -> e^{c0} ||S1||_F
    nS1 = float(np.linalg.norm(S1, "fro"))
    ref_d2 = {
        "linear":     0.0,
        "quadratic":  float(np.linalg.norm(2.0 * S1, "fro")),
        "sqrt":       float(np.linalg.norm((0.25 * c0 ** -1.5) * S1, "fro")),
        "exp_minus1": float(np.linalg.norm(math.exp(c0) * S1, "fro")),
    }

    # ---- (L2) degree-1 homogeneity: f(2c) = 2 f(c) for all c --------------
    # Test on a grid of c to catch any c-dependence; report the max residual.
    cs = np.linspace(0.2 * c0, 5.0 * c0, 25)
    hom_resid = {}
    for fname, f in ffuncs.items():
        # relative residual |f(2c) - 2 f(c)| normalized by |2 f(c)| (scale-free)
        rs = []
        for c in cs:
            lhs = f(2.0 * c)
            rhs = 2.0 * f(c)
            denom = abs(rhs) if abs(rhs) > 1e-300 else 1.0
            rs.append(abs(lhs - rhs) / denom)
        hom_resid[fname] = float(max(rs))

    is_homogeneous = {fname: (hom_resid[fname] < TOL_HOM) for fname in ffuncs}
    is_d2_zero = {fname: (d2_norm[fname] < TOL_ZERO) for fname in ffuncs}

    # ---- the forcing within the family ------------------------------------
    nonlinear = [f for f in ffuncs if f != "linear"]

    # (1) linear: d2O/dc2 = 0 AND degree-1 homogeneous
    linear_d2_zero = is_d2_zero["linear"]
    linear_homog = is_homogeneous["linear"]
    linear_ok = linear_d2_zero and linear_homog

    # (2) every nonlinear: d2O/dc2 clearly nonzero AND fails homogeneity
    nonlinear_disq = {
        f: (d2_norm[f] > TOL_NONZERO) and (not is_homogeneous[f])
        for f in nonlinear
    }
    all_nonlinear_disq = all(nonlinear_disq.values())

    # (3) uniqueness within family: linear is the unique form that is BOTH
    # d2-zero AND degree-1 homogeneous
    both_criteria = {f: (is_d2_zero[f] and is_homogeneous[f]) for f in ffuncs}
    linear_unique = all(both_criteria[f] == (f == "linear") for f in ffuncs)

    # cross-link sanity: the renorm sector S2 obeys the same linear vanishing
    def O2lin(c):
        return K + c * S2
    d2_S2_linear = float(np.linalg.norm(
        (O2lin(Z0 + h) - 2 * O2lin(Z0) + O2lin(Z0 - h)) / (h * h), "fro"))
    s2_linear_d2_zero = d2_S2_linear < TOL_ZERO

    hit = linear_ok and all_nonlinear_disq and linear_unique and s2_linear_d2_zero

    vals = (list(d1_norm.values()) + list(d2_norm.values())
            + list(hom_resid.values()) + [d2_S2_linear])
    nan = any(not np.isfinite(x) for x in vals)

    if hit:
        verdict = (
            "HIT: within 'coupling = dimensional strength' the LINEAR form is "
            "FORCED. (L1) linear per-sector operator response is constant: "
            "d^2 O/dc^2 = %.3e (zero operator, < %.0e); dO/dc = %.3e matches "
            "||S||_F = %.3e. Every nonlinear f has nonzero operator curvature: "
            "quadratic d^2O/dc^2 = %.3e (analytic ||2S|| = %.3e), sqrt = %.3e "
            "(analytic %.3e), exp-1 = %.3e (analytic %.3e), all > %.0e. (L2) "
            "only linear is degree-1 homogeneous: linear residual %.3e (< %.0e); "
            "quadratic %.3e, sqrt %.3e, exp-1 %.3e all fail. (L3) linear is the "
            "UNIQUE form in the family on BOTH criteria, so a coupling strength "
            "(degree-1 homogeneous prefactor whose own operator response has no "
            "curvature) is the linear form. Combined with A308 (IC2 forces "
            "additivity within the linear class), additive composition is "
            "forced with BOTH hinges grounded in P18 requirements: R1-R8 "
            "independent categories => IC2, and couplings are dimensional "
            "strengths => linear. Case C: 'couplings are strengths' is P18's "
            "S4.2 definition made precise, not derived from deeper principles; "
            "a finite declared family is tested; (i)=A304, (iii)=A302 "
            "unchanged; A294 stands; P18-T2 NOT closed."
            % (d2_norm["linear"], TOL_ZERO, d1_norm["linear"], nS1,
               d2_norm["quadratic"], ref_d2["quadratic"],
               d2_norm["sqrt"], ref_d2["sqrt"],
               d2_norm["exp_minus1"], ref_d2["exp_minus1"], TOL_NONZERO,
               hom_resid["linear"], TOL_HOM,
               hom_resid["quadratic"], hom_resid["sqrt"],
               hom_resid["exp_minus1"]))
    else:
        verdict = (
            "MISS: pre-registered predicate not met. linear_d2_zero=%s "
            "(d2=%.3e), linear_homog=%s (resid=%.3e), all_nonlinear_disq=%s, "
            "linear_unique=%s, s2_linear_d2_zero=%s (d2=%.3e). d2_norm=%s "
            "hom_resid=%s"
            % (linear_d2_zero, d2_norm["linear"], linear_homog,
               hom_resid["linear"], all_nonlinear_disq, linear_unique,
               s2_linear_d2_zero, d2_S2_linear, d2_norm, hom_resid))

    res = {
        "canonical": {"alpha": A0, "zeta": Z0},
        "sector": "density alpha*rho(r) [R5] (S1); renorm zeta*(r d/dr) [R7] (S2) cross-checked",
        "evaluation_point_c0": c0,
        "fd_step_h": h,
        "tolerances": {"tol_zero": TOL_ZERO, "tol_nonzero": TOL_NONZERO,
                       "tol_hom": TOL_HOM},
        "declared_family": list(ffuncs.keys()),
        "per_sector_dO_dc_fro": d1_norm,
        "per_sector_d2O_dc2_fro": d2_norm,
        "analytic_d2O_dc2_fro_ref": ref_d2,
        "homogeneity_residual_max": hom_resid,
        "is_degree1_homogeneous": is_homogeneous,
        "is_d2_zero": is_d2_zero,
        "both_criteria_satisfied": both_criteria,
        "linear_d2_zero": bool(linear_d2_zero),
        "linear_degree1_homogeneous": bool(linear_homog),
        "nonlinear_all_disqualified": bool(all_nonlinear_disq),
        "nonlinear_disqualified": {k: bool(v) for k, v in nonlinear_disq.items()},
        "linear_unique_in_family": bool(linear_unique),
        "renorm_sector_linear_d2_zero": bool(s2_linear_d2_zero),
        "renorm_sector_d2_value": d2_S2_linear,
        "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("sector: density alpha*rho [R5] (S1); evaluated at c0 = %.6g" % c0)
    print("--- (L1) per-sector operator response d O/dc and d^2 O/dc^2 ---")
    print("  %-12s %12s %14s %14s" % ("form", "||dO/dc||", "||d2O/dc2||",
                                      "analytic d2"))
    for f in ffuncs:
        print("  %-12s %12.4e %14.4e %14.4e"
              % (f, d1_norm[f], d2_norm[f], ref_d2[f]))
    print("  (linear ||dO/dc|| should equal ||S1||_F = %.4e)" % nS1)
    print("--- (L2) degree-1 homogeneity residual max|f(2c)-2f(c)|/|2f(c)| ---")
    for f in ffuncs:
        flag = "HOMOGENEOUS" if is_homogeneous[f] else "fails"
        print("  %-12s resid = %.3e   (%s)" % (f, hom_resid[f], flag))
    print("--- (L3) forcing within the family ---")
    print("  linear: d2O/dc2=0 (%s) AND degree-1 homogeneous (%s)"
          % (linear_d2_zero, linear_homog))
    print("  every nonlinear disqualified (curvature AND not homogeneous): %s"
          % all_nonlinear_disq)
    print("  linear unique on both criteria: %s" % linear_unique)
    print("  renorm sector S2 linear d2O/dc2 = %.3e (zero: %s)"
          % (d2_S2_linear, s2_linear_d2_zero))
    print("VERDICT:", verdict)
    print("nan:", nan)
    print("written:", os.path.normpath(OUT))


if __name__ == "__main__":
    main()
