#!/usr/bin/env python3
# Copyright Léon Fernando Vlegels -- CC BY 4.0
"""probe_selection_dirac.py -- Probe for OI-287-1 (A318).

FIT-VS-FORCED probe. The entire A296-A317 thread reduced the master operator to a
SCALAR radial operator. A317 broke the A316 family degeneracy by realizing the
families as the omega^k fiber-momentum sectors, but BOTH corpus-faithful fiber
forms entered the phase as an ADDITIVE m^2 energy and OVERSHOT: the canonical map
gave mu = 2.29-2.43, tau ~ 21, L2 ~ 13, no peak. A317 named the untested entry:
the phase entering through the FIRST-ORDER (Dirac / off-diagonal) structure rather
than as additive energy -- a different reduction of the master operator (the full
B^4 spinor / Dirac problem of P18 S2 / S3.1) rather than the scalar radial
reduction the whole thread used.

THE CONSTRUCTION THIS PROBE BUILDS (the point of A318).

P18 S2.2 / S3 puts the master operator on SECTIONS OF A SPINOR BUNDLE
Gamma(S (x) E) over B^4 with boundary S^3 (P18 S2.2). The scalar radial operator
of A311-A317 is the square D_{B^4}^2 reduced to a single component; the genuine
object is the FIRST-ORDER Dirac operator D_{B^4} acting on a multi-component
spinor. We build the B^4 RADIAL DIRAC SYSTEM, not the scalar radial operator.

EXACT RADIAL DIRAC MATRIX (declared here, grounded in P18 S2.2/S3.1 + A314).
On a 4D ball the Dirac operator in radial form couples an upper and a lower radial
function (f(r), g(r)) through a FIRST-ORDER system. With kappa the boundary-Dirac
eigenvalue on S^3 (A314: the genuine S^3 Dirac spectrum is +-(j+3/2), j=0,1,2,...)
and c = 3/2 the fixed 4D dimension constant (A312/A314: c = Delta_psi = (d-1)/2,
d=4, the same R/4 = 3/2 Lichnerowicz constant A314 machine-verified), the radial
Dirac operator on the two-component spinor Psi = (f, g)^T is

    D ~ [ [  m(r) ,  -d/dr - (kappa + c)/r ],
          [  d/dr - (kappa - c)/r ,  -m(r)  ] ]

m(r) is the mass / potential term built from the canonical density: we take
m(r) = sqrt(V_eff(r)), V_eff(r) = rho'(r)^2/(2 Omega^2) + E_self r^2 (1-r)^2, the
SAME V_eff the scalar thread used (so D^2 returns the scalar V_eff on the
diagonal, the corpus-faithful link to A311-A317). The KEY structural difference
from A317: the family / phase label enters through kappa -- the boundary-Dirac
eigenvalue, which carries the omega^k / Z_3 fiber structure (A314: the boundary
Dirac tower +-(j+3/2) is the spinor lift of the omega^k layer sectors) -- and it
couples to the density through the FIRST-ORDER OFF-DIAGONAL terms -(kappa+c)/r and
(kappa-c)/r and the spinor bilinear overlap, NOT as an additive m^2 energy. The
phase enters the matrix element GENTLY through the upper/lower coupling and the
cross-term f*g, not as the additive m^2 that overshot in A317.

THE THREE FAMILIES (canonical map, declared before computing).
The three families are the three lowest boundary-Dirac sectors carrying the
omega^k labels. On S^3 the Dirac spectrum is +-(j+3/2), j = 0,1,2,... (A314); the
positive tower is kappa in {3/2, 5/2, 7/2, ...} for j = 0,1,2. Family n
(n = 1,2,3; lepton generation) -> the n-th boundary-Dirac sector,

    kappa_n = (n-1) + 3/2  =  3/2, 5/2, 7/2   for n = 1, 2, 3.

This is the corpus-canonical omega^k / Z_3 mapping: the boundary-Dirac tower
+-(j+3/2) is the spinor realization of the layer-cycle sectors (A314 established
the +-(j+3/2) tower as the S^3 Dirac spectrum and tied R/4=3/2 to the Z_3 layer
constant), so the n-th sector j=n-1 carries omega^{n-1}. These give three
GENUINELY DISTINCT radial Dirac problems (distinct kappa -> distinct off-diagonal
coupling -> distinct bound states), breaking the A316 degeneracy by construction
through the FIRST-ORDER structure rather than an additive m^2.

DENSITY OVERLAP AS SPINOR MATRIX ELEMENT.
For family n solve the radial Dirac system for its lowest bound state (f_n, g_n)
(lowest |eigenvalue| of the first-order operator). Form the degree-n density
overlap as the SPINOR BILINEAR in the r^3 measure:

    R[n] = <Psi_n | C[n-1] pi^n r^n | Psi_n>_{r^3} / <Psi_n|Psi_n>_{r^3} / |E_n|

with the spinor density. DECLARED CHOICE: the density is the Dirac NORM density
f^2 + g^2 (the standard positive spinor probability density Psi^dagger Psi), in
the r^3 measure; <Psi|.|Psi>_{r^3} = integral (f^2 + g^2)(.) r^3 dr. (The
alternative f*g cross-term density is reported alongside as a declared
comparison, NOT the canonical one; the canonical overlap is f^2+g^2.) C = (2,3,16),
target normalized (1, 1.5, 8). The first-order Dirac structure means the
upper/lower coupling carries the phase (kappa) into the matrix element through the
RELATIVE WEIGHT of f vs g (set by the off-diagonal (kappa +- c)/r terms), gently,
unlike the m^2 additive energy of A317.

Compute normalized diagonal, mu = ratio[1], tau = ratio[2], L2 to (1,1.5,8),
valid_peak (column-standardized argmax of column n == n for all n AND a free-Dirac
/ sine-box control fails the same test in the same measure).

BASE CONSTANTS (verbatim from A311/A317): Omega = 4 pi^3 + pi^2 + pi, E_self =
13.177, rho'(r) = 48 pi^3 r^2 + 6 pi^2 r + 2 pi, rho(r) = 16 pi^3 r^3 + 3 pi^2 r^2
+ 2 pi r, C = (2,3,16), N ~ 800. r^3 measure throughout. scipy unavailable; numpy
only (np.linalg.eig / eigh). Staggered grid for the first derivative so the
discretized Dirac operator is stable and its free-field boundary eigenvalues
approximate the +-(j+3/2) tower (A314 sanity check).

PRE-REGISTERED VERDICT (declared before any number is computed):
  HIT / FORCED (CLOSES the magnitude piece of OI-287-1) iff the CANONICAL family
    map (family n -> n-th boundary-Dirac sector kappa_n = (n-1)+3/2, the omega^k
    eigenspace) reproduces (1, 1.5, 8) within 20% on all three entries AND
    restores a valid diagonal peak, with the Dirac structure corpus-faithful and
    NO fitting.
  STRENGTHENED-PARTIAL iff the Dirac / first-order coupling moves mu into
    [1.3, 1.7] (gentler than A317's overshoot, the predicted improvement) and/or
    gives three distinct entries closer to target than the scalar radial
    L2 = 0.399, but does not fully hit / peak.
  NULL iff the Dirac reduction does not improve on the scalar radial closest
    (L2 = 0.399, mu = 1.222) -- then the magnitude is not in this reduction either.

CRITICAL HONESTY: the canonical map is family n -> n-th boundary-Dirac sector
kappa_n (the omega^k / Z_3 eigenspace). A HIT requires THAT map. A tuned kappa
assignment or a tuned m(r) is a FIT and is flagged. tau and peak are reported, not
just mu. The central comparison is stated plainly: is the first-order Dirac
coupling GENTLER than A317's m^2 overshoot (did mu come down from 2.3 toward 1.5)?
The canonical density choice is f^2 + g^2 (declared up front; f*g reported as a
comparison). No NaN; the radial Dirac discretization is stable (staggered grid);
the free-Dirac boundary eigenvalues are checked against the +-(j+3/2) tower (A314).
"""
import os
import json
import math

import numpy as np


PI = math.pi
OMEGA = 4 * PI**3 + PI**2 + PI            # alpha^{-1}
E_SELF = 13.177
C = (2.0, 3.0, 16.0)                     # raw target coefficients
TARGET = [c / C[0] for c in C]           # normalized: (1, 1.5, 8)
N = 800                                   # grid size (kept in 600-1200)
C_DIM = 1.5                               # 4D dimension constant c = (d-1)/2 = 3/2
                                          # = Delta_psi = R/4 (A312/A314)

# CANONICAL family-to-sector map: family n (n=1,2,3) -> n-th boundary-Dirac
# sector kappa_n = (n-1) + 3/2 (A314 tower +-(j+3/2), j = n-1), the omega^k
# eigenspace / Z_3 spinor sector.
CANONICAL_KAPPA = [1.5, 2.5, 3.5]        # j = 0,1,2 -> n = 1,2,3

HERE = os.path.dirname(os.path.abspath(__file__))
OUT = os.path.join(HERE, "..", "..", "..", "..", "Outputs",
                   "selection_rule_dirac.json")


def rho_p(x):
    return 48 * PI**3 * x**2 + 6 * PI**2 * x + 2 * PI


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


def Vpot(x):
    return rho_p(x)**2 / (2 * OMEGA**2) + E_SELF * x**2 * (1 - x)**2


def mass(x):
    """m(r) = sqrt(V_eff): the mass/potential term so that D^2 returns V_eff on
    the diagonal (corpus-faithful link to the scalar A311-A317 V_eff)."""
    return np.sqrt(Vpot(x))


# --- staggered radial Dirac system -----------------------------------------
# Upper component f on the primary grid r_i (i=0..n-1), lower component g on the
# staggered grid r_{i+1/2}. The first derivative d/dr is a centered difference
# between staggered nodes, which makes the discrete first-order Dirac operator
# stable (no fermion-doubling spurious mode at this resolution) and self-adjoint
# in the staggered inner product. Dirichlet f(1)=0; g vanishes by the staggered
# placement at the wall.

def grids():
    """Primary grid for f, staggered grid for g, both interior of (0,1)."""
    dr = 1.0 / N
    rf = np.array([(i + 1) * dr for i in range(N - 1)])     # f-nodes, r in (0,1)
    rg = np.array([(i + 0.5) * dr for i in range(N - 1)])   # g-nodes, staggered
    return rf, rg, dr


def dirac_matrix(kappa, m_on_f, m_on_g):
    """Assemble the (2n-1) x (2n-1) Hermitian radial Dirac matrix.

    Block structure on (f_0..f_{n-1}, g_0..g_{n-1}):
      D ~ [ [  diag(m_f) ,  -B - diag((kappa+c)/r_f) P_fg ],
            [  B^T - diag((kappa-c)/r_g) P_gf ,  -diag(m_g) ] ]
    with B the centered first-derivative coupling f<->g on the staggered grid.
    We build the full matrix and symmetrize (the radial Dirac operator is
    self-adjoint in r^3 dr; the staggered first-derivative is anti-symmetric so
    the off-diagonal blocks are transposes up to the 1/r terms, which we place
    symmetrically as the average of the two placements to keep Hermiticity).
    """
    rf, rg, dr = grids()
    n = N - 1
    M = np.zeros((2 * n, 2 * n))
    # diagonal mass blocks
    for i in range(n):
        M[i, i] = m_on_f[i]
        M[n + i, n + i] = -m_on_g[i]
    # first-derivative coupling: d/dr maps g (staggered) -> f (primary).
    # (d/dr g)(r_f_i) ~ (g_i - g_{i-1}) / dr ; its adjoint maps f -> g.
    # off-diagonal entries connect f_i with g_i and g_{i-1}.
    for i in range(n):
        # -d/dr acting from lower block: coefficient on g_i and g_{i-1}
        # centered staggered derivative at f-node i uses g_i (right) and
        # g_{i-1} (left).
        c_kp = (kappa + C_DIM) / rf[i]          # (kappa+c)/r on f-node
        # d/dr term: + g_i/dr - g_{i-1}/dr  (so -d/dr is the negative)
        # 1/r term split symmetrically onto the two staggered neighbours
        # off-diagonal block top-right: row i (f), col n+j (g)
        # contribution from -d/dr - (kappa+c)/r
        coupl_right = -(1.0 / dr) - 0.5 * c_kp      # to g_i
        coupl_left = (1.0 / dr) - 0.5 * c_kp        # to g_{i-1}
        M[i, n + i] += coupl_right
        if i - 1 >= 0:
            M[i, n + (i - 1)] += coupl_left
    # symmetrize: the lower-left block is the transpose (self-adjointness),
    # which automatically encodes +d/dr - (kappa-c)/r up to the symmetric 1/r
    # placement. Adding the transpose then halving the diagonal-mass double-count.
    diag_keep = np.diag(M).copy()
    M = M + M.T
    for i in range(2 * n):
        M[i, i] = diag_keep[i]      # restore single diagonal (transpose doubled)
    return M, rf, rg, dr


def solve_dirac(kappa):
    """Solve the radial Dirac eigenproblem for sector kappa. Returns the lowest
    |E| bound state as (f, g) on the f-grid (g interpolated to f-grid for the
    overlap), the eigenvalue E, and the f/g grids."""
    rf, rg, dr = grids()
    m_on_f = mass(rf)
    m_on_g = mass(rg)
    M, rf, rg, dr = dirac_matrix(kappa, m_on_f, m_on_g)
    w, v = np.linalg.eigh(M)
    # lowest-|E| state: the Dirac spectrum is +-; the bound state of interest is
    # the smallest-magnitude positive eigenvalue (the ground charged-lepton-like
    # mode). Select smallest |w|.
    idx = int(np.argmin(np.abs(w)))
    E = float(w[idx])
    vec = v[:, idx]
    n = N - 1
    f = vec[:n]
    g = vec[n:]
    # interpolate g (staggered) onto the f-grid for the bilinear overlap
    g_on_f = np.empty(n)
    g_on_f[1:] = 0.5 * (g[1:] + g[:-1])
    g_on_f[0] = g[0]
    return f, g_on_f, E, rf, dr


def overlap_spinor(f, g, rf, dr, E, density):
    """Degree-k normalized spinor density overlap, divided by |E|, in r^3 dr.

    density = 'norm' -> f^2 + g^2 (CANONICAL, Psi^dagger Psi)
    density = 'cross' -> 2 f g    (declared comparison, the Dirac cross-term)
    """
    w3 = rf**3
    if density == "norm":
        d = f**2 + g**2
    elif density == "cross":
        d = 2.0 * f * g
    else:
        raise ValueError(density)
    norm = float(np.sum(d * w3) * dr)
    if abs(norm) < 1e-300:
        return [float("nan")] * 3
    out = []
    for k in (1, 2, 3):
        num = float(np.sum(d * C[k - 1] * PI**k * rf**k * w3) * dr)
        out.append(num / norm / abs(E))
    return out


def node_count(col):
    s = np.sign(col[np.abs(col) > 1e-9])
    return int(np.sum(s[1:] != s[:-1]))


def diag_peaks(M):
    """Column-standardized argmax of column k == k for all k."""
    Z = np.zeros_like(M)
    for k in range(3):
        c = M[:, k]
        Z[:, k] = (c - c.mean()) / c.std() if c.std() > 0 else 0.0
    return all(int(np.argmax(Z[:, k])) == k for k in range(3))


def norm_diag(M):
    d = [M[n, n] for n in range(3)]
    return [v / d[0] for v in d]


def l2(ratio):
    return math.sqrt(sum((ratio[i] - TARGET[i])**2 for i in range(3)))


def free_dirac_boundary_eigs(kappa, n_show=4):
    """Free radial Dirac (m=0) boundary eigenvalues for sector kappa, as the
    A314 sanity check: with no mass the lowest |E| modes should approximate the
    boundary-Dirac tower spacing near +-(kappa) = +-(j+3/2). We report the
    smallest-magnitude positive eigenvalues of the FREE operator (V=0)."""
    rf, rg, dr = grids()
    zero_f = np.zeros(N - 1)
    zero_g = np.zeros(N - 1)
    M, rf, rg, dr = dirac_matrix(kappa, zero_f, zero_g)
    w = np.linalg.eigvalsh(M)
    pos = np.sort(w[w > 1e-9])
    return [round(float(x), 4) for x in pos[:n_show]]


def build_family_matrix(density):
    """Build the 3x3 B-matrix for the canonical family map under one density
    choice. Family n -> kappa_n; lowest bound state; degree-k overlap row."""
    states_f = []
    states_g = []
    Es = []
    nodes = []
    for n in range(3):
        kap = CANONICAL_KAPPA[n]
        f, g, E, rf, dr = solve_dirac(kap)
        states_f.append(f)
        states_g.append(g)
        Es.append(E)
        nodes.append(node_count(f))
    M = np.zeros((3, 3))
    for n in range(3):
        row = overlap_spinor(states_f[n], states_g[n], rf, dr, Es[n], density)
        M[n, :] = row
    return M, Es, nodes, rf, dr


def box_states(rf):
    """Free-Dirac / sine-box control: sqrt(2) sin((n+1) pi r) as the upper
    component with a vanishing lower component (the trivial / non-Dirac basis),
    in the same r^3 measure. The control must FAIL the peak test."""
    psi = np.column_stack([math.sqrt(2.0) * np.sin((n + 1) * PI * rf)
                           for n in range(3)])
    lam = np.array([(n + 1)**2 * PI**2 for n in range(3)])
    return psi, lam


def box_matrix(psi, lam, rf, dr):
    M = np.zeros((3, 3))
    w3 = rf**3
    for n in range(3):
        d = psi[:, n]**2
        norm = float(np.sum(d * w3) * dr)
        for k in (1, 2, 3):
            num = float(np.sum(d * C[k - 1] * PI**k * rf**k * w3) * dr)
            M[n, k - 1] = num / norm / lam[n]
    return M


def evaluate(density, rf_box, psi_box, lam_box):
    M, Es, nodes, rf, dr = build_family_matrix(density)
    B_box = box_matrix(psi_box, lam_box, rf, dr)
    ratio = norm_diag(M)
    peak = bool(diag_peaks(M))
    box_peak = bool(diag_peaks(B_box))
    within20 = all(abs(ratio[i] - TARGET[i]) <= 0.20 * TARGET[i]
                   for i in range(3))
    diag_entries = [M[n, n] for n in range(3)]
    distinct = (abs(diag_entries[0] - diag_entries[1]) > 1e-6
                and abs(diag_entries[0] - diag_entries[2]) > 1e-6
                and abs(diag_entries[1] - diag_entries[2]) > 1e-6)
    # degeneracy guard (the A316 failure mode)
    degenerate = False
    for a in range(3):
        for b in range(a + 1, 3):
            if (abs(Es[a] - Es[b]) < 1e-6
                    and all(abs(M[a, q] - M[b, q]) < 1e-6 for q in range(3))):
                degenerate = True
    return {
        "density_choice": density,
        "canonical_kappa": CANONICAL_KAPPA,
        "diagonal_ratio": [round(float(v), 4) for v in ratio],
        "mu": round(float(ratio[1]), 4),
        "tau": round(float(ratio[2]), 4),
        "peaks_diagonal": peak,
        "box_control_peaks": box_peak,
        "valid_peak": bool(peak and not box_peak),
        "within_20pct": bool(within20),
        "hits_target": bool(within20 and peak and not box_peak),
        "L2_to_target": round(float(l2(ratio)), 4),
        "dirac_eigenvalues": [round(float(v), 4) for v in Es],
        "node_counts": nodes,
        "degenerate_selection": degenerate,
        "three_distinct_entries": bool(distinct),
        "mu_in_partial_band": bool(1.3 <= ratio[1] <= 1.7),
    }


def main():
    rf, rg, dr = grids()
    psi_box, lam_box = box_states(rf)

    # A314 sanity check: free-Dirac boundary eigenvalues per sector
    free_eigs = {f"kappa_{k}": free_dirac_boundary_eigs(k)
                 for k in CANONICAL_KAPPA}

    densities = {}
    for dens in ("norm", "cross"):
        densities[dens] = evaluate(dens, rf, psi_box, lam_box)

    # the A317 / scalar-radial bar to beat
    A317_overshoot_mu = 2.29           # F-Hopf mu (closest of A317's two)
    A317_overshoot_mu_max = 2.43       # F-lap mu
    SCALAR_closest_L2 = 0.399          # A316/A311 scalar-radial closest
    SCALAR_closest_mu = 1.222

    canonical = densities["norm"]      # canonical density is f^2 + g^2
    cross = densities["cross"]

    # central comparison: is the first-order coupling gentler than A317's m^2?
    gentler_than_A317 = bool(canonical["mu"] < A317_overshoot_mu)
    # did mu come down toward 1.5 from A317's overshoot?
    mu_came_down = bool(abs(canonical["mu"] - 1.5)
                        < abs(A317_overshoot_mu - 1.5))
    improves_scalar_L2 = bool(canonical["L2_to_target"]
                              < SCALAR_closest_L2 - 1e-9)
    moves_mu_vs_scalar = bool(abs(canonical["mu"] - 1.5)
                              < abs(SCALAR_closest_mu - 1.5) - 1e-9)

    degeneracy_broken = bool(canonical["three_distinct_entries"]
                             and not canonical["degenerate_selection"])

    # NaN guard
    nan_found = False
    for d in densities.values():
        for val in d["diagonal_ratio"] + d["dirac_eigenvalues"]:
            if val != val:
                nan_found = True

    # ---- VERDICT (canonical density f^2+g^2, canonical kappa map) -----------
    if canonical["hits_target"]:
        verdict = "FORCED"
    elif (degeneracy_broken
          and (canonical["mu_in_partial_band"]
               or improves_scalar_L2)):
        verdict = "STRENGTHENED_PARTIAL"
    elif not (gentler_than_A317 or improves_scalar_L2):
        verdict = "NULL"
    else:
        # gentler than A317 but not into band and not beating scalar L2: an
        # honest weak partial. Report as STRENGTHENED_PARTIAL only if it beats
        # the scalar closest OR lands in band; otherwise NULL on the magnitude.
        verdict = "NULL"

    summary = {
        "addendum": "A318",
        "addresses": "OI-287-1",
        "builds_on": ["A311", "A314", "A316", "A317"],
        "thesis": ("the phase enters through the FIRST-ORDER Dirac / off-diagonal "
                   "structure, not as additive m^2 energy. Build the B^4 radial "
                   "Dirac system (P18 S2.2/S3.1) on the two-component spinor "
                   "(f,g); family n -> n-th boundary-Dirac sector kappa_n = "
                   "(n-1)+3/2 (A314 tower +-(j+3/2), the omega^k / Z_3 spinor "
                   "sector); density overlap = spinor bilinear f^2+g^2 in r^3 dr "
                   "/ |E|. The off-diagonal (kappa+-c)/r coupling carries the "
                   "phase into the matrix element gently, unlike A317's m^2."),
        "target_normalized": TARGET,
        "constants": {"OMEGA_alpha_inv": OMEGA, "E_self": E_SELF,
                      "c_dim_constant": C_DIM},
        "canonical_kappa_map": CANONICAL_KAPPA,
        "canonical_density_choice": "norm (f^2 + g^2, Psi^dagger Psi)",
        "comparison_density_choice": "cross (2 f g, Dirac cross-term)",
        "free_dirac_boundary_eigenvalues": free_eigs,
        "A314_boundary_tower": "+-(j+3/2) = +-(3/2, 5/2, 7/2, ...)",
        "densities": densities,
        # canonical (f^2+g^2) reported entries
        "canonical_diagonal": canonical["diagonal_ratio"],
        "canonical_mu": canonical["mu"],
        "canonical_tau": canonical["tau"],
        "canonical_L2": canonical["L2_to_target"],
        "canonical_valid_peak": canonical["valid_peak"],
        "canonical_three_distinct": canonical["three_distinct_entries"],
        "canonical_eigenvalues": canonical["dirac_eigenvalues"],
        # cross-term comparison
        "cross_diagonal": cross["diagonal_ratio"],
        "cross_mu": cross["mu"],
        "cross_L2": cross["L2_to_target"],
        # the central comparisons
        "A317_overshoot_mu_F_Hopf": A317_overshoot_mu,
        "A317_overshoot_mu_F_lap": A317_overshoot_mu_max,
        "scalar_radial_closest_L2": SCALAR_closest_L2,
        "scalar_radial_closest_mu": SCALAR_closest_mu,
        "first_order_gentler_than_A317_m2": gentler_than_A317,
        "mu_came_down_toward_1.5_vs_A317": mu_came_down,
        "improves_on_scalar_radial_L2": improves_scalar_L2,
        "moves_mu_toward_1.5_vs_scalar": moves_mu_vs_scalar,
        "A316_degeneracy_broken": degeneracy_broken,
        "nan_found": nan_found,
        "closes_OI_287_1_magnitude": bool(verdict == "FORCED"),
        "verdict": verdict,
    }

    os.makedirs(os.path.dirname(OUT), exist_ok=True)
    with open(OUT, "w") as fp:
        json.dump(summary, fp, indent=2)

    print("=" * 74)
    print("A318 PROBE -- OI-287-1: the FIRST-ORDER B^4 radial DIRAC system. Does")
    print("   the phase entering through the off-diagonal (kappa+-c)/r coupling")
    print("   (rather than A317's additive m^2) supply (1,1.5,8) with a peak?")
    print("=" * 74)
    print("target (normalized):", TARGET)
    print("canonical map: family n -> n-th boundary-Dirac sector "
          "kappa_n = (n-1)+3/2 = [1.5, 2.5, 3.5] (A314 +-(j+3/2) tower)")
    print("canonical density: f^2 + g^2 (Psi^dagger Psi); comparison: 2 f g")
    print(f"A317 m^2 overshoot: mu = {A317_overshoot_mu} (F-Hopf) / "
          f"{A317_overshoot_mu_max} (F-lap)")
    print(f"scalar-radial closest: L2 = {SCALAR_closest_L2}, "
          f"mu = {SCALAR_closest_mu}")
    print("-" * 74)
    print("A314 sanity -- free radial Dirac boundary eigenvalues (lowest |E|+):")
    for k, ev in free_eigs.items():
        print(f"   {k}: {ev}   (tower target near +-(j+3/2))")
    print("-" * 74)
    for dens in ("norm", "cross"):
        v = densities[dens]
        tag = "CANON" if dens == "norm" else "compr"
        print(f"[{tag}] density={dens:6s} diag={v['diagonal_ratio']} "
              f"mu={v['mu']} tau={v['tau']} peak={v['valid_peak']} "
              f"<=20%={v['within_20pct']} HIT={v['hits_target']} "
              f"L2={v['L2_to_target']}")
        print(f"        E={v['dirac_eigenvalues']} nodes={v['node_counts']} "
              f"3distinct={v['three_distinct_entries']} "
              f"degen={v['degenerate_selection']}")
    print("-" * 74)
    print(f"A316 degeneracy broken (3 distinct entries, canonical): "
          f"{degeneracy_broken}")
    print(f"first-order coupling GENTLER than A317's m^2 (mu < 2.29): "
          f"{gentler_than_A317}")
    print(f"mu came down toward 1.5 vs A317 overshoot: {mu_came_down}")
    print(f"improves on scalar-radial L2 (<0.399): {improves_scalar_L2}")
    print(f"moves mu toward 1.5 vs scalar (1.222): {moves_mu_vs_scalar}")
    print(f"NaN found: {nan_found}")
    print(f"\nVERDICT: {verdict}")
    print("-" * 74)
    print("CENTRAL COMPARISON (first-order Dirac vs A317 m^2):")
    print(f"   A317 m^2 additive energy: mu = 2.29-2.43 (overshoot)")
    print(f"   A318 first-order coupling: mu = {canonical['mu']}  "
          f"({'GENTLER' if gentler_than_A317 else 'NOT gentler'})")
    if verdict == "FORCED":
        print("  -> the CANONICAL Dirac map (family n -> kappa_n, f^2+g^2) "
              "reproduces")
        print("     (1,1.5,8) within 20% with a valid peak, Dirac structure")
        print("     corpus-faithful, NO fitting. OI-287-1 magnitude CLOSES.")
    elif verdict == "STRENGTHENED_PARTIAL":
        print("  -> the first-order Dirac coupling moves mu gentler than A317's")
        print("     m^2 and/or beats the scalar-radial closest, but does NOT")
        print("     fully hit (1,1.5,8) with a peak under the canonical map.")
    else:
        print("  -> the Dirac first-order reduction does NOT improve on the")
        print("     scalar-radial closest (L2=0.399, mu=1.222) on the magnitude.")
        print("     The magnitude is not in this reduction either; redirect.")
    print(f"wrote {OUT}")


if __name__ == "__main__":
    main()
