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

FIT-VS-FORCED probe. The full master operator (P18 S2, S4.3) acts on sections of
B^4 = (radial r in [0,1]) x (boundary S^3). The radial reduction A296-A318 dropped
the boundary harmonic structure (worked at l = 0); the pure-boundary reduction
A319 dropped the radial bound state. This probe builds the ONE arena not yet
tested: the COUPLED bulk+boundary problem, keeping BOTH sectors.

THE SEVEN-ARENA PATTERN (the central finding the probe is built against).
Across seven arenas the mu entry of the charged-lepton selection rule splits in
two. The channels that do NOT distinguish the three families UNDERSHOOT:
  pure-radial mu = 1.222, r^3 measure mu = 1.337.
The channels that DO give three genuinely distinct families OVERSHOOT:
  fiber m^2 energy 2.3, Dirac first-order 3.14, boundary-lens 3.86,
  boundary-harm 6.28.
The target mu = 1.5 sits in a GAP [1.34, 2.0] reached by NO single arena: neither
the radial reduction (A300-A318) nor the pure-boundary reduction (A319). The
hypothesis: the radial sector undershoots, the boundary sector overshoots, so the
COUPLED contribution may interpolate mu into the gap.

THE COUPLED OBJECT (P18 S4.3, the corpus O_{l,k}, declared verbatim).
P18 S4.3 gives the explicit radial form of the master operator in the sector with
fixed boundary angular momentum l and family index k:

    O_{l,k} = -d^2/dr^2 - (3/r) d/dr + l(l+2)/r^2 + V_eff(r) + omega^k gamma,

with gamma = 3/4 (P18 S4.2, five-fold convergence, Paper 10), omega = exp(2 pi i / 3),
V_eff(r) = V_self(r) + alpha rho(r) (the canonical A311/A317 effective potential
held verbatim). THIS is the coupled object: the SAME radial operator the radial
reduction used, but now carrying BOTH the boundary-harmonic centrifugal term
l(l+2)/r^2 (the S^3 sector entering the radial equation, A319's L(L+2) Casimir made
a centrifugal term) AND the layer-cycle family term omega^k gamma (the Z_3 sector).
A state is psi_{n,L}(r) x Y_L(S^3): the radial bound state of the l = L sector
TIMES the boundary harmonic of degree L. The radial reduction set l = 0 and dropped
omega^k gamma; the boundary reduction kept L(L+2) and dropped the radial bound
state. The coupled problem keeps both.

THE COUPLED FAMILY MAPS (a small declared set, NO fitting).

  (C1) CANONICAL. Family n -> boundary harmonic degree L_n = n - 1 coupled to the
    radial GROUND state of the l = L_n sector of O_{l,k} (with the family phase
    k = n - 1, so omega^k gamma enters too). This is the corpus O_{l,k} structure
    declared up front: the boundary degree and the family index are BOTH n - 1.
    The full B^4 overlap is taken over the genuine measure r^3 dr (radial volume)
    times the S^3 harmonic normalization, the multiplicity-floored (L+1)^2 factor
    that the boundary harmonic carries (A319's boundary spatial-overlap analog).
    The degree-n density coupling is the bulk density moment
    <psi | C[n-1] pi^n r^n | psi>_{r^3} (the SAME radial moment as every prior
    probe) TIMES the boundary harmonic factor (the S^3 piece), divided by the FULL
    coupled eigenvalue lambda_{n,L} = (radial eigenvalue of the l=L_n sector) +
    L_n(L_n+2) (the boundary Casimir added to the radial spectrum, the genuine
    coupled spectrum). Both the radial moment AND the boundary Casimir enter; this
    is the genuine coupled object. (This is the CANONICAL map; a HIT requires it.)

  (C2) LINKED. Family n -> the coupled mode where the radial node index and the
    boundary degree are LINKED as O_{l,k} specifies the family/angular pairing:
    the family index k = n - 1 selects BOTH the family phase omega^k gamma AND the
    radial NODE index (the n-th radial mode, node count n - 1) within the l = L_n
    sector. So family n is the (n-1)-th excited radial state of the l = n-1 sector,
    not the ground state. The overlap and eigenvalue are both coupled exactly as in
    C1 (full B^4 measure, lambda_{n,L} = radial + L(L+2)), but the radial state is
    the linked excited mode. This is the alternative coupled linkage stated up front.

  (C3) CONTROLS (the two limits the coupling should interpolate between).
    (C3a) PURE RADIAL: L = 0 for all families (no boundary harmonic, no omega^k
      gamma phase distinction beyond the real part), recovers the A311 pure-radial
      mu = 1.222 -- the undershoot limit.
    (C3b) PURE BOUNDARY: radial ground fixed (the l = 0 ground state for all
      families), the boundary Casimir L(L+2) and multiplicity (L+1)^2 carrying the
      family distinction, recovers the A319 boundary overshoot -- the overshoot
      limit. Both controls must FAIL the valid-peak test (degenerate / no genuine
      coupled structure).

THE KEY QUESTION. Does any coupled map land mu in the GAP [1.34, 2.0] -- the band
NEITHER the pure radial (undershoots, <= 1.34) NOR the pure boundary (overshoots,
>= 2.0) could reach -- or hit (1, 1.5, 8) with a valid peak? If the coupling
interpolates the undershooting radial and the overshooting boundary into the gap,
the coupled problem is the right arena.

PRE-REGISTERED VERDICT (declared before any number is computed).
  HIT / FORCED (CLOSES the magnitude piece of OI-287-1) iff the CANONICAL coupled
    map C1 (family n -> boundary degree L = n-1 coupled to its radial ground state,
    the corpus O_{l,k} structure, declared up front) reproduces (1, 1.5, 8) within
    20% on all three entries AND restores a valid column-standardized peak, full
    B^4 measure, NO fitting.
  STRENGTHENED-PARTIAL iff a coupled map reaches the GAP mu in [1.34, 2.0] (which
    NO single reduction could) and/or gets L2 below the standing best 0.399, but
    does not fully hit / peak. This would confirm the coupled problem is the right
    arena.
  NULL iff the coupled problem does not reach the gap and does not improve on
    L2 = 0.399 -- then the magnitude is in NONE of radial, boundary, or their
    coupling, and the framing itself is the likely issue (companion probe A321).

CRITICAL HONESTY: the canonical coupled map is the corpus O_{l,k} (family n ->
degree L = n-1 coupled to its radial ground, declared before computing). A HIT
requires THAT map (C1). A tuned (radial, boundary) pairing is a FIT and is flagged.
tau and peak are reported, not just mu. The probe states plainly whether the gap
[1.34, 2.0] is reached. No NaN; the coupled eigenvalue lambda_{n,L} > 0. N = 1000
radial; the boundary Casimir L(L+2) is exact (so(4) Casimir, A314/A319), no S^3
discretization.
"""
import os
import json
import math

import numpy as np


def gen_eigh(A, Bm):
    """Generalized symmetric-definite eig A v = lambda W v, W SPD, via numpy
    Cholesky reduction (scipy unavailable). Verbatim from A311/A317.
        W = L L^T ; C = L^{-1} A L^{-T} ; C y = lambda y ; v = L^{-T} y.
    """
    L = np.linalg.cholesky(Bm)
    Linv = np.linalg.inv(L)
    Cmat = Linv @ A @ Linv.T
    Cmat = 0.5 * (Cmat + Cmat.T)
    w, y = np.linalg.eigh(Cmat)
    v = Linv.T @ y
    return w, v


PI = math.pi
OMEGA = 4 * PI**3 + PI**2 + PI            # alpha^{-1}
E_SELF = 13.177
GAMMA = 0.75                              # P18 S4.2: gamma = 3/4 (five-fold conv.)
C = (2.0, 3.0, 16.0)                      # raw target coefficients
TARGET = [c / C[0] for c in C]           # normalized: (1, 1.5, 8)
N = 1000

# the standing bar (OPEN_HEARTS_STATUS, A311/A316): node-count map closest L2
RADIAL_BEST_L2 = 0.399
RADIAL_BEST_MU = 1.222
# the GAP no single reduction reached (frontier after A319)
GAP_LO, GAP_HI = 1.34, 2.0

# CANONICAL map: family n -> boundary degree L = n-1, family index k = n-1.
FAMILY_L = [0, 1, 2]                      # L_n = n - 1
FAMILY_K = [0, 1, 2]                      # k_n = n - 1 (omega^k phase)

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


# --- corpus density / potential (verbatim from A311/A317) --------------------
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 grid():
    dx = 1.0 / N
    x = np.linspace(dx, 1 - dx, N - 1)
    return x, dx


def radial4d_matrices(x, dx, V_on_grid, weight):
    """Symmetric 4D-radial operator A and mass W, self-adjoint in weight dr.
    Verbatim from A313/A316/A317 (radial4d_matrices). V_on_grid carries the FULL
    coupled potential V_eff(r) + l(l+2)/r^2 (the boundary centrifugal term)."""
    n = N - 1
    w_mid = 0.5 * (weight[:-1] + weight[1:])
    w_mid_left0 = max(0.5 * (0.0 + weight[0]), 0.0)
    A = np.zeros((n, n))
    for i in range(n):
        wl = w_mid[i - 1] if i > 0 else w_mid_left0
        wr = w_mid[i] if i + 1 < n else 0.5 * (weight[i] + weight[i])
        A[i, i] = (wl + wr) / dx**2
        if i + 1 < n:
            A[i, i + 1] = -w_mid[i] / dx**2
            A[i + 1, i] = -w_mid[i] / dx**2
    A += np.diag(weight * V_on_grid)
    A = 0.5 * (A + A.T)
    return A, np.diag(weight)


def solve_coupled_sector(L, n_modes=3):
    """Solve the COUPLED sector operator O_{l,k} (P18 S4.3) in r^3 dr:

        O_L = -d^2/dr^2 - (3/r) d/dr + L(L+2)/r^2 + V_eff(r)

    self-adjoint in the r^3 dr measure. The boundary harmonic degree L enters the
    radial equation through the centrifugal l(l+2)/r^2 term (P18 S4.3). The
    omega^k gamma term is a constant shift on each sector (family phase); it is
    added to the eigenvalue at coupling time (Re(omega^k) gamma), not to the
    radial operator that determines the eigenSTATE (a constant shift commutes with
    the operator and does not change the eigenfunctions, only the eigenvalues).

    Returns x, dx, radial eigenvalues, r^3-normalized eigenvectors (n_modes lowest).
    """
    x, dx = grid()
    w_arr = x**3
    centrifugal = L * (L + 2) / x**2 if L > 0 else np.zeros_like(x)
    Vfull = Vpot(x) + centrifugal
    A, Wm = radial4d_matrices(x, dx, Vfull, w_arr)
    w, v = gen_eigh(A, Wm)
    order = np.argsort(w)
    w = w[order]
    v = v[:, order]
    vv = v[:, :n_modes].copy()
    for k in range(n_modes):
        nrm = math.sqrt(float(np.sum(vv[:, k]**2 * w_arr) * dx))
        vv[:, k] = vv[:, k] / nrm
    return x, dx, w[:n_modes], vv


def boundary_harmonic_factor(L):
    """The S^3 piece of the full B^4 overlap: the boundary harmonic normalization
    factor the degree-L harmonic carries, the multiplicity (L+1)^2 (A319's
    boundary spatial-overlap analog; the harmonic "volume" the density couples to,
    multiplicity-floored so L=0 carries (0+1)^2 = 1, no zero). Genuine A314 datum."""
    return float((L + 1)**2)


def coupled_eigenvalue(radial_ev, L, k):
    """The FULL coupled eigenvalue lambda_{n,L} = (radial eigenvalue of the l=L
    sector) + L(L+2) (boundary Casimir added to the radial spectrum) + the family
    phase Re(omega^k) gamma (the layer-cycle shift, P18 S4.3 omega^k gamma). Both
    the radial spectrum AND the boundary Casimir enter; the genuine coupled
    spectrum. Re(omega^k) = cos(2 pi k / 3) = (1, -1/2, -1/2)."""
    casimir = L * (L + 2)
    phase = math.cos(2 * PI * k / 3.0) * GAMMA
    lam = radial_ev + casimir + phase
    return lam


def overlap_row(col, x, dx, w_arr, lam, harm_factor):
    """The full B^4 degree-k density overlap row for one coupled mode:
        R[n][k] = <psi | C[k-1] pi^k r^k | psi>_{r^3} * harm_factor / lam
    the bulk density moment in r^3 dr (the radial moment, every prior probe) TIMES
    the boundary harmonic factor (the S^3 piece) divided by the FULL coupled
    eigenvalue lam = radial + L(L+2) + Re(omega^k) gamma. Both sectors enter."""
    norm = float(np.sum(col**2 * w_arr) * dx)
    return [float(np.sum(col**2 * C[k - 1] * PI**k * x**k * w_arr) * dx)
            / norm * harm_factor / lam for k in (1, 2, 3)]


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)]
    if abs(d[0]) < 1e-300:
        return [float("nan")] * 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 box_states(x):
    psi_box = np.column_stack([math.sqrt(2.0) * np.sin((n + 1) * PI * x)
                               for n in range(3)])
    lam_box = np.array([(n + 1)**2 * PI**2 for n in range(3)])
    return psi_box, lam_box


# --- C1 CANONICAL: family n -> degree L=n-1 ground state, coupled -------------
def build_C1(x, dx, w_arr):
    """Family n -> boundary degree L=n-1 coupled to the radial GROUND state of
    the l=L sector. Full B^4 overlap: radial moment * (L+1)^2 / (radial + L(L+2)
    + Re(omega^k) gamma)."""
    states = np.zeros((N - 1, 3))
    radial_evs = np.zeros(3)
    lams = np.zeros(3)
    harm = np.zeros(3)
    nodes = [0, 0, 0]
    for n in range(3):
        L = FAMILY_L[n]
        k = FAMILY_K[n]
        _, _, w_L, v_L = solve_coupled_sector(L, n_modes=3)
        states[:, n] = v_L[:, 0]            # GROUND mode of the l=L sector
        radial_evs[n] = w_L[0]
        lams[n] = coupled_eigenvalue(w_L[0], L, k)
        harm[n] = boundary_harmonic_factor(L)
        nodes[n] = node_count(v_L[:, 0])
    M = np.zeros((3, 3))
    for n in range(3):
        M[n, :] = overlap_row(states[:, n], x, dx, w_arr, lams[n], harm[n])
    return M, radial_evs, lams, harm, nodes


# --- C2 LINKED: family n -> (n-1)-th excited radial mode of l=n-1 sector ------
def build_C2(x, dx, w_arr):
    """Family n -> the LINKED coupled mode: the radial node index and boundary
    degree are linked, family n is the (n-1)-th radial mode (node count n-1) of
    the l=L=n-1 sector. Overlap and eigenvalue coupled as C1."""
    states = np.zeros((N - 1, 3))
    radial_evs = np.zeros(3)
    lams = np.zeros(3)
    harm = np.zeros(3)
    nodes = [0, 0, 0]
    for n in range(3):
        L = FAMILY_L[n]
        k = FAMILY_K[n]
        mode_idx = n                        # (n-1)-th excited = index n-1; n=0->0
        _, _, w_L, v_L = solve_coupled_sector(L, n_modes=3)
        states[:, n] = v_L[:, mode_idx]     # LINKED excited mode
        radial_evs[n] = w_L[mode_idx]
        lams[n] = coupled_eigenvalue(w_L[mode_idx], L, k)
        harm[n] = boundary_harmonic_factor(L)
        nodes[n] = node_count(v_L[:, mode_idx])
    M = np.zeros((3, 3))
    for n in range(3):
        M[n, :] = overlap_row(states[:, n], x, dx, w_arr, lams[n], harm[n])
    return M, radial_evs, lams, harm, nodes


# --- C3a PURE RADIAL control: L=0 all families (recovers A311 mu=1.222) -------
def build_C3a(x, dx, w_arr):
    """Control: L = 0 for all three families, the three lowest radial modes of the
    l=0 sector, no boundary harmonic factor (harm = 1), no coupled Casimir, no
    omega^k phase distinction (k effectively folded). Recovers the A311 pure-radial
    diagonal mu = 1.222. Must FAIL the peak test (the undershoot limit)."""
    _, _, w0, v0 = solve_coupled_sector(0, n_modes=3)
    M = np.zeros((3, 3))
    radial_evs = w0.copy()
    nodes = [node_count(v0[:, n]) for n in range(3)]
    for n in range(3):
        # pure radial: lam = radial eigenvalue only, harm = 1 (no S^3 piece)
        M[n, :] = overlap_row(v0[:, n], x, dx, w_arr, w0[n], 1.0)
    return M, radial_evs, nodes


# --- C3b PURE BOUNDARY control: radial ground fixed (recovers A319 overshoot) -
def build_C3b(x, dx, w_arr):
    """Control: the l=0 radial GROUND state fixed for all families, the boundary
    Casimir L(L+2) and multiplicity (L+1)^2 carrying the family distinction (the
    A319 pure-boundary reduction, no radial bound-state variation). Must overshoot
    and FAIL the peak test (the overshoot limit)."""
    _, _, w0, v0 = solve_coupled_sector(0, n_modes=1)
    ground = v0[:, 0]
    M = np.zeros((3, 3))
    lams = np.zeros(3)
    harm = np.zeros(3)
    for n in range(3):
        L = FAMILY_L[n]
        cas = L * (L + 2)
        w = cas if cas > 1e-9 else 1.0      # L=0 multiplicity floor (A319)
        harm[n] = boundary_harmonic_factor(L)
        lams[n] = w
        # boundary reduction: radial moment of the SAME ground state, boundary
        # spectral weight + multiplicity carry the distinction
        M[n, :] = overlap_row(ground, x, dx, w_arr, w, harm[n])
    return M, lams, harm


def evaluate(M, M_ctrl):
    ratio = norm_diag(M)
    peak = bool(diag_peaks(M))
    ctrl_peak = bool(diag_peaks(M_ctrl))
    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-9
                and abs(diag_entries[0] - diag_entries[2]) > 1e-9
                and abs(diag_entries[1] - diag_entries[2]) > 1e-9)
    L2v = l2(ratio)
    return {
        "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,
        "control_peaks": ctrl_peak,
        "valid_peak": bool(peak and not ctrl_peak),
        "within_20pct": bool(within20),
        "hits_target": bool(within20 and peak and not ctrl_peak),
        "L2_to_target": round(float(L2v), 4),
        "three_distinct_entries": bool(distinct),
        "mu_in_gap": bool(GAP_LO <= ratio[1] <= GAP_HI),
        "beats_radial_L2": bool(L2v < RADIAL_BEST_L2 - 1e-9),
    }


def main():
    x, dx = grid()
    w_arr = x**3
    psi_box, lam_box = box_states(x)

    # the two controls' raw matrices (used as the peak-control for each map)
    M_c3a, evs_c3a, nodes_c3a = build_C3a(x, dx, w_arr)
    M_c3b, lams_c3b, harm_c3b = build_C3b(x, dx, w_arr)

    # the coupled maps
    M_c1, evs_c1, lams_c1, harm_c1, nodes_c1 = build_C1(x, dx, w_arr)
    M_c2, evs_c2, lams_c2, harm_c2, nodes_c2 = build_C2(x, dx, w_arr)

    # box control matrix for the peak test (same r^3 measure)
    M_box = np.zeros((3, 3))
    for n in range(3):
        M_box[n, :] = overlap_row(psi_box[:, n], x, dx, w_arr, lam_box[n], 1.0)

    # evaluate each map against the box control peak
    res_c1 = evaluate(M_c1, M_box)
    res_c2 = evaluate(M_c2, M_box)
    res_c3a = evaluate(M_c3a, M_box)
    res_c3b = evaluate(M_c3b, M_box)

    # canonical = C1 (corpus O_{l,k}, family n -> degree L=n-1 ground, declared)
    canonical = res_c1

    # NaN guard
    nan_found = False
    for r in (res_c1, res_c2, res_c3a, res_c3b):
        for v in r["diagonal_ratio"] + [r["mu"], r["tau"], r["L2_to_target"]]:
            if v != v:
                nan_found = True
    # coupled eigenvalues must be > 0
    lam_positive = all(float(v) > 0 for v in list(lams_c1) + list(lams_c2))

    # does ANY coupled map (C1 or C2) reach the gap [1.34, 2.0]?
    c1_gap = res_c1["mu_in_gap"]
    c2_gap = res_c2["mu_in_gap"]
    reaches_gap = bool(c1_gap or c2_gap)
    # does any coupled map beat the standing-best L2?
    beats_radial = bool(res_c1["beats_radial_L2"] or res_c2["beats_radial_L2"])

    # control sanity: C3a recovers the radial undershoot (mu ~ 1.22), C3b
    # overshoots (mu >= 2.0). Recorded as a sanity flag, not a hard gate here
    # (the verifier asserts the limits exactly).
    c3a_undershoots = res_c3a["mu"] <= GAP_LO + 0.1
    c3b_overshoots = res_c3b["mu"] >= GAP_HI - 0.1

    # ---- VERDICT (canonical = C1, corpus O_{l,k}) --------------------------
    if canonical["hits_target"]:
        verdict = "FORCED"
    elif reaches_gap or beats_radial:
        verdict = "STRENGTHENED_PARTIAL"
    else:
        verdict = "NULL"

    summary = {
        "addendum": "A320",
        "addresses": "OI-287-1",
        "builds_on": ["A311", "A317", "A319"],
        "thesis": ("the COUPLED bulk+boundary B^4 problem (P18 S4.3 O_{l,k}: the "
                   "radial operator with BOTH the boundary-harmonic centrifugal "
                   "term l(l+2)/r^2 AND the layer-cycle family term omega^k gamma "
                   "together). The radial reduction (A300-A318) undershoots, the "
                   "pure-boundary reduction (A319) overshoots; KEY QUESTION: does "
                   "the coupled contribution interpolate mu into the GAP [1.34, "
                   "2.0] no single reduction could reach, or hit (1,1.5,8) with a "
                   "peak?"),
        "target_normalized": TARGET,
        "constants": {"OMEGA_alpha_inv": OMEGA, "E_self": E_SELF,
                      "gamma": GAMMA, "C_raw": list(C)},
        "coupled_operator": ("O_{l,k} = -d^2/dr^2 - (3/r)d/dr + l(l+2)/r^2 + "
                             "V_eff(r) + omega^k gamma  (P18 S4.3)"),
        "family_L_assignment": FAMILY_L,
        "family_k_assignment": FAMILY_K,
        "radial_best_L2": RADIAL_BEST_L2,
        "radial_best_mu": RADIAL_BEST_MU,
        "gap_band": [GAP_LO, GAP_HI],
        # C1 CANONICAL
        "C1_radial_eigenvalues": [round(float(v), 4) for v in evs_c1],
        "C1_coupled_eigenvalues": [round(float(v), 4) for v in lams_c1],
        "C1_harmonic_factors": [round(float(v), 4) for v in harm_c1],
        "C1_node_counts": nodes_c1,
        "C1_diagonal": res_c1["diagonal_ratio"],
        "C1_mu": res_c1["mu"],
        "C1_tau": res_c1["tau"],
        "C1_L2": res_c1["L2_to_target"],
        "C1_valid_peak": res_c1["valid_peak"],
        "C1_three_distinct": res_c1["three_distinct_entries"],
        "C1_mu_in_gap": res_c1["mu_in_gap"],
        "C1_beats_radial_L2": res_c1["beats_radial_L2"],
        "C1_hits_target": res_c1["hits_target"],
        # C2 LINKED
        "C2_radial_eigenvalues": [round(float(v), 4) for v in evs_c2],
        "C2_coupled_eigenvalues": [round(float(v), 4) for v in lams_c2],
        "C2_node_counts": nodes_c2,
        "C2_diagonal": res_c2["diagonal_ratio"],
        "C2_mu": res_c2["mu"],
        "C2_tau": res_c2["tau"],
        "C2_L2": res_c2["L2_to_target"],
        "C2_valid_peak": res_c2["valid_peak"],
        "C2_three_distinct": res_c2["three_distinct_entries"],
        "C2_mu_in_gap": res_c2["mu_in_gap"],
        "C2_beats_radial_L2": res_c2["beats_radial_L2"],
        # C3 controls (the two limits)
        "C3a_pure_radial_diagonal": res_c3a["diagonal_ratio"],
        "C3a_pure_radial_mu": res_c3a["mu"],
        "C3a_pure_radial_valid_peak": res_c3a["valid_peak"],
        "C3a_recovers_radial_undershoot": bool(c3a_undershoots),
        "C3b_pure_boundary_diagonal": res_c3b["diagonal_ratio"],
        "C3b_pure_boundary_mu": res_c3b["mu"],
        "C3b_pure_boundary_valid_peak": res_c3b["valid_peak"],
        "C3b_recovers_boundary_overshoot": bool(c3b_overshoots),
        # the central questions
        "canonical_map": "C1 (family n -> degree L=n-1 ground, corpus O_{l,k})",
        "C1_reaches_gap": c1_gap,
        "C2_reaches_gap": c2_gap,
        "coupled_reaches_gap_1.34_2.0": reaches_gap,
        "coupled_beats_radial_L2": beats_radial,
        "coupled_eigenvalues_positive": bool(lam_positive),
        "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("A320 PROBE -- OI-287-1: the COUPLED bulk+boundary B^4 problem (P18 S4.3")
    print("   O_{l,k}: radial operator with BOTH l(l+2)/r^2 AND omega^k gamma). The")
    print("   radial reduction undershoots, the boundary reduction overshoots. Does")
    print("   the coupling interpolate mu into the GAP [1.34,2.0] or hit (1,1.5,8)?")
    print("=" * 74)
    print("target (normalized):", TARGET)
    print(f"standing bar: L2 = {RADIAL_BEST_L2}, radial mu = {RADIAL_BEST_MU}")
    print(f"the GAP no single reduction reached: mu in [{GAP_LO}, {GAP_HI}]")
    print(f"coupled operator: {summary['coupled_operator']}")
    print("-" * 74)
    print("[CONTROLS, the two limits]")
    print(f"  C3a pure radial (L=0 all): diag={res_c3a['diagonal_ratio']} "
          f"mu={res_c3a['mu']} peak={res_c3a['valid_peak']} "
          f"(recovers A311 undershoot: {c3a_undershoots})")
    print(f"  C3b pure boundary (ground fixed): diag={res_c3b['diagonal_ratio']} "
          f"mu={res_c3b['mu']} peak={res_c3b['valid_peak']} "
          f"(recovers A319 overshoot: {c3b_overshoots})")
    print("-" * 74)
    print("[COUPLED MAPS]")
    print(f"  C1 CANONICAL (degree L=n-1 ground, O_{{l,k}}):")
    print(f"     radial_ev={[round(float(v),3) for v in evs_c1]} "
          f"coupled_lam={[round(float(v),3) for v in lams_c1]} "
          f"harm={[round(float(v),1) for v in harm_c1]} nodes={nodes_c1}")
    print(f"     diag={res_c1['diagonal_ratio']} mu={res_c1['mu']} "
          f"tau={res_c1['tau']} peak={res_c1['valid_peak']} "
          f"L2={res_c1['L2_to_target']} mu_in_gap={res_c1['mu_in_gap']} "
          f"beats_radial={res_c1['beats_radial_L2']} HIT={res_c1['hits_target']}")
    print(f"  C2 LINKED (excited mode n-1 of l=n-1):")
    print(f"     coupled_lam={[round(float(v),3) for v in lams_c2]} "
          f"nodes={nodes_c2}")
    print(f"     diag={res_c2['diagonal_ratio']} mu={res_c2['mu']} "
          f"tau={res_c2['tau']} peak={res_c2['valid_peak']} "
          f"L2={res_c2['L2_to_target']} mu_in_gap={res_c2['mu_in_gap']} "
          f"beats_radial={res_c2['beats_radial_L2']}")
    print("-" * 74)
    print(f"coupled reaches the [1.34,2.0] GAP no single reduction could: "
          f"{reaches_gap} (C1 {c1_gap}, C2 {c2_gap})")
    print(f"coupled beats the standing-best L2 (<0.399): {beats_radial}")
    print(f"coupled eigenvalues all > 0: {lam_positive}")
    print(f"NaN found: {nan_found}")
    print(f"\nVERDICT: {verdict}")
    print("-" * 74)
    if verdict == "FORCED":
        print("  -> the CANONICAL coupled map C1 (corpus O_{l,k}, family n ->")
        print("     degree L=n-1 ground) reproduces (1,1.5,8) within 20% with a")
        print("     valid peak, full B^4 measure, NO fitting. The coupled problem")
        print("     interpolates the radial undershoot and boundary overshoot into")
        print("     the target. OI-287-1 magnitude CLOSES.")
    elif verdict == "STRENGTHENED_PARTIAL":
        print("  -> a coupled map REACHES the [1.34,2.0] gap no single reduction")
        print("     could, and/or beats the standing-best L2, short of a full")
        print("     hit/peak. The COUPLED problem is the right arena: the magnitude")
        print("     lives in the interpolation, not in either reduction alone.")
    else:
        print("  -> the coupled problem does NOT reach the gap and does NOT beat")
        print("     the standing-best L2. The magnitude is in NONE of radial,")
        print("     boundary, or their coupling. The framing itself is the likely")
        print("     issue: (2,3,16) may not be an overlap-diagonal quantity at all")
        print("     (companion probe A321).")
    print(f"wrote {OUT}")


if __name__ == "__main__":
    main()
