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

DERIVATION probe, not a fit. A311 established that the inner-product measure
that carries the A307 lift is GEOMETRICALLY FORCED: it is the radial volume
element r^3 dr of the B^4 arena (P18 S2.1 base metric ds^2 = dr^2 + r^2
dOmega_3^2), and the density rho is r^3 to leading order (0.933 of rho(1)). A311
numbers: the flat-dx control recovers the A299 base (1, 1.12, 7.50); the r^3
measure on the canonical flat-1D operator gives mu=1.337, diagonal
(1, 1.337, 8.706), L2=0.724; the corpus-faithful 4D-radial operator
-d^2/dr^2 - (3/r) d/dr + V, self-adjoint in its r^3 dr measure, gives mu=1.222,
diagonal (1, 1.222, 8.286), L2=0.399 (the closest of all, sharpening tau toward
8). Two limits survived A311: mu stops short of 1.5, and no candidate restored a
column-standardized diagonal peak. A311's hypothesis: the residual is carried by
the ANGULAR (l) and FIBER (S^1) sectors that the pure-radial reduction drops.

This probe tests exactly that hypothesis. P18 S2.2 gives the bulk radial
equation with the angular-momentum term

    [-d^2/dr^2 - (3/r) d/dr + l(l+2)/r^2] R = E R   (P18 S2.2, eq 165)

and P18 S4.3 the explicit radial form of the assembled operator with the
matching l(l+2)/r^2 angular term and the fiber family term omega^k gamma:

    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, eq 367)

with omega = exp(2 pi i / 3) the cube root of unity (P18 S2.5), gamma = 3/4 the
layer-cycle coefficient (P18 S4.2), k = 0,1,2 the family index. P18 S2.2 also
gives the unperturbed bulk spectrum E_{n,l} = (2n + l + 2)^2: the radial node
index n and the angular momentum l enter the operator together. The corpus
LINKS the family/radial index and l in S4.3 (O_{l,k}); the three charged-lepton
families are the three layer-cycle eigenspaces (P18 S5.3, T psi_n^{(k)} =
omega^k psi_n^{(k)}). The geometrically forced reading of "family n" in the full
arena is therefore: family n carries angular momentum l = n (n = 0,1,2 -> the
three S^3 angular sectors with degeneracy (l+1)^2) AND fiber family index k = n
(the n-th cube-root eigenspace). This is the full-arena object A311 pointed at.

PRE-REGISTERED CANDIDATES (declared before running; FORCED by the arena vs
exploratory; no free fitting):

  CONTROLS (reproduce A311, anchor the verifier):
  (C_flat)    flat 1D operator, flat dx measure. Must recover (1, 1.12, 7.50).
              [A311 control]                                      [FORCED-anchor]
  (C_r3)      flat 1D operator, r^3 dr measure (generalized
              H psi = lambda r^3 psi). [A311: mu=1.337]                  [FORCED]
  (C_4Dself)  THE GEOMETRICALLY CONSISTENT PURE-RADIAL OBJECT: the 4D radial
              operator -d^2/dr^2 - (3/r) d/dr + V_eff, self-adjoint in r^3 dr,
              eigenstates AND overlap both in r^3 dr. [A311: mu=1.222, L2=0.399]
                                                                        [FORCED]

  FULL-ARENA (the A313 thesis -- add the dropped sectors):
  (C_ang)     ADD THE ANGULAR SECTOR. Family n is solved in the l = n angular
              sector: the radial operator gains the centrifugal term
              l(l+2)/r^2 with l = n, still self-adjoint in r^3 dr (the radial
              part of the full r^3 dr dOmega_3 volume; the S^3 angular harmonic
              of degree l is unit-normalized on S^3, so the angular integral
              contributes a fixed harmonic normalization that cancels in the
              ratio R[n][k]/<psi_n|psi_n>). Each family's mode is the ground
              radial mode of its own l-sector operator O_{l=n}. Overlap formed
              in r^3 dr.                                                 [FORCED]
  (C_fiber)   ADD THE HOPF FIBER. On top of C_ang, the declared fiber measure
              factor: the Hopf lapse m = sqrt(1 - u^2), u = cos(2 beta) the
              Hopf latitude (Paper 30 / kernel Hopf lapse). Read radially with
              the canonical chart beta = (pi/2) r so that u = cos(pi r) sweeps
              the fiber latitude across the ball and m = sqrt(1 - cos^2(pi r)) =
              |sin(pi r)|; the full-arena measure is then r^3 |sin(pi r)| dr.
              The fiber family term omega^k gamma (P18 S4.3) adds a real
              diagonal shift Re(omega^k) gamma = gamma cos(2 pi k / 3) to the
              l = n radial operator for family k = n. Eigenstates and overlap
              both in the lapse-weighted measure.                       [FORCED]

For each candidate: solve the appropriate (generalized) eigenproblem by the
numpy Cholesky reduction (scipy unavailable), form the A299-style diagonal
overlap IN THAT MEASURE, report normalized diagonal, valid_peak (column-
standardized argmax over rows of column k == k for all k AND the box control
psi_box[n] = sqrt(2) sin((n+1) pi x), formed in the SAME measure, FAILS the same
test), L2 to (1, 1.5, 8), mu = ratio[1], tau = ratio[2].

For the angular/fiber candidates the three families are no longer the lowest
three modes of one operator; they are the GROUND radial mode of three different
l-sector operators (l = 0,1,2), since "family n carries l = n". The diagonal
entry n is then <phi_n | c_k pi^k r^k | phi_n>_W / lambda_n with phi_n the ground
mode of O_{l=n} (and, for C_fiber, the omega^k gamma shift on the family
operator). This is the corpus' explicit S4.3 O_{l,k} structure, not a free
choice.

PRE-REGISTERED VERDICT (declared before running):
  HIT iff some geometrically-forced full-arena candidate (C_ang or C_fiber, no
       fitting) yields BOTH valid_peak AND all three normalized entries within
       20% of (1, 1.5, 8). This would CLOSE the magnitude piece of OI-287-1.
  STRENGTHENED-PARTIAL iff adding the angular and/or fiber sector moves mu
       materially closer to 1.5 than the pure-radial value (mu > 1.40) OR
       restores a valid peak without full magnitude. Report which sector
       (angular vs fiber) supplies the motion.
  NULL iff the full-arena measure does not improve on the pure-radial A311
       result (mu stays ~1.33, no peak). Then the residual is NOT in the dropped
       angular/fiber sectors and the gap to 1.5 is elsewhere; recorded as the
       finding (it redirects the search).

Multiple-comparison caveat: 5 candidates declared (all geometrically FORCED or
FORCED-anchor controls; no exploratory free knobs), counted. Honesty over a hit.
W > 0 on (0,1) for every generalized problem; no NaN; angular/fiber assignments
read from P18 S2.2 / S4.3, not fitted.
"""
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 in the sandbox).
        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 layer-cycle coefficient 3/4
C = (2.0, 3.0, 16.0)                     # raw target coefficients
TARGET = [c / C[0] for c in C]           # normalized: (1, 1.5, 8)
N = 1500

HERE = os.path.dirname(os.path.abspath(__file__))
OUT = os.path.join(HERE, "..", "..", "..", "..", "Outputs",
                   "selection_rule_fullarena.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 grid():
    dx = 1.0 / N
    x = np.linspace(dx, 1 - dx, N - 1)
    return x, dx


def hamiltonian(x, dx, potential):
    """Canonical flat-1D operator H = -d^2/dx^2 + V, Dirichlet, A299/A307 form."""
    off = -1.0 / dx**2 * np.ones(N - 2)
    return (np.diag(2.0 / dx**2 + potential(x))
            + np.diag(off, 1) + np.diag(off, -1))


def eigfuns_flat(potential):
    """Standard eigenproblem H psi = lambda psi, flat-dx normalized."""
    x, dx = grid()
    H = hamiltonian(x, dx, potential)
    w, v = np.linalg.eigh(H)
    return x, dx, w[:3], v[:, :3] / math.sqrt(dx)


def eigfuns_weighted(potential, weight):
    """Generalized eigenproblem H psi = lambda W psi, W = diag(weight(x)),
    H the canonical flat-1D operator. psi W-normalized."""
    x, dx = grid()
    H = hamiltonian(x, dx, potential)
    Wdiag = weight(x)
    if np.any(Wdiag <= 0):
        raise ValueError("weight must be strictly positive on the interior")
    w, v = gen_eigh(H, np.diag(Wdiag))
    order = np.argsort(w)
    w = w[order]
    v = v[:, order]
    vv = v[:, :3].copy()
    for n in range(3):
        nrm = math.sqrt(float(np.sum(vv[:, n]**2 * Wdiag) * dx))
        vv[:, n] = vv[:, n] / nrm
    return x, dx, w[:3], vv


def radial4d_matrices(x, dx, V_on_grid, weight):
    """Assemble the symmetric 4D-radial operator A and mass W for the weight
    W = weight (the radial part of the arena measure on the grid):

        A_ij = sum W R_i' R_j' + W V R_i R_j ,  W_ij = diag(W),

    self-adjoint in the W dr measure (Sturm-Liouville -1/W d/dr(W d/dr) + V).
    For weight = r^3 this is P18 S2.2's -d^2/dr^2 - (3/r) d/dr + V; a general
    positive weight keeps the operator self-adjoint in W dr.

    Midpoint weights keep the discrete operator symmetric.
    """
    n = N - 1
    w_mid = 0.5 * (weight[:-1] + weight[1:])           # weight at interior mids
    w_mid_left0 = max(0.5 * (0.0 + weight[0]), 0.0)    # mid between r=0 and x[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)
    Wm = np.diag(weight)
    return A, Wm


def eigfuns_radial4d(potential, weight=None):
    """The corpus-faithful 4D-radial operator (P18 S2.2) in its r^3 dr measure
    (or a general positive weight for the fiber-lapse measure). Returns the
    lowest three W-orthonormal radial modes and their eigenvalues."""
    x, dx = grid()
    w_arr = x**3 if weight is None else weight(x)
    V = potential(x)
    A, Wm = radial4d_matrices(x, dx, V, w_arr)
    w, v = gen_eigh(A, Wm)
    order = np.argsort(w)
    w = w[order]
    v = v[:, order]
    vv = v[:, :3].copy()
    for k in range(3):
        nrm = math.sqrt(float(np.sum(vv[:, k]**2 * w_arr) * dx))
        vv[:, k] = vv[:, k] / nrm
    return x, dx, w[:3], vv


def eigfuns_angular(potential, weight, l_shift_fn):
    """FULL-ARENA angular candidate: family n is the GROUND radial mode of the
    l = n sector operator O_{l=n} = -d^2/dr^2 - (3/r) d/dr + l(l+2)/r^2 + V(r) +
    shift(n), self-adjoint in the weight dr measure.

    l_shift_fn(n) returns an additive constant diagonal shift for family n (used
    by the fiber candidate for the omega^k gamma term; 0.0 for the pure angular
    candidate). The centrifugal term l(l+2)/r^2 with l = n is folded into the
    potential of the n-th sector. Returns the three family ground modes
    (one per l-sector) and their ground eigenvalues, all W-orthonormal in
    weight dr.
    """
    x, dx = grid()
    w_arr = weight(x)
    if np.any(w_arr <= 0):
        raise ValueError("weight must be strictly positive on the interior")
    psis = np.zeros((N - 1, 3))
    lams = np.zeros(3)
    for n in range(3):
        l = n
        cent = l * (l + 2) / x**2                       # centrifugal term
        V_n = potential(x) + cent + l_shift_fn(n)
        A, Wm = radial4d_matrices(x, dx, V_n, w_arr)
        w, v = gen_eigh(A, Wm)
        order = np.argsort(w)
        ground = v[:, order[0]]
        nrm = math.sqrt(float(np.sum(ground**2 * w_arr) * dx))
        psis[:, n] = ground / nrm
        lams[n] = w[order[0]]
    return x, dx, lams, psis


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


def density_matrix(psi, x, dx, measure):
    """R[n][k] = <psi_n | c_k pi^k x^k | psi_n>_W / <psi_n|psi_n>_W."""
    M = np.zeros((3, 3))
    for n in range(3):
        norm = float(np.sum(psi[:, n]**2 * measure) * dx)
        for k in (1, 2, 3):
            num = float(np.sum(psi[:, n]**2 * C[k - 1] * PI**k * x**k
                               * measure) * dx)
            M[n, k - 1] = num / norm
    return M


def diag_peaks(M):
    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 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


def evaluate(psi, lam, x, dx, psi_box, lam_box, measure):
    """Diagonal / peak / L2 / mu for one measure W = measure(array)."""
    Dm = density_matrix(psi, x, dx, measure)
    B = Dm / np.array(lam)[:, None]
    Dm_box = density_matrix(psi_box, x, dx, measure)
    B_box = Dm_box / lam_box[:, None]
    ratio = norm_diag(B)
    peak = bool(diag_peaks(B))
    box_peak = bool(diag_peaks(B_box))
    valid_peak = peak and not box_peak
    dist = l2(ratio)
    within20 = all(abs(ratio[i] - TARGET[i]) <= 0.20 * TARGET[i]
                   for i in range(3))
    mu = ratio[1]
    return {
        "diagonal_ratio": [round(float(v), 4) for v in ratio],
        "tau": round(float(ratio[2]), 4),
        "peaks_diagonal": peak,
        "box_control_peaks": box_peak,
        "valid_peak": valid_peak,
        "L2_to_target": round(dist, 4),
        "within_20pct": bool(within20),
        "mu": round(float(mu), 4),
        "eigenvalues": [round(float(v), 4) for v in lam],
        "node_counts": [node_count(psi[:, n]) for n in range(3)],
    }


# All candidates are geometrically forced (controls anchor A311; the full-arena
# angular/fiber candidates read their l and omega^k gamma from P18 S2.2/S4.3).
FORCED = {"C_flat", "C_r3", "C_4Dself", "C_ang", "C_fiber"}
FULL_ARENA = {"C_ang", "C_fiber"}
PURE_RADIAL = {"C_r3", "C_4Dself"}


def main():
    x, dx = grid()
    _, _, lam_flat, psi_flat = eigfuns_flat(Vpot)
    psi_box, lam_box = box_states(x)
    ones = np.ones_like(x)
    r3 = x**3
    # Hopf lapse measure: r^3 |sin(pi r)| (the radial sweep of the lapse
    # m = sqrt(1 - cos^2(pi r)) with beta = (pi/2) r the Hopf latitude chart)
    lapse = np.abs(np.sin(PI * x))
    fiber_w = r3 * lapse

    results = {}

    # (C_flat) control: W = 1, canonical H -> must recover A299 (1,1.12,7.50)
    results["C_flat"] = evaluate(psi_flat, lam_flat, x, dx,
                                 psi_box, lam_box, ones)

    # (C_r3) FORCED: W = r^3, canonical flat-1D H, generalized problem [A311]
    _, _, lam_r3, psi_r3 = eigfuns_weighted(Vpot, lambda xx: xx**3)
    results["C_r3"] = evaluate(psi_r3, lam_r3, x, dx, psi_box, lam_box, r3)

    # (C_4Dself) FORCED: the corpus-faithful 4D-radial operator in r^3 [A311]
    _, _, lam_4d, psi_4d = eigfuns_radial4d(Vpot)
    results["C_4Dself"] = evaluate(psi_4d, lam_4d, x, dx,
                                   psi_box, lam_box, r3)

    # (C_ang) FULL-ARENA: add the angular sector, family n -> l = n
    _, _, lam_ang, psi_ang = eigfuns_angular(
        Vpot, lambda xx: xx**3, lambda n: 0.0)
    results["C_ang"] = evaluate(psi_ang, lam_ang, x, dx,
                                psi_box, lam_box, r3)

    # (C_fiber) FULL-ARENA: add the Hopf fiber -- lapse measure r^3|sin(pi r)|
    #   plus the omega^k gamma family shift Re(omega^n) gamma on the l=n sector
    def fiber_shift(n):
        return GAMMA * math.cos(2 * PI * n / 3.0)
    _, _, lam_fib, psi_fib = eigfuns_angular(
        Vpot, lambda xx: xx**3 * np.abs(np.sin(PI * xx)), fiber_shift)
    results["C_fiber"] = evaluate(psi_fib, lam_fib, x, dx,
                                  psi_box, lam_box, fiber_w)

    # ---- aggregate -------------------------------------------------------
    flat_mu = float(results["C_flat"]["mu"])
    r3_mu = float(results["C_r3"]["mu"])
    self_mu = float(results["C_4Dself"]["mu"])
    pure_radial_best_mu = max(r3_mu, self_mu, key=lambda m: -abs(m - 1.5))
    # closest pure-radial mu to 1.5 (for the "materially closer" comparison)
    pr_closest_to_15 = min((r3_mu, self_mu), key=lambda m: abs(m - 1.5))

    non_control = {k: v for k, v in results.items() if k != "C_flat"}
    closest = min(non_control, key=lambda k: non_control[k]["L2_to_target"])
    closest_ratio = results[closest]["diagonal_ratio"]
    closest_mu = results[closest]["mu"]
    closest_l2 = results[closest]["L2_to_target"]

    for k, v in results.items():
        v["mu_moved_toward_1.5"] = bool(
            abs(v["mu"] - 1.5) < abs(flat_mu - 1.5))
        v["is_full_arena"] = k in FULL_ARENA
        v["is_forced"] = k in FORCED

    any_valid = any(v["valid_peak"] for v in results.values())
    full_arena_valid = any(results[k]["valid_peak"] for k in FULL_ARENA)
    hit = any(results[k]["valid_peak"] and results[k]["within_20pct"]
              for k in FULL_ARENA)

    # which full-arena candidate moves mu closest to 1.5, and by how much
    fa_mus = {k: float(results[k]["mu"]) for k in FULL_ARENA}
    fa_best = min(fa_mus, key=lambda k: abs(fa_mus[k] - 1.5))
    fa_best_mu = fa_mus[fa_best]

    # STRENGTHENED-PARTIAL: a full-arena candidate gives mu > 1.40 OR restores
    # a valid peak (materially closer than the pure-radial value)
    fa_mu_gt_140 = any(fa_mus[k] > 1.40 for k in FULL_ARENA)
    moved_closer_than_radial = any(
        abs(fa_mus[k] - 1.5) < abs(pr_closest_to_15 - 1.5) - 1e-9
        for k in FULL_ARENA)
    strengthened_partial = (fa_mu_gt_140 or full_arena_valid)

    # which sector supplies the motion (angular vs fiber)
    ang_mu = float(results["C_ang"]["mu"])
    fib_mu = float(results["C_fiber"]["mu"])
    # the angular sector's contribution: C_ang mu vs the pure-radial C_r3 mu
    angular_motion = ang_mu - r3_mu
    # the fiber sector's contribution: C_fiber mu vs C_ang mu
    fiber_motion = fib_mu - ang_mu
    if abs(angular_motion) >= abs(fiber_motion):
        sector_supplying_motion = "angular"
    else:
        sector_supplying_motion = "fiber"

    # NULL: the full-arena measure does not improve on the pure-radial result
    # (mu stays ~1.33, no peak, no material move toward 1.5)
    full_arena_no_improve = (
        (not full_arena_valid)
        and (not fa_mu_gt_140)
        and (not moved_closer_than_radial))

    if hit:
        verdict = "HIT"
    elif strengthened_partial:
        verdict = "STRENGTHENED_PARTIAL"
    elif full_arena_no_improve:
        verdict = "NULL"
    else:
        # full-arena moved mu closer to 1.5 than radial but stayed <= 1.40 and
        # no peak: a genuine but weak partial -- report as STRENGTHENED_PARTIAL
        # only on the mu>1.40-or-peak criterion above, else NULL-leaning.
        verdict = "NULL"

    nan_found = False
    for v in results.values():
        for val in v["diagonal_ratio"] + v["eigenvalues"]:
            if val != val:
                nan_found = True

    summary = {
        "addendum": "A313",
        "addresses": "OI-287-1",
        "builds_on": ["A307", "A311"],
        "thesis": ("does the full B^4 arena measure -- adding the angular l "
                   "and Hopf fiber sectors the pure-radial reduction drops -- "
                   "close mu to 1.5 and restore the diagonal peak?"),
        "target_normalized": TARGET,
        "flat_control_diagonal": results["C_flat"]["diagonal_ratio"],
        "flat_control_mu": flat_mu,
        "constants": {"OMEGA_alpha_inv": OMEGA, "E_self": E_SELF,
                      "gamma_layer_cycle": GAMMA},
        "geometric_assignment": (
            "family n carries angular momentum l = n (P18 S2.2 l(l+2)/r^2; "
            "S4.3 O_{l,k}) and fiber family index k = n (P18 S2.5 omega^k, "
            "S5.3 T psi_n^{(k)} = omega^k psi_n^{(k)}); each family is the "
            "ground radial mode of its own l-sector operator"),
        "candidates": results,
        "n_candidates_declared": len(results),
        "forced_candidates": sorted(FORCED),
        "full_arena_candidates": sorted(FULL_ARENA),
        "pure_radial_candidates": sorted(PURE_RADIAL),
        "C_flat_diagonal": results["C_flat"]["diagonal_ratio"],
        "C_r3_mu": r3_mu,
        "C_4Dself_diagonal": results["C_4Dself"]["diagonal_ratio"],
        "C_4Dself_mu": self_mu,
        "C_4Dself_L2": results["C_4Dself"]["L2_to_target"],
        "C_ang_diagonal": results["C_ang"]["diagonal_ratio"],
        "C_ang_mu": ang_mu,
        "C_ang_L2": results["C_ang"]["L2_to_target"],
        "C_fiber_diagonal": results["C_fiber"]["diagonal_ratio"],
        "C_fiber_mu": fib_mu,
        "C_fiber_L2": results["C_fiber"]["L2_to_target"],
        "pure_radial_closest_mu_to_1.5": pr_closest_to_15,
        "full_arena_best_candidate": fa_best,
        "full_arena_best_mu": fa_best_mu,
        "angular_sector_mu_motion": round(angular_motion, 4),
        "fiber_sector_mu_motion": round(fiber_motion, 4),
        "sector_supplying_motion": sector_supplying_motion,
        "full_arena_mu_exceeds_1.40": bool(fa_mu_gt_140),
        "full_arena_moved_closer_than_pure_radial": bool(
            moved_closer_than_radial),
        "full_arena_restores_valid_peak": bool(full_arena_valid),
        "closest_noncontrol": closest,
        "closest_diagonal": closest_ratio,
        "closest_L2": closest_l2,
        "closest_mu": closest_mu,
        "any_valid_peak": any_valid,
        "nan_found": nan_found,
        "verdict": verdict,
    }

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

    print("=" * 72)
    print("A313 PROBE -- OI-287-1: does the FULL B^4 arena (angular l + Hopf")
    print("              fiber) close mu -> 1.5 and restore the diagonal peak?")
    print("=" * 72)
    print("target (normalized):", TARGET)
    print("C-flat control diagonal:", results["C_flat"]["diagonal_ratio"],
          " (W=1, recovers A299)")
    print("geometric assignment: family n -> angular l=n, fiber k=n "
          "(P18 S2.2/S4.3)")
    print("-" * 72)
    for name, v in results.items():
        tag = "FULL-ARENA" if v["is_full_arena"] else (
            "pure-radial" if name in PURE_RADIAL else "control    ")
        print(f"{name:9s} [{tag}] diag={v['diagonal_ratio']}  "
              f"valid={v['valid_peak']} L2={v['L2_to_target']} "
              f"mu={v['mu']} muMoved={v['mu_moved_toward_1.5']} "
              f"<=20%={v['within_20pct']}")
    print("-" * 72)
    print(f"candidates declared: {len(results)} (all geometrically forced)")
    print(f"pure-radial mu closest to 1.5: {pr_closest_to_15}")
    print(f"angular-sector mu motion (C_ang - C_r3): {angular_motion:+.4f}")
    print(f"fiber-sector mu motion   (C_fiber - C_ang): {fiber_motion:+.4f}")
    print(f"sector supplying the motion: {sector_supplying_motion}")
    print(f"full-arena best mu: {fa_best_mu} ({fa_best})")
    print(f"full-arena mu > 1.40: {fa_mu_gt_140}")
    print(f"full-arena restores valid peak: {full_arena_valid}")
    print(f"closest non-control: {closest} diag={closest_ratio} "
          f"L2={closest_l2} mu={closest_mu}")
    print(f"any valid diagonal peak: {any_valid}")
    print(f"NaN found: {nan_found}")
    print(f"\nVERDICT: {verdict}")
    if verdict == "HIT":
        print("  -> a forced full-arena candidate restores BOTH a valid peak")
        print("     and the diagonal within 20%. The magnitude piece of")
        print("     OI-287-1 CLOSES: the residual was the dropped sectors.")
    elif verdict == "STRENGTHENED_PARTIAL":
        print("  -> adding the angular/fiber sector moves mu materially closer")
        print(f"     to 1.5 (or restores a peak); the {sector_supplying_motion}")
        print("     sector supplies the motion. Short of full closure.")
    else:
        print("  -> the full-arena measure does NOT improve on the pure-radial")
        print("     A311 result. The residual to 1.5 is NOT in the dropped")
        print("     angular/fiber sectors; the gap is elsewhere. This")
        print("     REDIRECTS the OI-287-1 magnitude search.")
    print(f"wrote {OUT}")


if __name__ == "__main__":
    main()
