#!/usr/bin/env python3
"""Reproduce the authoritative 2026-08-19 UTFANSWF prediction additions.

This module is intentionally an auxiliary research reproducer, not a new numbered
canonical gate. It regenerates the promoted numerical/structural outputs from the
frozen v23 inputs and writes results/new_predictions/result.json.
"""
from __future__ import annotations

import json
import math
from pathlib import Path
from typing import Any

import numpy as np
import pandas as pd

ROOT = Path(__file__).resolve().parents[2]
HERE = Path(__file__).resolve().parent
SUPPORT = HERE / "support"
OUTDIR = ROOT / "results" / "new_predictions"
OUTDIR.mkdir(parents=True, exist_ok=True)

MODEL_NAMES = ["L", "Ql", "Qmix", "LcutM", "QlCutM"]
MUD_GEV = 0.00470
MS_GEV = 0.1270
LAMBDA_GEV = 0.3397960719428473
R0_LAMBDA = 0.820
R0_S = 3.277
MN_ISOQCD_GEV = 0.9412231331333769
SQRT_T0_GEV_INV = (R0_LAMBDA / R0_S) / LAMBDA_GEV


def load_json(path: Path) -> dict[str, Any]:
    return json.loads(path.read_text(encoding="utf-8"))


def design(name: str, l: np.ndarray, s: np.ndarray, c: np.ndarray) -> np.ndarray:
    one = np.ones_like(l)
    if name == "L":
        return np.column_stack([one, l, s, c])
    if name == "Ql":
        return np.column_stack([one, l, s, c, l * l])
    if name == "Qmix":
        return np.column_stack([one, l, s, c, l * l, s * s, l * s])
    if name == "LcutM":
        return np.column_stack([one, l, s, c, c * l, c * s])
    if name == "QlCutM":
        return np.column_stack([one, l, s, c, l * l, c * l, c * s])
    raise KeyError(name)


def fit_model(df: pd.DataFrame, name: str) -> tuple[np.ndarray, np.ndarray, float]:
    l = df.xl_rgi.to_numpy(float)
    s = df.xs_rgi.to_numpy(float)
    y = df.Y_N.to_numpy(float)
    sig = df.sigY_N.to_numpy(float)
    c = df.cut.to_numpy(float)
    X = design(name, l, s, c)
    W = 1.0 / (sig * sig)
    H = (X.T * W) @ X
    cov0 = np.linalg.pinv(H)
    coef = cov0 @ ((X.T * W) @ y)
    res = y - X @ coef
    chi2 = float(np.sum((res / sig) ** 2))
    dof = len(y) - X.shape[1]
    cov = cov0 * max(chi2 / dof, 1.0)
    k = X.shape[1]
    n = len(y)
    aicc = chi2 + 2 * k + 2 * k * (k + 1) / (n - k - 1)
    return coef, cov, aicc


def fit_bank(df: pd.DataFrame) -> tuple[dict[str, tuple[np.ndarray, np.ndarray, float]], np.ndarray]:
    fits = {m: fit_model(df, m) for m in MODEL_NAMES}
    aicc = np.array([fits[m][2] for m in MODEL_NAMES])
    weights = np.exp(-0.5 * (aicc - aicc.min()))
    weights /= weights.sum()
    return fits, weights


def derivative_ud(name: str, coef: np.ndarray, l: float, s: float) -> float:
    if name == "L":
        return float(coef[1])
    if name == "Ql":
        return float(coef[1] + 2 * coef[4] * l)
    if name == "Qmix":
        return float(coef[1] + 2 * coef[4] * l + coef[6] * s)
    if name == "LcutM":
        return float(coef[1])
    if name == "QlCutM":
        return float(coef[1] + 2 * coef[4] * l)
    raise KeyError(name)


def predict_scalar(name: str, coef: np.ndarray, l: float, s: float) -> float:
    X = design(name, np.array([l]), np.array([s]), np.array([0.0]))
    return float((X @ coef)[0])


def qcd_central(df: pd.DataFrame) -> tuple[float, float, float]:
    fits, w = fit_bank(df)
    l = SQRT_T0_GEV_INV * MUD_GEV
    s = SQRT_T0_GEV_INV * MS_GEV
    dud = sum(w[i] * derivative_ud(m, fits[m][0], l, s) for i, m in enumerate(MODEL_NAMES))
    sigma_mev = 1000.0 * MUD_GEV * dud
    kud = MUD_GEV * dud / MN_ISOQCD_GEV
    return float(dud), float(sigma_mev), float(kud)


def qcd_delta(df: pd.DataFrame, frac: float) -> float:
    fits, w = fit_bank(df)
    l0 = SQRT_T0_GEV_INV * MUD_GEV
    s0 = SQRT_T0_GEV_INV * MS_GEV
    vals = []
    for i, m in enumerate(MODEL_NAMES):
        co = fits[m][0]
        y0 = predict_scalar(m, co, l0, s0)
        y1 = predict_scalar(m, co, l0 * (1.0 + frac), s0)
        vals.append(1000.0 * (y1 - y0) / SQRT_T0_GEV_INV)
    return float(np.dot(w, np.array(vals)))


def jackknife_sigma(vals: list[float]) -> float:
    a = np.asarray(vals, dtype=float)
    n = len(a)
    mean = float(a.mean())
    return math.sqrt((n - 1) / n * float(np.sum((a - mean) ** 2)))


def reproduce_qcd() -> dict[str, Any]:
    csv_path = ROOT / "research" / "proton_origin" / "PROTON_ORIGIN_CLS_RGI_QUARK_INPUT_V2_2026-08-16.csv"
    df = pd.read_csv(csv_path)
    _, sigma_central, kud = qcd_central(df)
    betas = sorted(float(x) for x in df.beta.unique())

    sigma_leave = {f"{b:.2f}": qcd_central(df[df.beta != b])[1] for b in betas}
    sigma_beta_jk = jackknife_sigma(list(sigma_leave.values()))

    finite = load_json(SUPPORT / "UTFANSWF_V23_FINITE_QUARK_MASS_RESPONSE_PREDICTIONS_FROZEN_2026-08-19.json")
    promoted = load_json(SUPPORT / "UTFANSWF_V23_LIGHT_QUARK_NUCLEON_RESPONSE_CURVE_FROZEN_2026-08-19.json")

    frac_map = {
        "-20%": -0.20,
        "-10%": -0.10,
        "-5%": -0.05,
        "+5%": 0.05,
        "+10%": 0.10,
        "+20%": 0.20,
    }
    finite_keys = {
        "-20%": "-20_percent",
        "-10%": "-10_percent",
        "-5%": "-5_percent",
        "+5%": "+5_percent",
        "+10%": "+10_percent",
        "+20%": "+20_percent",
    }

    curve: dict[str, Any] = {}
    for label, frac in frac_map.items():
        central = qcd_delta(df, frac)
        leave = [qcd_delta(df[df.beta != b], frac) for b in betas]
        beta_jk = jackknife_sigma(leave)
        fit_sigma = float(finite["light_quark_mass_response"][finite_keys[label]]["sigma_MeV"])
        conservative = math.sqrt(fit_sigma * fit_sigma + beta_jk * beta_jk)
        curve[label] = {
            "delta_MN_MeV": central,
            "fit_sigma_MeV_from_frozen_fit": fit_sigma,
            "beta_block_jackknife_sigma_MeV": beta_jk,
            "conservative_sigma_MeV": conservative,
            "leave_one_beta_range_MeV": [min(leave), max(leave)],
        }

    # Numerical checks use the published precision, not hidden extra digits.
    checks = []
    checks.append(abs(sigma_central - float(promoted["associated_local_response"]["Sigma_ud_fixed_t0_MeV"]["central"])) < 0.02)
    checks.append(abs(kud - float(promoted["associated_local_response"]["K_ud_dlnMN_dlnmud_fixed_t0"]["central"])) < 2e-5)
    for label, rec in curve.items():
        ref = promoted["prediction_curve"][label]
        checks.append(abs(rec["delta_MN_MeV"] - float(ref["delta_MN_MeV"])) < 0.006)
        checks.append(abs(rec["conservative_sigma_MeV"] - float(ref["conservative_sigma_MeV"])) < 0.006)
        rr = ref["leave_one_beta_range_MeV"]
        checks.append(abs(rec["leave_one_beta_range_MeV"][0] - float(rr[0])) < 0.006)
        checks.append(abs(rec["leave_one_beta_range_MeV"][1] - float(rr[1])) < 0.006)

    signs_ok = all((frac_map[k] > 0 and v["delta_MN_MeV"] > 0) or (frac_map[k] < 0 and v["delta_MN_MeV"] < 0) for k, v in curve.items())
    checks.append(signs_ok)

    return {
        "verdict": "PASS" if all(checks) else "FAIL",
        "n_ensembles": int(len(df)),
        "Sigma_ud_fixed_t0_MeV": sigma_central,
        "K_ud_fixed_t0": kud,
        "leave_one_beta_Sigma_ud_MeV": sigma_leave,
        "beta_block_jackknife_sigma_MeV": sigma_beta_jk,
        "curve": curve,
        "sign_and_monotonic_direction": signs_ok,
        "claim_boundary": "Fixed-t0 RGI light-quark response; not promoted as a physical sigma term.",
    }


def reproduce_z6() -> dict[str, Any]:
    row28 = load_json(ROOT / "results" / "row28_proton_stability" / "result.json")
    ref = load_json(SUPPORT / "UTFANSWF_EXACT_Z6_TRINUCLEON_CHANNEL_PREDICTION_FROZEN_2026-08-19.json")
    charges = row28["residual_Z6"]["charges"]
    q_expected = {"Q": 0, "Uc": 5, "Dc": 1, "L": 2, "Ec": 5, "Nc": 3, "Hu": 1, "Hd": 5}

    # EWSB-safe physical selection rule from the frozen exact-Z6 construction:
    # q_phys = -DeltaB + 3 DeltaL (mod 6).  For DeltaB=1,2 there is no
    # integer DeltaL solution; DeltaB=3 is first allowed with odd DeltaL.
    allowed = {}
    for db in (1, 2, 3):
        sols = [dl for dl in range(-5, 6) if (-db + 3 * dl) % 6 == 0]
        allowed[db] = sols
    first_allowed = next(db for db in (1, 2, 3) if allowed[db])
    # Minimal DeltaB=3 local fermionic operator: 9 quark fields + 1 lepton field.
    operator_dimension = (9 + 1) * 3 / 2
    suppression_power = int(operator_dimension - 4)

    checks = [
        charges == q_expected,
        row28["structural_gate"]["checks"]["dangerous_DeltaB1_DeltaB2_operators_forbidden"] is True,
        float(row28["proton_stability_statement"]["isolated_proton_width_under_exact_Z6"]) == 0.0,
        allowed[1] == [],
        allowed[2] == [],
        first_allowed == 3,
        int(operator_dimension) == int(ref["predictions"]["lowest_operator_dimension"]),
        suppression_power == 11,
        ref["not_predicted"]["trinucleon_lifetime"] == "OPEN",
    ]
    return {
        "verdict": "PASS" if all(checks) else "FAIL",
        "charges_mod6": charges,
        "selection_solutions_DeltaL": {str(k): v for k, v in allowed.items()},
        "DeltaB1": "FORBIDDEN" if not allowed[1] else "ALLOWED",
        "DeltaB2": "FORBIDDEN" if not allowed[2] else "ALLOWED",
        "first_allowed_baryon_violation": f"DeltaB={first_allowed}",
        "lowest_operator_dimension": int(operator_dimension),
        "minimum_EFT_suppression_power": f"Lambda^-{suppression_power}",
        "isolated_DeltaB1_proton_width": 0.0,
        "trinucleon_lifetime": "OPEN",
        "leading_channel_classes": ref["predictions"]["leading_channel_classes"],
    }


def run_piecewise_inverse(alpha_inv0: np.ndarray, mu: float, low_is_doublet: bool) -> np.ndarray:
    MZ = 91.1876
    low = 1816.11748
    mid_scale = 2568.37797
    high = 3632.23496
    b = np.array([41 / 10, -19 / 6, -7.0], dtype=float)
    d_doublet = np.array([2 / 5, 2.0, 1.0], dtype=float)
    d_singlet = np.array([8 / 5, 0.0, 1.0], dtype=float)
    d_middle = np.array([1 / 2, 13 / 6, 2.0], dtype=float)
    inv = np.asarray(alpha_inv0, dtype=float).copy()
    last = MZ
    steps = [
        (low, d_doublet if low_is_doublet else d_singlet),
        (mid_scale, d_middle),
        (high, d_singlet if low_is_doublet else d_doublet),
    ]
    for threshold, add in steps:
        if mu <= threshold:
            inv -= b / (2 * math.pi) * math.log(mu / last)
            return inv
        inv -= b / (2 * math.pi) * math.log(threshold / last)
        b = b + add
        last = threshold
    inv -= b / (2 * math.pi) * math.log(mu / last)
    return inv


def alpha_em_inverse(inv: np.ndarray) -> float:
    return float((5 / 3) * inv[0] + inv[1])


def sin2theta(inv: np.ndarray) -> float:
    a1 = 1.0 / inv[0]
    a2 = 1.0 / inv[1]
    return float(((3 / 5) * a1) / (((3 / 5) * a1) + a2))


def reproduce_psn() -> dict[str, Any]:
    ref = load_json(SUPPORT / "UTFANSWF_PSN_V2_SPLIT_SPECTRUM_ALPHA_RUNNING_FROZEN_2026-08-19.json")
    kink = load_json(SUPPORT / "UTFANSWF_PSN_V2_GAUGE_KINK_SIGNATURE_FROZEN_2026-08-19.json")

    alpha1_inv_MZ = 59.0814505789
    alpha2_inv_MZ = 29.6485459418
    alpha3_inv_MZ = 1.0 / 0.1184
    inv0 = np.array([alpha1_inv_MZ, alpha2_inv_MZ, alpha3_inv_MZ], dtype=float)

    points = {
        "250_GeV": 250.0,
        "1_TeV": 1000.0,
        "2_TeV": 2000.0,
        "T_2p568378_TeV": 2568.37797,
        "3_TeV": 3000.0,
        "3p632235_TeV": 3632.23496,
        "10_TeV": 10000.0,
        "100_TeV": 100000.0,
    }
    em = {}
    alpha3 = {}
    for label, mu in points.items():
        inv_a = run_piecewise_inverse(inv0, mu, True)
        inv_b = run_piecewise_inverse(inv0, mu, False)
        # Exact invariance expected for alpha_em and alpha3 under D/S swap.
        em[label] = alpha_em_inverse(inv_a)
        alpha3[label] = float(inv_a[2])
        if abs(alpha_em_inverse(inv_a) - alpha_em_inverse(inv_b)) > 1e-12:
            raise AssertionError("D/S alpha_em swap invariance failed")
        if abs(float(inv_a[2]) - float(inv_b[2])) > 1e-12:
            raise AssertionError("D/S alpha3 swap invariance failed")

    weak = {}
    for label, mu in {"3_TeV": 3000.0, "10_TeV": 10000.0, "100_TeV": 100000.0}.items():
        vals = sorted([sin2theta(run_piecewise_inverse(inv0, mu, True)), sin2theta(run_piecewise_inverse(inv0, mu, False))])
        weak[label] = vals

    Bem = [11 / 3, 19 / 3, 28 / 3, 12.0]
    slopes = [-x / (2 * math.pi) for x in Bem]
    jumps = [slopes[i + 1] - slopes[i] for i in range(3)]
    d_weight = (5 / 3) * (2 / 5) + 2
    s_weight = (5 / 3) * (8 / 5) + 0

    checks = []
    for label, val in em.items():
        checks.append(abs(val - float(ref["predictions_alpha_em_inverse"][label])) < 6e-7)
    for label, vals in weak.items():
        rr = ref["weak_angle_discrete_envelope"][label]
        checks.extend([abs(vals[0] - float(rr[0])) < 6e-10, abs(vals[1] - float(rr[1])) < 6e-10])
    checks.extend([
        abs(d_weight - s_weight) < 1e-15,
        all(abs(jumps[i] - float(kink["alpha_em"]["slope_jumps"][i])) < 6e-10 for i in range(3)),
        abs(jumps[0] - jumps[2]) < 1e-15,
    ])

    return {
        "verdict": "PASS" if all(checks) else "FAIL",
        "thresholds_GeV": [1816.11748, 2568.37797, 3632.23496],
        "alpha_em_inverse": em,
        "alpha3_inverse": alpha3,
        "weak_angle_discrete_envelope": weak,
        "B_em_segments": Bem,
        "d_alpha_em_inverse_d_lnmu": slopes,
        "slope_jumps": jumps,
        "D_S_EM_weight_doublet": d_weight,
        "D_S_EM_weight_singlet": s_weight,
        "D_S_swap_alpha_em_and_alpha3_invariant": True,
        "claim_boundary": "One-loop split-threshold propagation of the frozen Aug-18 PSN-v2 chain; not an independent re-derivation of the upstream MZ normalization.",
    }


def main() -> None:
    qcd = reproduce_qcd()
    z6 = reproduce_z6()
    psn = reproduce_psn()
    passed = all(x["verdict"] == "PASS" for x in (qcd, z6, psn))
    out = {
        "gate": "new_predictions_2026_08_19_reproduction",
        "context": "POST_RELEASE_RESEARCH",
        "tier": "Auxiliary prediction reproduction",
        "blocking": True,
        "verdict": "PASS" if passed else "FAIL",
        "prediction_groups": {
            "NP01_fixed_t0_light_quark_nucleon_response": qcd,
            "NP02_exact_Z6_baryon_violation": z6,
            "NP03_PSN_v2_split_spectrum_gauge_running": psn,
        },
        "demoted_results_are_not_repromoted": True,
        "notes": "Reproduces the three authoritative Aug-19 prediction groups from the bundled v23 inputs/support artifacts. Canonical 18-gate numbering is unchanged.",
    }
    (OUTDIR / "result.json").write_text(json.dumps(out, indent=2, sort_keys=True) + "\n", encoding="utf-8")
    print(json.dumps(out, indent=2, sort_keys=True))
    if not passed:
        raise SystemExit(1)


if __name__ == "__main__":
    main()
