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

DERIVATION probe, not a fit. A307 localized the missing factor of the A299
selection rule to the INNER-PRODUCT MEASURE: with the density measure rho dx and
the rho-weighted generalized eigenproblem H psi = lambda rho psi, mu lifts from
the A299 flat-measure base 1.12 to 1.34. A309 then showed the exact measure is
NOT pinned by fit (rho dx is not the L2-closest weight; x^1 dx is, at L2=0.33;
rho dx is only a stationary point of the lift over its power family). So the open
question is no longer "which measure fits best" but "is a particular measure
FORCED by the geometry."

This probe tests exactly that. Paper 18 S2 places the master operator on
sections over the 4-ball B^4 with boundary S^3, Hopf-fibered. The base metric is
ds^2 = dr^2 + r^2 dOmega_3^2 (P18 S2.1), so the 4-volume element is r^3 dr dOmega_3
and the GEOMETRICALLY FORCED radial measure for rotationally-symmetric overlaps
is W_geom(r) = r^3 (the surface of S^3 of radius r scales as r^3), NOT flat dr.
P18 S2.2 gives the radial operator explicitly with the -(3/r) d/dr term:

    D_{B^4}^2 = -d^2/dr^2 - (3/r) d/dr + L^2/r^2 + R/4   (P18 eq, S2.2)

The radial operator -d^2/dr^2 - (3/r) d/dr is the genuine 4D radial Laplacian; it
is self-adjoint in the r^3 dr measure, where it takes the Sturm-Liouville form
-1/r^3 d/dr(r^3 d/dr). The thesis under test: the correct inner product on the
radial reduction is this geometric arena volume r^3 dr, determined not free, and
it is what lifts mu. The A299/A307 canonical H is the flat-1D reduction; the
corpus-faithful object is the 4D-radial operator with its r^3 measure.

Note the analytic relation between the forced geometric measure and A307's
density measure: rho(r) = 16 pi^3 r^3 + 3 pi^2 r^2 + 2 pi r, so rho(r) ~ 16 pi^3
r^3 to leading order near the dominant cubic term, and rho(r)/r^3 -> 16 pi^3 as
r grows on [0,1]. The density measure and the 4D radial measure r^3 SHARE the
leading r^3 behaviour: rho dx is, to leading order, the arena volume r^3 dr times
a constant. A307's lift and the geometric measure are the same object at leading
order.

PRE-REGISTERED CANDIDATES (declared before running; all geometrically motivated,
no free fitting; FORCED candidates distinguished from comparison candidates):

  FORCED by the B^4 arena (the thesis):
  (G-r3)   W = r^3 on the canonical A299 H. The 4D radial volume element of B^4
           applied as the Sturm-Liouville weight to the same operator A299/A307
           use, solved as the generalized H psi = lambda r^3 psi and the overlap
           formed in the r^3 measure. This is the FORCED measure if the thesis
           holds.
  (G-r3w)  The genuine 4D-radial operator H_4D = -d^2/dr^2 - (3/r) d/dr + V(r) of
           P18 S2.2, self-adjoint in the r^3 dr measure (Sturm-Liouville form
           -1/r^3 d/dr(r^3 d/dr) + V). Built from the weak form so the discrete
           generalized eigenproblem A psi = lambda W psi is symmetric-definite
           with W = r^3 dr. This is the corpus-faithful arena object: the right
           operator AND the right measure together.

  COMPARISON (declared, geometrically motivated, not the forced arena volume):
  (G-flat) W = 1, canonical A299 H. CONTROL, the un-geometric choice; must
           recover (1, 1.12, 7.50).
  (G-rho)  W = rho(r), canonical A299 H. A307's density measure; expected
           mu ~ 1.34. The leading-r^3 cousin of G-r3.
  (G-Hopf) W = r^3 * sqrt(1 - r^2), canonical A299 H. A declared Hopf-motivated
           weight: the r^3 arena volume times the Hopf-fiber lapse m = sqrt(1-u^2)
           read radially (the S^1 fiber contributes a lapse factor on the
           boundary). Strictly positive on (0,1). Tests whether the fiber measure
           sharpens or degrades the arena lift.

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

Operator (canonical A299/A307, the flat-1D reduction):
  H = -d^2/dx^2 + [rho'(x)]^2/(2 Omega^2) + E_self x^2 (1-x)^2
  rho'(x) = 48 pi^3 x^2 + 6 pi^2 x + 2 pi
  Omega   = 4 pi^3 + pi^2 + pi
  E_self  = 13.177
  rho(x)  = 16 pi^3 x^3 + 3 pi^2 x^2 + 2 pi x
  Dirichlet [0,1], N=1500.  C = (2, 3, 16), target normalized (1, 1.5, 8).
  Box control psi_box[n] = sqrt(2) sin((n+1) pi x).

PRE-REGISTERED VERDICT (declared before running):
  HIT iff the geometrically-forced measure (r^3, or the 4D-radial operator in its
       r^3 measure) yields BOTH valid_peak AND all three normalized entries within
       20% of (1, 1.5, 8).
  STRENGTHENED-PARTIAL iff r^3 (or the 4D-radial G-r3w) lifts mu materially toward
       1.5 (mu in [1.3, 1.7]) AND is geometrically forced (not a fit parameter).
       This upgrades A307 from "factor lives in some measure" to "the measure is
       the forced 4D arena volume," even short of full peak.
  NULL iff the geometric measure does not lift mu (the measure-locus of A307 is
       real but its origin is not the arena volume).

Multiple-comparison caveat: 5 candidates declared (2 FORCED, 3 COMPARISON),
counted. Honesty over a hit. W > 0 on (0,1) for every generalized problem.
"""
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
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_geomeasure.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 is the canonical flat-1D operator; W carries the measure. psi
    W-normalized: sum(psi^2 * weight) dx = 1.
    """
    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")
    Wm = np.diag(Wdiag)
    w, v = gen_eigh(H, Wm)
    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 eigfuns_radial4d(potential):
    """The corpus-faithful 4D-radial operator (P18 S2.2):

        H_4D = -d^2/dr^2 - (3/r) d/dr + V(r),

    self-adjoint in the r^3 dr measure, Sturm-Liouville form
    -1/r^3 d/dr(r^3 d/dr) + V. Built from the weak (bilinear) form on the same
    grid so the discrete generalized problem A psi = lambda W psi is symmetric
    positive-definite with W the r^3 dr mass matrix. Dirichlet at r=1; the
    natural r^3 vanishing handles the r=0 endpoint.

      A_ij = integral r^3 [ R_i' R_j' + V R_i R_j ] dr   (stiffness + potential)
      W_ij = integral r^3 R_i R_j dr                     (mass / measure)

    On the uniform interior grid we use the standard finite-difference stencil
    with the r^3 weight folded into both A and W, giving a symmetric generalized
    eigenproblem whose eigenfunctions are W-orthonormal in r^3 dr. This is the
    genuine 4D radial Laplacian problem of the B^4 arena, not the flat-1D
    reduction.
    """
    x, dx = grid()                       # interior nodes r in (0,1)
    w3 = x**3                             # r^3 measure on the nodes
    V = potential(x)
    n = N - 1

    # Stiffness for -1/r^3 d/dr(r^3 d/dr): symmetric form sum r^3 R' S' dr.
    # Edge weights r^3 evaluated at the midpoints between consecutive nodes
    # keep the operator self-adjoint in r^3 dr. Midpoint r_{i+1/2}^3.
    r_mid = x + 0.5 * dx                  # midpoint to the right of node i
    w3_mid = r_mid**3                     # r^3 at midpoints (length n; last is
    # the midpoint between node n-1 and the Dirichlet wall at r=1)
    # Left midpoint of node 0 is between r=0 wall and node 0:
    r_mid_left0 = x[0] - 0.5 * dx
    w3_mid_left0 = max(r_mid_left0, 0.0)**3

    A = np.zeros((n, n))
    # diagonal: (w3_mid[i] + w3_mid[i-1]) / dx^2  + r^3 V dx folded
    for i in range(n):
        wl = w3_mid[i - 1] if i > 0 else w3_mid_left0
        wr = w3_mid[i]
        A[i, i] = (wl + wr) / dx**2
        if i + 1 < n:
            A[i, i + 1] = -wr / dx**2
            A[i + 1, i] = -wr / dx**2
    # potential term: r^3 V on the diagonal (mass-lumped)
    A += np.diag(w3 * V)
    A = 0.5 * (A + A.T)

    # mass / measure matrix W = diag(r^3) (mass-lumped r^3 dr)
    Wm = np.diag(w3)

    # eigenvalues of A psi = lambda (1/dx^2 scaling already in A;
    # multiply A potential/mass by dx^2? Keep consistent: A has 1/dx^2 in
    # stiffness and r^3 V (no 1/dx^2) in potential. To match, scale W by dx^2
    # so that lambda has the right units: A_total psi = lambda W psi with
    # W including the dx^2 from the kinetic normalization.)
    # Cleaner: divide A's stiffness already carries 1/dx^2; the mass matrix in a
    # finite-difference (not FEM) reading is W = diag(r^3) and the eigenvalue is
    # lambda = <psi|H_4D|psi>_{r^3} / <psi|psi>_{r^3}. The stiffness sum
    # approximates integral r^3 R'^2 dr / dx (one dx absorbed); fold the
    # remaining dx by scaling W by dx^... -> instead solve the Rayleigh form
    # directly below. Here use W = diag(r^3) and stiffness already ~ /dx^2 *
    # (per-cell), consistent with hamiltonian() which also omits the dx volume.
    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 * w3) * dx))
        vv[:, k] = vv[:, k] / nrm
    return x, dx, w[:3], vv


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.

    measure is the weight array W(x) on the grid. W = ones reproduces the A299
    flat overlap exactly.
    """
    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(v, 4) for v in ratio],
        "tau": round(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(mu, 4),
        "mu_in_partial_band": bool(1.3 <= mu <= 1.7),
        "eigenvalues": [round(float(v), 4) for v in lam],
        "node_counts": [node_count(psi[:, n]) for n in range(3)],
    }


# FORCED-by-arena candidates: the 4D radial volume r^3 and the 4D-radial
# operator in its r^3 measure. The thesis stands or falls on these.
FORCED = {"G_r3", "G_r3w_radial4d"}


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
    rho_x = rho(x)
    hopf = x**3 * np.sqrt(1.0 - x**2)     # r^3 * lapse m = sqrt(1-r^2)

    results = {}

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

    # (G-r3) FORCED: W = r^3, canonical H, generalized H psi = lambda r^3 psi
    _, _, lam_r3, psi_r3 = eigfuns_weighted(Vpot, lambda xx: xx**3)
    results["G_r3"] = evaluate(psi_r3, lam_r3, x, dx, psi_box, lam_box, r3)

    # (G-r3w) FORCED: the corpus-faithful 4D-radial operator in its r^3 measure
    _, _, lam_r3w, psi_r3w = eigfuns_radial4d(Vpot)
    results["G_r3w_radial4d"] = evaluate(psi_r3w, lam_r3w, x, dx,
                                         psi_box, lam_box, r3)

    # (G-rho) COMPARISON: A307's density measure
    _, _, lam_rho, psi_rho = eigfuns_weighted(Vpot, rho)
    results["G_rho"] = evaluate(psi_rho, lam_rho, x, dx,
                                psi_box, lam_box, rho_x)

    # (G-Hopf) COMPARISON: r^3 * Hopf lapse sqrt(1-r^2)
    _, _, lam_h, psi_h = eigfuns_weighted(
        Vpot, lambda xx: xx**3 * np.sqrt(1.0 - xx**2))
    results["G_Hopf"] = evaluate(psi_h, lam_h, x, dx, psi_box, lam_box, hopf)

    # ---- analytic leading-term fact: rho(r) ~ 16 pi^3 r^3 ----------------
    # rho(r)/r^3 = 16 pi^3 + 3 pi^2 / r + 2 pi / r^2; the leading constant is
    # 16 pi^3. Measure the leading coefficient of rho as a polynomial in r:
    lead_coeff_rho = 16 * PI**3
    # ratio rho(r)/r^3 at the right endpoint r -> 1: rho(1)/1 = 16pi^3+3pi^2+2pi
    rho_over_r3_at_1 = rho(1.0) / 1.0
    # fraction of rho(1) carried by the r^3 term:
    r3_fraction_at_1 = (16 * PI**3) / rho(1.0)

    # ---- aggregate -------------------------------------------------------
    flat_mu = results["G_flat"]["mu"]
    non_control = {k: v for k, v in results.items() if k != "G_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_forced_arena_measure"] = k in FORCED

    any_valid = any(v["valid_peak"] for v in results.values())
    hit_forced = any(results[k]["valid_peak"] and results[k]["within_20pct"]
                     for k in FORCED)
    hit_any = any(v["valid_peak"] and v["within_20pct"]
                  for kk, v in non_control.items())

    # STRENGTHENED-PARTIAL: a FORCED candidate lifts mu into [1.3,1.7]
    forced_in_band = [k for k in FORCED
                      if results[k]["mu_in_partial_band"]]
    forced_lifts = len(forced_in_band) > 0

    # NULL: no forced candidate lifts mu toward 1.5 at all
    forced_moved = any(results[k]["mu_moved_toward_1.5"] for k in FORCED)

    if hit_forced:
        verdict = "HIT"
    elif forced_lifts:
        verdict = "STRENGTHENED_PARTIAL"
    elif not forced_moved:
        verdict = "NULL"
    else:
        # forced candidate moved mu toward 1.5 but not into the band: an honest
        # weak-partial, reported as STRENGTHENED_PARTIAL only if in band, else
        # NULL-leaning. Keep distinct label.
        verdict = "WEAK_PARTIAL"

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

    summary = {
        "addendum": "A311",
        "addresses": "OI-287-1",
        "builds_on": ["A307", "A309"],
        "thesis": ("is the inner-product measure forced by the B^4/S^3/Hopf "
                   "radial geometry (r^3 dr), not searched?"),
        "target_normalized": TARGET,
        "flat_control_diagonal": results["G_flat"]["diagonal_ratio"],
        "flat_control_mu": flat_mu,
        "constants": {"OMEGA_alpha_inv": OMEGA, "E_self": E_SELF},
        "candidates": results,
        "n_candidates_declared": len(results),
        "forced_arena_candidates": sorted(FORCED),
        "comparison_candidates": sorted(
            set(results) - FORCED - {"G_flat"}) + ["G_flat (control)"],
        "G_r3_diagonal": results["G_r3"]["diagonal_ratio"],
        "G_r3_mu": results["G_r3"]["mu"],
        "G_r3_L2": results["G_r3"]["L2_to_target"],
        "G_r3w_radial4d_diagonal": results["G_r3w_radial4d"]["diagonal_ratio"],
        "G_r3w_radial4d_mu": results["G_r3w_radial4d"]["mu"],
        "G_r3w_radial4d_L2": results["G_r3w_radial4d"]["L2_to_target"],
        "G_rho_diagonal": results["G_rho"]["diagonal_ratio"],
        "G_rho_mu": results["G_rho"]["mu"],
        "closest_noncontrol": closest,
        "closest_diagonal": closest_ratio,
        "closest_L2": closest_l2,
        "closest_mu": closest_mu,
        "any_valid_peak": any_valid,
        "forced_candidates_in_partial_band": forced_in_band,
        "forced_lifts_mu": bool(forced_lifts),
        "forced_moved_toward_1.5": bool(forced_moved),
        "analytic_rho_vs_r3": {
            "rho_leading_coefficient": round(lead_coeff_rho, 6),
            "note": ("rho(r) = 16 pi^3 r^3 + 3 pi^2 r^2 + 2 pi r; leading term "
                     "16 pi^3 r^3 shares the r^3 behaviour of the 4D arena "
                     "volume measure"),
            "rho_over_r3_at_r1": round(rho_over_r3_at_1, 6),
            "r3_term_fraction_of_rho_at_r1": round(r3_fraction_at_1, 6),
        },
        "x1_fit_reference_A309": 0.3252,
        "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("A311 PROBE -- OI-287-1: is the inner-product measure FORCED by the")
    print("              B^4/S^3/Hopf radial geometry (r^3 dr)?")
    print("=" * 72)
    print("target (normalized):", TARGET)
    print("G-flat control diagonal:", results["G_flat"]["diagonal_ratio"],
          " (W=1, recovers A299)")
    print("-" * 72)
    for name, v in results.items():
        tag = "FORCED" if v["is_forced_arena_measure"] else "compar"
        print(f"{name:18s} [{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)} "
          f"(forced arena: {sorted(FORCED)})")
    print(f"closest non-control: {closest} diag={closest_ratio} "
          f"L2={closest_l2} mu={closest_mu}")
    print(f"FORCED candidates lifting mu into [1.3,1.7]: {forced_in_band}")
    print(f"any valid diagonal peak: {any_valid}")
    print("analytic: rho(r) ~ 16 pi^3 r^3 leading; r^3 term carries "
          f"{r3_fraction_at_1:.3f} of rho(1)")
    print(f"NaN found: {nan_found}")
    print(f"\nVERDICT: {verdict}")
    if verdict == "HIT":
        print("  -> the forced 4D arena measure r^3 restores BOTH a valid peak")
        print("     and the diagonal within 20%. The measure is GEOMETRICALLY")
        print("     FORCED; OI-287-1 magnitude reduces to the arena volume.")
    elif verdict == "STRENGTHENED_PARTIAL":
        print("  -> the forced 4D arena measure r^3 lifts mu into [1.3,1.7].")
        print("     A307 upgrades from 'factor lives in some measure' to 'the")
        print("     measure is the forced 4D arena volume r^3 dr', short of")
        print(f"     full peak. Forced lifters: {forced_in_band}.")
    elif verdict == "WEAK_PARTIAL":
        print("  -> the forced arena measure moves mu toward 1.5 but not into")
        print("     the band [1.3,1.7]. Geometric origin plausible, not strong.")
    else:
        print("  -> the forced arena measure does NOT lift mu. The measure-locus")
        print("     of A307 is real but its origin is not the arena volume r^3.")
    print(f"wrote {OUT}")


if __name__ == "__main__":
    main()
