from __future__ import annotations

import datetime
import json
import math
from fractions import Fraction as F
from pathlib import Path
from typing import Any

import numpy as np
import mpmath as mp
from scipy.integrate import solve_ivp
from scipy.optimize import brentq, minimize_scalar

ROOT = Path(__file__).resolve().parents[1]
OUT = ROOT / "results" / "utfanswf_embedding" / "result.json"
SPEC = ROOT / "research" / "utfanswf_embedding" / "UTFANSWF_EMBEDDING_SPEC.md"
SANSWF_RESULT = ROOT / "results" / "sanswf_symmetry" / "result.json"
ACTIVE_REGISTRY = ROOT / "datasets" / "v23_active_registry.json"
LEGACY_UNIFICATION = ROOT / "datasets" / "unification_tau.yaml"

ABS_TOL = 2.0e-10
MZ = 91.1876
MSUSY = 3000.0
ALPHA_MZ_PINNED = np.array([0.01695, 0.03377, 0.1184], dtype=float)
ALPHA3_MZ = 0.1184
GI2 = 0.5
GI = math.sqrt(GI2)
ALPHAU_INV_PRIMITIVE = 4.0 * math.pi / GI2  # 8*pi
R_SIGMA = 1.0 / (5.0 * GI)


def utc_now_z() -> str:
    return datetime.datetime.now(datetime.UTC).isoformat(timespec="seconds").replace("+00:00", "Z")


def suN_generators(n: int) -> np.ndarray:
    gs: list[np.ndarray] = []
    for i in range(n):
        for j in range(i + 1, n):
            s = np.zeros((n, n), dtype=complex)
            s[i, j] = s[j, i] = 0.5
            gs.append(s)
            a = np.zeros((n, n), dtype=complex)
            a[i, j] = -0.5j
            a[j, i] = 0.5j
            gs.append(a)
    for k in range(1, n):
        d = np.zeros(n, dtype=float)
        d[:k] = 1.0
        d[k] = -float(k)
        gs.append((np.diag(d) / math.sqrt(2.0 * k * (k + 1.0))).astype(complex))
    return np.asarray(gs)


def wedge2_representation(gens: np.ndarray) -> np.ndarray:
    n = gens.shape[1]
    cols: list[np.ndarray] = []
    for i in range(n):
        for j in range(i + 1, n):
            v = np.zeros(n * n, dtype=complex)
            v[i * n + j] = 1.0 / math.sqrt(2.0)
            v[j * n + i] = -1.0 / math.sqrt(2.0)
            cols.append(v)
    p = np.stack(cols, axis=1)
    eye = np.eye(n, dtype=complex)
    return np.asarray([p.conj().T @ (np.kron(t, eye) + np.kron(eye, t)) @ p for t in gens])


def anomaly_tensor(gens: np.ndarray) -> np.ndarray:
    prod = np.einsum("bij,cjk->bcik", gens, gens)
    anti = prod + np.swapaxes(prod, 0, 1)
    return np.einsum("aij,bcji->abc", gens, anti).real


def su5_anomaly_probe() -> dict[str, Any]:
    g5 = suN_generators(5)
    g10 = wedge2_representation(g5)
    g5bar = -g5.conjugate()
    d5 = anomaly_tensor(g5)
    d10 = anomaly_tensor(g10)
    d5bar = anomaly_tensor(g5bar)
    denom = float(np.sum(d5 * d5))
    a10 = float(np.sum(d10 * d5) / denom)
    a5bar = float(np.sum(d5bar * d5) / denom)
    res10 = float(np.max(np.abs(d10 - a10 * d5)))
    res5bar = float(np.max(np.abs(d5bar - a5bar * d5)))
    cancel = float(np.max(np.abs(d10 + d5bar)))
    return {
        "A_10_over_A_5": a10,
        "A_5bar_over_A_5": a5bar,
        "tensor_residual_10": res10,
        "tensor_residual_5bar": res5bar,
        "10_plus_5bar_anomaly_tensor_residual": cancel,
        "three_generations_residual": 3.0 * cancel,
        "pass": bool(abs(a10 - 1.0) <= ABS_TOL and abs(a5bar + 1.0) <= ABS_TOL and cancel <= ABS_TOL),
    }


def sm_anomaly_probe() -> dict[str, Any]:
    # One left-handed Weyl generation: multiplicity, Y, SU2 dim/T2, SU3 dim/T3, SU3 cubic sign.
    fields = [
        ("Q", 6, F(1, 6), 2, F(1, 2), 3, F(1, 2), 1),
        ("Uc", 3, F(-2, 3), 1, F(0), 3, F(1, 2), -1),
        ("Dc", 3, F(1, 3), 1, F(0), 3, F(1, 2), -1),
        ("L", 2, F(-1, 2), 2, F(1, 2), 1, F(0), 0),
        ("Ec", 1, F(1), 1, F(0), 1, F(0), 0),
    ]
    grav_y = sum(F(mult) * y for _, mult, y, *_ in fields)
    y3 = sum(F(mult) * y ** 3 for _, mult, y, *_ in fields)
    su2_y = sum(F(cdim) * t2 * y for _, _, y, _, t2, cdim, _, _ in fields)
    su3_y = sum(F(wdim) * t3 * y for _, _, y, wdim, _, _, t3, _ in fields)
    su3_cubic = sum(F(wdim) * F(sign) for _, _, _, wdim, _, cdim, _, sign in fields if cdim == 3)
    n_su2_doublets = 3 + 1
    return {
        "grav2_Y": str(grav_y),
        "Y3": str(y3),
        "SU2_2_Y": str(su2_y),
        "SU3_2_Y": str(su3_y),
        "SU3_3_relative": str(su3_cubic),
        "SU2_Witten_doublets": n_su2_doublets,
        "pass": bool(grav_y == 0 and y3 == 0 and su2_y == 0 and su3_y == 0 and su3_cubic == 0 and n_su2_doublets % 2 == 0),
    }


def residual_z6_probe() -> dict[str, Any]:
    # Same exact parent/residual structure carried by Row 28, recomputed here so the active
    # unified embedding cannot accidentally keep conventional DeltaB=1 decay active.
    lx = [-4, 2, 2]
    ex = [5, -1, -1]
    nx = [3, -3, -3]
    fields = []
    for g in range(3):
        fields += [
            ("Q", 6, F(1, 6), 0, 2, F(1, 2), 3, F(1, 2)),
            ("Uc", 3, F(-2, 3), -1, 1, F(1, 2), 3, F(0)),
            ("Dc", 3, F(1, 3), 1, 1, F(1, 2), 3, F(0)),
            ("L", 2, F(-1, 2), lx[g], 2, F(0), 1, F(1, 2)),
            ("Ec", 1, F(1), ex[g], 1, F(0), 1, F(0)),
            ("Nc", 1, F(0), nx[g], 1, F(0), 1, F(0)),
        ]
    grav = sum(F(mult) * F(x) for _, mult, y, x, *_ in fields)
    cubic = sum(F(mult) * F(x) ** 3 for _, mult, y, x, *_ in fields)
    xy2 = sum(F(mult) * F(x) * y ** 2 for _, mult, y, x, *_ in fields)
    x2y = sum(F(mult) * F(x) ** 2 * y for _, mult, y, x, *_ in fields)
    su3 = F(0)
    su2 = F(0)
    for _, _, _, x, dweak, t3, dcolor, t2 in fields:
        if t3:
            su3 += F(dweak) * t3 * F(x)
        if t2:
            su2 += F(dcolor) * t2 * F(x)

    z6 = {"Q": 0, "Uc": 5, "Dc": 1, "L": 2, "Ec": 5, "Nc": 3, "Hu": 1, "Hd": 5}
    def q6(*names: str) -> int:
        return sum(z6[n] for n in names) % 6
    required = {
        "Q_Hu_Uc": q6("Q", "Hu", "Uc"),
        "Q_Hd_Dc": q6("Q", "Hd", "Dc"),
        "L_Hd_Ec": q6("L", "Hd", "Ec"),
        "L_Hu_Nc": q6("L", "Hu", "Nc"),
        "Nc_Nc": q6("Nc", "Nc"),
        "Hu_Hd": q6("Hu", "Hd"),
        "Weinberg_(LHu)^2": (2 * q6("L", "Hu")) % 6,
    }
    dangerous = {
        "Uc_Dc_Dc": q6("Uc", "Dc", "Dc"),
        "Q_L_Dc": q6("Q", "L", "Dc"),
        "L_L_Ec": q6("L", "L", "Ec"),
        "L_Hu": q6("L", "Hu"),
        "Q_Q_Q_L": q6("Q", "Q", "Q", "L"),
        "Uc_Uc_Dc_Ec": q6("Uc", "Uc", "Dc", "Ec"),
    }
    anomalies = [grav, cubic, xy2, x2y, su2, su3]
    return {
        "parent_anomalies": {
            "grav2_U1X": str(grav), "U1X3": str(cubic), "U1X_U1Y2": str(xy2),
            "U1X2_U1Y": str(x2y), "SU2_2_U1X": str(su2), "SU3_2_U1X": str(su3),
        },
        "parent_all_zero": all(v == 0 for v in anomalies),
        "z6_charges": z6,
        "required_operators": required,
        "dangerous_operators": dangerous,
        "required_all_allowed": all(v == 0 for v in required.values()),
        "dangerous_all_forbidden": all(v != 0 for v in dangerous.values()),
        "single_proton_DeltaB1_active": False,
        "DeltaB2_active": False,
        "selection_rule": "DeltaB = 3*k",
        "isolated_proton_width_under_exact_Z6": 0.0,
        "pass": bool(all(v == 0 for v in anomalies) and all(v == 0 for v in required.values()) and all(v != 0 for v in dangerous.values())),
    }


def rotation12(theta: float) -> np.ndarray:
    c, s = math.cos(theta), math.sin(theta)
    return np.array([[c, s, 0.0], [-s, c, 0.0], [0.0, 0.0, 1.0]])


def rotation13(theta: float) -> np.ndarray:
    c, s = math.cos(theta), math.sin(theta)
    return np.array([[c, 0.0, s], [0.0, 1.0, 0.0], [-s, 0.0, c]])


def rotation23(theta: float) -> np.ndarray:
    c, s = math.cos(theta), math.sin(theta)
    return np.array([[1.0, 0.0, 0.0], [0.0, c, s], [0.0, -s, c]])


def seesaw_probe() -> dict[str, Any]:
    """Close the type-I seesaw mathematics for a declared benchmark slice.

    This does *not* derive the flavor texture or the heavy-neutrino spectrum.
    Those remain declared benchmark inputs. What is derived here is the full
    6x6 Majorana spectrum and its agreement with the Schur-complement seesaw
    expansion, so the active gate no longer relies only on the approximate
    low-energy formula.
    """
    vu = 174.0
    mr = np.array([5.0e13, 2.0e14, 1.2e15], dtype=float)
    m_ev = np.array([0.001, 0.0086, 0.050], dtype=float)
    m_gev = m_ev * 1.0e-9
    u = rotation23(0.79) @ rotation13(0.15) @ rotation12(0.58)

    # Deterministic benchmark parameterization in the diagonal-M_R basis,
    # equivalent to a real Casas-Ibarra slice with R=I.
    y = np.diag(np.sqrt(mr)) @ np.diag(np.sqrt(m_gev)) @ u.T / vu
    m_d = vu * y.T
    mr_inv = np.diag(1.0 / mr)
    schur = -(m_d @ mr_inv @ m_d.T)
    target = -(u @ np.diag(m_gev) @ u.T)
    schur_residual = float(np.max(np.abs(schur - target)))
    y_singular = np.linalg.svd(y, compute_uv=False)
    y_max = float(np.max(y_singular))

    # Explicit high-precision diagonalization of the complete 6x6 Majorana
    # matrix. High precision is required because the eigenvalues span ~27
    # orders of magnitude; double precision would erase the light spectrum.
    mp.mp.dps = 80
    vu_mp = mp.mpf("174")
    mr_mp = [mp.mpf("5e13"), mp.mpf("2e14"), mp.mpf("1.2e15")]
    mgev_mp = [mp.mpf("0.001e-9"), mp.mpf("0.0086e-9"), mp.mpf("0.050e-9")]

    def mp_r12(theta: str) -> mp.matrix:
        t = mp.mpf(theta); c, ss = mp.cos(t), mp.sin(t)
        return mp.matrix([[c, ss, 0], [-ss, c, 0], [0, 0, 1]])

    def mp_r13(theta: str) -> mp.matrix:
        t = mp.mpf(theta); c, ss = mp.cos(t), mp.sin(t)
        return mp.matrix([[c, 0, ss], [0, 1, 0], [-ss, 0, c]])

    def mp_r23(theta: str) -> mp.matrix:
        t = mp.mpf(theta); c, ss = mp.cos(t), mp.sin(t)
        return mp.matrix([[1, 0, 0], [0, c, ss], [0, -ss, c]])

    u_mp = mp_r23("0.79") * mp_r13("0.15") * mp_r12("0.58")
    y_mp = mp.diag([mp.sqrt(x) for x in mr_mp]) * mp.diag([mp.sqrt(x) for x in mgev_mp]) * u_mp.T / vu_mp
    md_mp = vu_mp * y_mp.T
    full = mp.matrix(6, 6)
    for i in range(3):
        for j in range(3):
            full[i, j] = mp.mpf("0")
            full[i, j + 3] = md_mp[i, j]
            full[i + 3, j] = md_mp.T[i, j]
            full[i + 3, j + 3] = mr_mp[i] if i == j else mp.mpf("0")

    evals, _ = mp.eigsy(full)
    exact_light_ev = sorted([abs(evals[i]) * mp.mpf("1e9") for i in range(3)])
    target_light_ev = sorted([mp.mpf("0.001"), mp.mpf("0.0086"), mp.mpf("0.050")])
    exact_heavy_gev = [evals[i] for i in range(3, 6)]
    light_rel = [abs((a - b) / b) for a, b in zip(exact_light_ev, target_light_ev)]
    heavy_rel = [abs((a - b) / b) for a, b in zip(exact_heavy_gev, mr_mp)]

    theta = m_d @ mr_inv
    theta_singular = np.linalg.svd(theta, compute_uv=False)
    theta_max = float(np.max(theta_singular))

    one_e15 = math.sqrt(1.0e15 * 0.05e-9) / vu
    one_point_two_e15 = math.sqrt(1.2e15 * 0.05e-9) / vu

    exact_light_max_rel = float(max(light_rel))
    exact_heavy_max_rel = float(max(heavy_rel))
    return {
        "scope": "DECLARED_BENCHMARK_SLICE__FULL_6X6_MATHEMATICAL_CLOSURE",
        "derived_here": [
            "Dirac mass matrix from the declared R=I benchmark parameterization",
            "Schur-complement light-neutrino mass matrix",
            "full 6x6 Majorana eigen-spectrum at 80-digit precision",
            "heavy-light mixing singular values",
        ],
        "not_derived_here": [
            "origin of the three M_R eigenvalues",
            "origin of the light-neutrino masses/mixing angles",
            "general complex Casas-Ibarra R matrix / CP phases",
            "full flavor/Yukawa origin",
        ],
        "v_u_GeV": vu,
        "M_R_GeV": mr.tolist(),
        "m_nu_input_eV": m_ev.tolist(),
        "mixing_angles_rad": {"theta12": 0.58, "theta13": 0.15, "theta23": 0.79, "CP_phase": 0.0},
        "mixing_matrix_orthogonality_residual": float(np.max(np.abs(u.T @ u - np.eye(3)))),
        "schur_complement_reconstruction_Linf_GeV": schur_residual,
        "exact_full_6x6_light_masses_eV": [float(x) for x in exact_light_ev],
        "exact_full_6x6_heavy_masses_GeV": [float(x) for x in exact_heavy_gev],
        "exact_light_vs_declared_max_relative_error": exact_light_max_rel,
        "exact_heavy_vs_declared_max_relative_shift": exact_heavy_max_rel,
        "heavy_light_mixing_singular_values": theta_singular.tolist(),
        "heavy_light_mixing_max": theta_max,
        "Y_singular_values": y_singular.tolist(),
        "Y_max": y_max,
        "sqrt_4pi": math.sqrt(4.0 * math.pi),
        "single_scale_Y_at_MR_1e15_mnu_0p05": one_e15,
        "single_scale_Y_at_MR_1p2e15_mnu_0p05": one_point_two_e15,
        "pass": bool(
            schur_residual <= 1e-20
            and exact_light_max_rel <= 1e-20
            and exact_heavy_max_rel <= 1e-20
            and theta_max < 1e-10
            and y_max < math.sqrt(4.0 * math.pi)
        ),
    }


bSM = np.array([41 / 10, -19 / 6, -7.0], dtype=float)
BSM = np.array([[199 / 50, 27 / 10, 44 / 5], [9 / 10, 35 / 6, 12], [11 / 10, 9 / 2, -26]], dtype=float)
bM = np.array([33 / 5, 1.0, -3.0], dtype=float)
BM = np.array([[199 / 25, 27 / 5, 88 / 5], [9 / 5, 25, 24], [11 / 5, 9, 14]], dtype=float)
Cadj = np.array([0.0, 2.0, 3.0])
wSigma = np.array([0.0, 2.0, 3.0])


def beta(_t: float, a: np.ndarray, b: np.ndarray, B: np.ndarray) -> np.ndarray:
    return b / (2.0 * math.pi) * a * a + a * a * (B @ a) / (8.0 * math.pi ** 2)


def run_up_from_mz(a_mz: np.ndarray, mu: float) -> np.ndarray:
    lo = solve_ivp(lambda t, a: beta(t, a, bSM, BSM), [math.log(MZ), math.log(MSUSY)], a_mz, rtol=1e-11, atol=1e-13)
    a_ms = lo.y[:, -1]
    inv_dr = 1.0 / a_ms - Cadj / (12.0 * math.pi)
    a_dr = 1.0 / inv_dr
    hi = solve_ivp(lambda t, a: beta(t, a, bM, BM), [math.log(MSUSY), math.log(mu)], a_dr, rtol=1e-11, atol=1e-13)
    return hi.y[:, -1]


def protected_inv_at_boundary(logmu: float, ratio: float = R_SIGMA) -> np.ndarray:
    mu = math.exp(logmu)
    inv = 1.0 / run_up_from_mz(ALPHA_MZ_PINNED, mu)
    delta = wSigma * math.log(1.0 / ratio) / (2.0 * math.pi)
    return inv - delta


def unification_spread(logmu: float, ratio: float = R_SIGMA) -> float:
    vals = protected_inv_at_boundary(logmu, ratio)
    return float(np.max(vals) - np.min(vals))


def canonical_threshold_repair_probe() -> dict[str, Any]:
    """Reproduce the corrected minimal-SUSY-SU(5) heavy-threshold closure.

    This is a pinned-input closure diagnostic, not a blind prediction. The
    matching scale and rho=MSigma/MX=1 are the frozen August-10 repair slice.
    The threshold coefficients are derived from the heavy representation beta
    weights and checked before solving the spectrum.
    """
    mu = 9.8e15
    rho = 1.0

    # One-loop heavy beta weights in GUT-normalized (U1, SU2, SU3) order.
    # HC is the color-triplet Higgs pair, Sigma is Sigma_3 + Sigma_8, and X is
    # the massive X/Y vector supermultiplet contribution.
    b_hc = np.array([2.0 / 5.0, 0.0, 1.0])
    b_sigma = np.array([0.0, 2.0, 3.0])
    b_x = np.array([-10.0, -6.0, -4.0])
    c1_vec = np.array([-1.0, 3.0, -2.0])
    c2_vec = np.array([5.0, -3.0, -2.0])
    c1_proj = {"HC": float(c1_vec @ b_hc), "Sigma": float(c1_vec @ b_sigma), "X": float(c1_vec @ b_x)}
    c2_proj = {"HC": float(c2_vec @ b_hc), "Sigma": float(c2_vec @ b_sigma), "X": float(c2_vec @ b_x)}

    inv = 1.0 / run_up_from_mz(ALPHA_MZ_PINNED, mu)
    c1 = float(c1_vec @ inv)
    c2 = float(c2_vec @ inv)

    # Corrected canonical solutions for rho = MSigma/MX:
    # L_HC = -(5*pi/6) C1
    # L_X  = -(pi/18) C2 + (1/3) ln rho
    # L_S  = -(pi/18) C2 - (2/3) ln rho
    l_hc = -(5.0 * math.pi / 6.0) * c1
    l_x = -(math.pi / 18.0) * c2 + (1.0 / 3.0) * math.log(rho)
    l_sigma = -(math.pi / 18.0) * c2 - (2.0 / 3.0) * math.log(rho)
    mhc = mu / math.exp(l_hc)
    mx = mu / math.exp(l_x)
    msigma = mu / math.exp(l_sigma)

    delta = (b_hc * l_hc + b_sigma * l_sigma + b_x * l_x) / (2.0 * math.pi)
    inv_u = inv - delta
    spread = float(np.max(inv_u) - np.min(inv_u))
    alpha_u_inv = float(np.mean(inv_u))

    projection_pass = bool(
        abs(c1_proj["HC"] + 12.0 / 5.0) <= ABS_TOL
        and abs(c1_proj["Sigma"]) <= ABS_TOL
        and abs(c1_proj["X"]) <= ABS_TOL
        and abs(c2_proj["HC"]) <= ABS_TOL
        and abs(c2_proj["Sigma"] + 12.0) <= ABS_TOL
        and abs(c2_proj["X"] + 24.0) <= ABS_TOL
    )
    spectrum_reproduction_pass = bool(
        abs(mhc / 1.8700319286833188e16 - 1.0) < 1e-10
        and abs(mx / 1.0831613303122448e16 - 1.0) < 1e-10
        and abs(msigma / 1.0831613303122448e16 - 1.0) < 1e-10
        and abs(alpha_u_inv - 25.791003228) < 1e-8
    )

    return {
        "scope": "PINNED_INPUT_CORRECTED_HEAVY_THRESHOLD_CLOSURE",
        "inputs": {
            "alpha_MZ_pinned": ALPHA_MZ_PINNED.tolist(),
            "MZ_GeV": MZ,
            "MSUSY_GeV": MSUSY,
            "matching_scale_mu_GeV": mu,
            "rho_MSigma_over_MX": rho,
            "running": "two-loop gauge-only SM->MSSM with one-loop MSbar->DRbar conversion",
        },
        "heavy_beta_weights": {"HC": b_hc.tolist(), "Sigma3_plus_Sigma8": b_sigma.tolist(), "XY": b_x.tolist()},
        "C1_projected_weights": c1_proj,
        "C2_projected_weights": c2_proj,
        "corrected_C2": "-(1/(2*pi))*(24*L_X + 12*L_Sigma)",
        "corrected_rho_solution": {
            "L_X": "-(pi/18)*C2 + (1/3)*ln(rho)",
            "L_Sigma": "-(pi/18)*C2 - (2/3)*ln(rho)",
        },
        "inverse_couplings_at_mu_before_heavy_matching": inv.tolist(),
        "C1": c1,
        "C2": c2,
        "log_mu_over_MHC": l_hc,
        "log_mu_over_MX": l_x,
        "log_mu_over_MSigma": l_sigma,
        "MHC_GeV": mhc,
        "MX_GeV": mx,
        "MSigma_GeV": msigma,
        "threshold_delta_inverse": delta.tolist(),
        "inverse_unified_couplings": inv_u.tolist(),
        "alphaU_inverse": alpha_u_inv,
        "closure_spread_inverse": spread,
        "threshold_weight_projection_pass": projection_pass,
        "frozen_repair_spectrum_reproduced": spectrum_reproduction_pass,
        "top_yukawa_extension": {
            "status": "OPEN_REQUIRES_DECLARED_TANBETA_AND_SUSY_YUKAWA_MATCHING_BOUNDARY",
            "reason": "The final hardening rules prohibit fitting tan(beta); without a declared MSSM Yukawa matching boundary the top-Yukawa contribution is not uniquely specified.",
        },
        "pass": bool(projection_pass and spectrum_reproduction_pass and spread <= 1e-10),
    }

def protected_unification_probe() -> dict[str, Any]:
    res = minimize_scalar(lambda x: unification_spread(x, R_SIGMA), bounds=(math.log(1e15), math.log(1e18)), method="bounded", options={"xatol": 1e-12})
    mu = math.exp(res.x)
    msig = R_SIGMA * mu
    vals = protected_inv_at_boundary(res.x, R_SIGMA)
    spread = float(np.max(vals) - np.min(vals))
    rel = spread / float(np.mean(vals))

    # Hypercharge normalization from the rank-3/rank-2 split H=diag(2,2,2,-3,-3), Y=H/6.
    h = np.diag([2.0, 2.0, 2.0, -3.0, -3.0])
    y = h / 6.0
    tr_y2 = float(np.trace(y @ y))
    g1_over_gy_sq = 2.0 * tr_y2

    return {
        "inputs": {"alpha_MZ_pinned": ALPHA_MZ_PINNED.tolist(), "MZ_GeV": MZ, "MSUSY_GeV": MSUSY},
        "MU_GeV": mu,
        "MSigma_GeV": msig,
        "MSigma_over_MU": R_SIGMA,
        "inverse_couplings_after_matching": vals.tolist(),
        "spread_inverse": spread,
        "relative_spread": rel,
        "relative_spread_percent": 100.0 * rel,
        "Tr_Y_squared": tr_y2,
        "g1_squared_over_gY_squared": g1_over_gy_sq,
        "active_threshold_content": "Sigma_3 + Sigma_8 with weight (0,2,3)",
        "legacy_XY_HC_threshold_branch_active": False,
        "pass": bool(100.0 * rel < 0.1 and abs(tr_y2 - 5.0 / 6.0) <= ABS_TOL and abs(g1_over_gy_sq - 5.0 / 3.0) <= ABS_TOL),
    }


def blind_uv_diagnostic() -> dict[str, Any]:
    """Audit the primitive alpha_U^{-1}=8*pi route after the full threshold repair.

    The historical calculation below retained only the Sigma_3+Sigma_8 threshold.
    It is reproduced for provenance, but the corrected minimal-SUSY-SU(5) map also
    requires H_C and X/Y threshold locations.  With alpha_U fixed and only alpha_3
    allowed to determine a scale, those additional heavy ratios are not determined
    by the present UTFANSWF-native boundary data.  Therefore the old 1.1456% number
    is no longer treated as an active blind electromagnetic prediction.
    """
    def run_down_legacy_sigma_only(logmu: float) -> tuple[float, float, np.ndarray]:
        mu = math.exp(logmu)
        ms = R_SIGMA * mu
        delta = wSigma * math.log(mu / ms) / (2.0 * math.pi)
        inv_u = np.full(3, ALPHAU_INV_PRIMITIVE) + delta
        a_u = 1.0 / inv_u
        hi = solve_ivp(lambda t, a: beta(t, a, bM, BM), [math.log(mu), math.log(MSUSY)], a_u, rtol=1e-11, atol=1e-13)
        a_dr = hi.y[:, -1]
        inv_ms = 1.0 / a_dr + Cadj / (12.0 * math.pi)
        a_ms = 1.0 / inv_ms
        lo = solve_ivp(lambda t, a: beta(t, a, bSM, BSM), [math.log(MSUSY), math.log(MZ)], a_ms, rtol=1e-11, atol=1e-13)
        return mu, ms, lo.y[:, -1]

    def qcdres(logmu: float) -> float:
        return float(run_down_legacy_sigma_only(logmu)[2][2] - ALPHA3_MZ)

    root = brentq(qcdres, math.log(1e14), math.log(1e18), xtol=1e-13)
    mu, ms, a = run_down_legacy_sigma_only(root)
    inv = 1.0 / a
    em_inv = (5.0 / 3.0) * inv[0] + inv[1]
    em_hold = (5.0 / 3.0) / ALPHA_MZ_PINNED[0] + 1.0 / ALPHA_MZ_PINNED[1]
    rel = (em_inv - em_hold) / em_hold

    # Full corrected threshold structure needs L_HC plus the X/Sigma combination.
    # The native split supplies R_SIGMA, but no independent H_C/MU or X/MU ratio.
    # The pinned closure spectrum cannot be imported into a blind EM prediction,
    # because it was solved using alpha_1(MZ) and alpha_2(MZ).
    missing_blind_boundaries = [
        "M_HC/M_U (or equivalent L_HC)",
        "M_X/M_U (or equivalent L_X)",
    ]
    available_native_boundary = {
        "alphaU_inverse_primitive": ALPHAU_INV_PRIMITIVE,
        "MSigma_over_MU_candidate": R_SIGMA,
        "alpha3_MZ_allowed_to_fix_one_scale": ALPHA3_MZ,
    }

    return {
        "status": "BLOCKED_BY_HEAVY_THRESHOLD_BOUNDARY",
        "active_blind_prediction_available": False,
        "reason": (
            "The corrected HC+XY+Sigma threshold map contains more independent heavy-threshold information "
            "than the current UTFANSWF-native blind boundary fixes. Importing the pinned closure masses would "
            "reintroduce alpha1/alpha2 information and violate the blind electromagnetic contract."
        ),
        "available_UTFANSWF_native_boundary": available_native_boundary,
        "missing_independent_heavy_boundaries": missing_blind_boundaries,
        "pinned_corrected_threshold_closure_allowed_as_blind_input": False,
        "legacy_sigma_only_diagnostic": {
            "status": "HISTORICAL_INCOMPLETE_THRESHOLD_DIAGNOSTIC_ONLY",
            "alphaU_inverse_primitive": ALPHAU_INV_PRIMITIVE,
            "uses_only_alpha3_to_fix_scale": True,
            "MU_GeV": mu,
            "MSigma_GeV": ms,
            "alpha_MZ_prediction": a.tolist(),
            "alpha_em_inverse_MZ_prediction": em_inv,
            "withheld_v22_alpha_em_inverse_MZ": em_hold,
            "absolute_percent_error": 100.0 * abs(rel),
            "old_1pct_gate_pass": bool(abs(rel) < 0.01),
            "incomplete_threshold_content": "Sigma_3 + Sigma_8 only; H_C and X/Y omitted",
        },
    }


def main() -> None:
    sanswf_status = "MISSING"
    sanswf_pass = False
    if SANSWF_RESULT.exists():
        obj = json.loads(SANSWF_RESULT.read_text(encoding="utf-8"))
        sanswf_status = str(obj.get("status", "UNKNOWN"))
        sanswf_pass = bool(obj.get("verdict") == "PASS" and sanswf_status == "PASS_SANSWF_SYMMETRY_STABILITY_TO_UTFANSWF_HANDOFF")

    su5_anom = su5_anomaly_probe()
    sm_anom = sm_anomaly_probe()
    z6 = residual_z6_probe()
    seesaw = seesaw_probe()
    legacy_protected_unification = protected_unification_probe()
    unification = canonical_threshold_repair_probe()
    blind_uv = blind_uv_diagnostic()

    legacy_exists = LEGACY_UNIFICATION.exists()
    registry_ok = False
    registry_obj: dict[str, Any] = {}
    if ACTIVE_REGISTRY.exists():
        registry_obj = json.loads(ACTIVE_REGISTRY.read_text(encoding="utf-8"))
        legacy_state = registry_obj.get("historical_not_active", {}).get("datasets/unification_tau.yaml", {}).get("status")
        registry_ok = bool(legacy_state == "HISTORICAL_LEGACY_NOT_ACTIVE")

    checks = {
        "sanswf_handoff_passes": sanswf_pass,
        "su5_10_plus_5bar_cubic_anomaly_cancels": bool(su5_anom["pass"]),
        "sm_decomposition_anomalies_cancel_per_generation": bool(sm_anom["pass"]),
        "seesaw_full_6x6_benchmark_closure_and_perturbativity_passes": bool(seesaw["pass"]),
        "corrected_full_heavy_threshold_closure_passes": bool(unification["pass"]),
        "residual_exact_Z6_proton_stability_structure_passes": bool(z6["pass"]),
        "legacy_unification_tau_preserved_but_marked_inactive": bool(legacy_exists and registry_ok),
        "legacy_single_proton_finite_lifetime_not_active": bool(not z6["single_proton_DeltaB1_active"] and z6["isolated_proton_width_under_exact_Z6"] == 0.0),
    }

    claim_guards = {
        "corrected_threshold_closure_uses_pinned_low_energy_couplings_and_is_not_blind_prediction": True,
        "top_yukawa_extension_not_faked_without_tanbeta_boundary": bool(unification["top_yukawa_extension"]["status"].startswith("OPEN_")),
        "primitive_alphaU_boundary_diagnostic_does_not_block_embedding_existence": True,
        "blind_full_threshold_boundary_not_faked": bool(
            blind_uv["status"] == "BLOCKED_BY_HEAVY_THRESHOLD_BOUNDARY"
            and blind_uv["active_blind_prediction_available"] is False
            and blind_uv["pinned_corrected_threshold_closure_allowed_as_blind_input"] is False
        ),
        "legacy_sigma_only_near_fail_preserved_as_history": bool(
            blind_uv["legacy_sigma_only_diagnostic"]["status"] == "HISTORICAL_INCOMPLETE_THRESHOLD_DIAGNOSTIC_ONLY"
        ),
        "seesaw_full_diagonalization_is_still_a_declared_RI_benchmark_not_general_flavor_solution": True,
        "first_generation_yukawa_origin_not_claimed": True,
        "conventional_DeltaB1_proton_decay_branch_retained_only_as_history": True,
        "SO10_and_heterotic_alternatives_not_claimed_closed_by_SU5_gate": True,
        "historical_rows_not_renumbered": True,
    }

    verdict = "PASS" if all(checks.values()) and all(claim_guards.values()) else "FAIL"
    status = "PASS_UTFANSWF_UNIFIED_EMBEDDING__BLIND_UV_BOUNDARY_OPEN" if verdict == "PASS" else "FAIL_UTFANSWF_UNIFIED_EMBEDDING"

    result = {
        "row": 0.4,
        "stage": "UTFANSWF",
        "gate": "utfanswf_unified_embedding",
        "signature": "UTFANSWF_UNIFIED_EMBEDDING_PREFLIGHT_v2_FULL_THRESHOLD_BLIND_BOUNDARY_AUDIT",
        "context": "LIVE",
        "tier": "Critical",
        "blocking": True,
        "timestamp_utc": utc_now_z(),
        "verdict": verdict,
        "status": status,
        "upstream_dependency": {
            "required": "PASS_SANSWF_SYMMETRY_STABILITY_TO_UTFANSWF_HANDOFF",
            "observed": sanswf_status,
            "pass": sanswf_pass,
        },
        "checks": checks,
        "claim_guards": claim_guards,
        "metrics": {
            "SU5_representation_anomaly": su5_anom,
            "SM_decomposition_anomaly": sm_anom,
            "seesaw": seesaw,
            "corrected_full_heavy_threshold_closure": unification,
            "legacy_sigma_only_protected_unification": legacy_protected_unification,
            "blind_UV_boundary_audit": blind_uv,
            "exact_Z6_proton_stability": z6,
            "active_registry": registry_obj,
        },
        "active_prediction_propagation": {
            "single_proton_DeltaB1": "FORBIDDEN under exact Z6 candidate",
            "DeltaB2": "FORBIDDEN under exact Z6 candidate",
            "selection_rule": "DeltaB=3*k",
            "isolated_proton_width": 0.0,
            "old_finite_proton_lifetime_forecast": "HISTORICAL_SUPERSEDED_BRANCH",
            "active_unification_branch": "corrected minimal-SUSY-SU5 HC+XY+Sigma heavy-threshold closure with pinned low-energy couplings",
            "legacy_unification_tau_yaml": "HISTORICAL_LEGACY_NOT_ACTIVE",
        },
        "corrections_to_legacy_writeup": [
            "The v23 active proton branch cannot simultaneously advertise a conventional finite DeltaB=1 single-proton lifetime; that forecast is retained only as a historical branch after the exact-Z6 candidate is activated.",
            "The v23 active unification closure uses the corrected full HC+XY+Sigma threshold map; the prior Sigma-only shortcut is retained only as a legacy diagnostic.",
            "For M_R=1.2e15 GeV, m_nu=0.05 eV, v_u=174 GeV, the single-scale Yukawa is %.12g; the older ~1.44 shorthand is not the exact value for those declared inputs." % seesaw["single_scale_Y_at_MR_1p2e15_mnu_0p05"],
        ],
        "open_boundaries": [
            "The corrected HC+XY+Sigma heavy-threshold result is a closure/consistency diagnostic using pinned low-energy couplings, not a blind derivation of the electroweak couplings.",
            "The top-Yukawa extension remains open because no non-fitted tan(beta)/MSSM Yukawa matching boundary is declared in the current bundle.",
            "The primitive alphaU^-1=8pi candidate is not a source-complete blind prediction after the corrected HC+XY+Sigma threshold repair: independent H_C/M_U and X/M_U information is still missing. The older ~1.1456% near-fail is retained only as a Sigma-only historical diagnostic.",
            "The full 6x6 seesaw diagonalization closes the mathematics only for the declared benchmark slice; the M_R spectrum, light-neutrino inputs, CP phases, and full flavor/Yukawa origin remain open.",
            "SO(10) and heterotic alternatives remain possible embedding branches but are not numerically closed by this SU(5)-focused preflight.",
            "The DeltaB=3 Wilson coefficient/scale remains unfixed, so no trinucleon lifetime is produced here.",
        ],
        "specification": str(SPEC.relative_to(ROOT)),
        "notes": "Critical unified-embedding preflight. Activates corrected full heavy-threshold closure + exact-Z6 proton-stability propagation while preserving incomplete/legacy branches as historical records; blind electromagnetic UV closure remains blocked by missing independent heavy-threshold boundary data.",
    }

    OUT.parent.mkdir(parents=True, exist_ok=True)
    OUT.write_text(json.dumps(result, indent=2, sort_keys=True), encoding="utf-8")
    print(json.dumps({"row": 0.4, "stage": "UTFANSWF", "verdict": verdict, "status": status, "MU_GeV": unification["inputs"]["matching_scale_mu_GeV"], "unification_spread_percent": 100.0 * unification["closure_spread_inverse"] / unification["alphaU_inverse"], "seesaw_Ymax": seesaw["Y_max"], "blind_UV_status": blind_uv["status"], "legacy_sigma_only_error_percent": blind_uv["legacy_sigma_only_diagnostic"]["absolute_percent_error"]}, indent=2))
    if verdict != "PASS":
        raise SystemExit(1)


if __name__ == "__main__":
    main()
