#!/usr/bin/env python3
"""
UTFANSWF v21 — Row 20 Thin-Slice Attribution Diagnostic

Native reproducibility extension only:
- reads the existing Row 20 dataset/theory/covariance files;
- reproduces the locked Row 20 block chi^2 accounting;
- reports block-level and slice-level residual attribution;
- adds an optional SH0ES overlay diagnostic;
- changes no physics, thresholds, datasets, parameters, or gate pass/fail rules.

Run from project root:
    python tools/row20_thin_slice_attribution.py
"""
from __future__ import annotations

import csv
import json
import math
import os
from datetime import datetime, timezone
from pathlib import Path
from typing import Dict, List, Tuple

import numpy as np

ROOT = Path(__file__).resolve().parents[1]
DATA = ROOT / "datasets"
OUTDIR = ROOT / "results" / "row20_h0_attribution"
OUTDIR.mkdir(parents=True, exist_ok=True)

THRESHOLD = 1.25
SH0ES_H0 = 73.04
SH0ES_SIGMA = 1.04


def load_vec(name: str) -> np.ndarray:
    return np.loadtxt(DATA / name, delimiter=",", dtype=float).reshape(-1)


def load_cov(name: str) -> np.ndarray:
    C = np.loadtxt(DATA / name, delimiter=",", dtype=float)
    C = np.atleast_2d(C).astype(float)
    return C


def chi2(d: np.ndarray, t: np.ndarray, C: np.ndarray) -> float:
    r = d - t
    return float(r @ np.linalg.inv(C) @ r)


def sn_fit(d: np.ndarray, t: np.ndarray, C: np.ndarray) -> Tuple[float, float, np.ndarray]:
    """Return chi2 after fitting only the existing standard SN additive offset."""
    invC = np.linalg.inv(C)
    ones = np.ones_like(d)
    raw = d - t
    denom = float(ones @ invC @ ones)
    if denom == 0:
        raise ValueError("SN offset denominator is zero")
    m_fit = float(ones @ invC @ raw) / denom
    r_post = d - (t + m_fit)
    x2 = float(r_post @ invC @ r_post)
    return x2, m_fit, r_post


def block_compute(tag: str, data_fn: str, theory_fn: str, cov_fn: str, labels: List[str], sn: bool = False) -> Dict:
    d, t, C = load_vec(data_fn), load_vec(theory_fn), load_cov(cov_fn)
    if sn:
        x2, m_fit, r = sn_fit(d, t, C)
        dof = max(len(d) - 1, 1)
        nuisance = {"SN_M_fit_existing_offset": m_fit}
    else:
        r = d - t
        x2 = chi2(d, t, C)
        dof = len(d)
        nuisance = {}
    sig = np.sqrt(np.maximum(np.diag(C), 0.0))
    diag_pull = np.divide(r, sig, out=np.full_like(r, np.nan, dtype=float), where=sig > 0)

    # Cholesky-whitened residuals exactly sum to chi2 for positive-definite C.
    try:
        L = np.linalg.cholesky(C)
        white = np.linalg.solve(L, r)
        white_contrib = white * white
        white_note = "Cholesky-whitened contributions; order-dependent for correlated covariance but sum exactly to block chi2."
    except np.linalg.LinAlgError:
        invC = np.linalg.inv(C)
        white = np.full_like(r, np.nan)
        white_contrib = np.full_like(r, np.nan)
        white_note = "Covariance not Cholesky factorizable."

    # Leave-one-slice attribution using covariance submatrices.
    loo = []
    for i, lab in enumerate(labels):
        if len(d) <= 1:
            continue
        mask = np.ones(len(d), dtype=bool)
        mask[i] = False
        d2, t2, C2 = d[mask], t[mask], C[np.ix_(mask, mask)]
        if sn:
            x2_without, _, _ = sn_fit(d2, t2, C2)
            dof_without = max(len(d2) - 1, 1)
        else:
            x2_without = chi2(d2, t2, C2)
            dof_without = len(d2)
        loo.append({
            "slice": lab,
            "chi2_without_slice": x2_without,
            "dof_without_slice": dof_without,
            "delta_chi2_removed": x2 - x2_without,
            "chi2_per_dof_without_slice": x2_without / max(dof_without, 1),
        })

    node_rows = []
    for i, lab in enumerate(labels):
        node_rows.append({
            "slice": lab,
            "data": float(d[i]),
            "theory": float(t[i]),
            "residual_after_gate_handling": float(r[i]),
            "sigma_diag": float(sig[i]),
            "diag_pull_sigma": float(diag_pull[i]),
            "white_residual": float(white[i]) if np.isfinite(white[i]) else None,
            "white_chi2_contribution": float(white_contrib[i]) if np.isfinite(white_contrib[i]) else None,
        })

    return {
        "tag": tag,
        "chi2": x2,
        "dof": dof,
        "chi2_per_dof": x2 / max(dof, 1),
        "n": len(d),
        "nuisance": nuisance,
        "nodes": node_rows,
        "leave_one_slice_out": sorted(loo, key=lambda x: x["delta_chi2_removed"], reverse=True),
        "white_note": white_note,
    }


def fmt(x: float) -> str:
    return f"{x:.6f}"


def main() -> None:
    params = json.loads((DATA / "cosmo_params.json").read_text(encoding="utf-8"))
    h0 = float(params.get("H0", 67.4))

    blocks = [
        block_compute("CMB", "cmb_data.csv", "cmb_theory.csv", "cmb_cov.csv", ["R", "ell_A", "omega_b_h2"]),
        block_compute("BAO_DV", "bao_data_DV.csv", "bao_theory_DV.csv", "bao_cov_DV.csv", ["D_V/r_d z=0.295"]),
        block_compute("BAO_DM_H", "bao_data_DM_H.csv", "bao_theory_DM_H.csv", "bao_cov_DM_H.csv", ["D_M z=0.510", "H z=0.510", "D_M z=0.706", "H z=0.706"]),
        block_compute("SN", "sn_data.csv", "sn_theory.csv", "sn_cov.csv", ["SN mu z=0.10", "SN mu z=0.35", "SN mu z=0.70"], sn=True),
        block_compute("GW", "gw_data.csv", "gw_theory.csv", "gw_cov.csv", [f"GW proxy H0 node {i+1}" for i in range(len(load_vec("gw_data.csv")))])
    ]

    by_tag = {b["tag"]: b for b in blocks}
    x2_cosmo_no_gw = sum(by_tag[k]["chi2"] for k in ["CMB", "BAO_DV", "BAO_DM_H", "SN"])
    dof_cosmo_no_gw = sum(by_tag[k]["dof"] for k in ["CMB", "BAO_DV", "BAO_DM_H", "SN"])
    x2_full = x2_cosmo_no_gw + by_tag["GW"]["chi2"]
    dof_full = dof_cosmo_no_gw + by_tag["GW"]["dof"]

    sh0es_pull = (h0 - SH0ES_H0) / SH0ES_SIGMA
    sh0es_chi2 = sh0es_pull * sh0es_pull

    summary = {
        "signature": "ROW20_THIN_SLICE_ATTRIBUTION_v1_NATIVE_DIAGNOSTIC",
        "context": "DIAGNOSTIC_NON_GATING",
        "rules": {
            "physics_changed": False,
            "thresholds_changed": False,
            "new_parameters_added": False,
            "datasets_changed": False,
            "gate_verdict_changed": False,
        },
        "timestamp_utc": datetime.now(timezone.utc).isoformat(timespec="seconds").replace("+00:00", "Z"),
        "H0_used_by_v21": h0,
        "threshold_reference_chi2_per_dof": THRESHOLD,
        "block_summary": [
            {"tag": b["tag"], "chi2": b["chi2"], "dof": b["dof"], "chi2_per_dof": b["chi2_per_dof"]}
            for b in blocks
        ],
        "totals": {
            "CMB_BAO_SN_no_GW": {"chi2": x2_cosmo_no_gw, "dof": dof_cosmo_no_gw, "chi2_per_dof": x2_cosmo_no_gw / dof_cosmo_no_gw, "verdict_against_reference_threshold": "PASS" if x2_cosmo_no_gw / dof_cosmo_no_gw <= THRESHOLD else "FAIL"},
            "locked_v21_with_GW": {"chi2": x2_full, "dof": dof_full, "chi2_per_dof": x2_full / dof_full, "verdict_against_reference_threshold": "PASS" if x2_full / dof_full <= THRESHOLD else "FAIL"},
            "SH0ES_overlay": {"H0_SH0ES": SH0ES_H0, "sigma": SH0ES_SIGMA, "pull_sigma": sh0es_pull, "chi2": sh0es_chi2},
            "CMB_BAO_SN_plus_SH0ES_no_GW": {"chi2": x2_cosmo_no_gw + sh0es_chi2, "dof": dof_cosmo_no_gw + 1, "chi2_per_dof": (x2_cosmo_no_gw + sh0es_chi2)/(dof_cosmo_no_gw + 1), "verdict_against_reference_threshold": "PASS" if (x2_cosmo_no_gw + sh0es_chi2)/(dof_cosmo_no_gw + 1) <= THRESHOLD else "FAIL"},
            "locked_v21_plus_SH0ES_overlay_with_GW": {"chi2": x2_full + sh0es_chi2, "dof": dof_full + 1, "chi2_per_dof": (x2_full + sh0es_chi2)/(dof_full + 1), "verdict_against_reference_threshold": "PASS" if (x2_full + sh0es_chi2)/(dof_full + 1) <= THRESHOLD else "FAIL"},
        },
        "thin_slice_targets": [],
        "blocks": blocks,
        "interpretation": []
    }

    # Identify strongest leave-one-slice improvements across non-GW cosmology blocks.
    targets = []
    for b in blocks:
        if b["tag"] == "GW":
            continue
        for row in b["leave_one_slice_out"]:
            targets.append({"block": b["tag"], **row})
    targets.sort(key=lambda x: x["delta_chi2_removed"], reverse=True)
    summary["thin_slice_targets"] = targets[:10]

    # Conservative interpretation based only on diagnostics.
    summary["interpretation"] = [
        "The locked v21 Row 20 total remains PASS when the GW block is included, matching the released harness accounting.",
        "The CMB+BAO+SN subset without GW is stressed and fails against the reference threshold, so the H0 attribution cannot yet be reduced to SH0ES alone.",
        "The largest thin-slice residuals are localized rather than broad: the diagnostic ranks the specific CMB/SN/BAO nodes responsible for most chi2 reduction when removed.",
        "This is an attribution result only, not a resolution. A v21.1 candidate would require native reproducible resolution or formal attribution of these nodes without physics, threshold, dataset, or parameter changes.",
    ]

    (OUTDIR / "row20_thin_slice_attribution.json").write_text(json.dumps(summary, indent=2), encoding="utf-8")

    # Markdown report
    lines = []
    lines.append("# UTFANSWF v21 — Row 20 Thin-Slice Attribution Diagnostic")
    lines.append("")
    lines.append("**Context:** DIAGNOSTIC_NON_GATING. This extension changes no physics, thresholds, datasets, parameters, or gate verdicts.")
    lines.append("")
    lines.append("## Block summary")
    lines.append("")
    lines.append("| Block | chi2 | dof | chi2/dof |")
    lines.append("|---|---:|---:|---:|")
    for b in blocks:
        lines.append(f"| {b['tag']} | {fmt(b['chi2'])} | {b['dof']} | {fmt(b['chi2_per_dof'])} |")
    lines.append("")
    lines.append("## Totals")
    lines.append("")
    lines.append("| Test | chi2 | dof | chi2/dof | verdict vs 1.25 |")
    lines.append("|---|---:|---:|---:|---|")
    for name, t in summary["totals"].items():
        if "chi2_per_dof" in t:
            lines.append(f"| {name} | {fmt(t['chi2'])} | {t['dof']} | {fmt(t['chi2_per_dof'])} | {t['verdict_against_reference_threshold']} |")
        else:
            lines.append(f"| {name} | {fmt(t['chi2'])} | 1 | {fmt(t['chi2'])} | diagnostic pull {fmt(t['pull_sigma'])} sigma |")
    lines.append("")
    lines.append("## Strongest thin-slice targets")
    lines.append("")
    lines.append("| Rank | Block | Slice removed | delta chi2 removed | chi2/dof without slice |")
    lines.append("|---:|---|---|---:|---:|")
    for i, r in enumerate(summary["thin_slice_targets"][:10], 1):
        lines.append(f"| {i} | {r['block']} | {r['slice']} | {fmt(r['delta_chi2_removed'])} | {fmt(r['chi2_per_dof_without_slice'])} |")
    lines.append("")
    lines.append("## Node residuals")
    for b in blocks:
        if b["tag"] == "GW":
            continue
        lines.append("")
        lines.append(f"### {b['tag']}")
        if b.get("nuisance"):
            lines.append(f"Existing nuisance handling: `{b['nuisance']}`")
        lines.append("| Slice | data | theory | residual | diag pull sigma | whitened chi2 contribution |")
        lines.append("|---|---:|---:|---:|---:|---:|")
        for n in b["nodes"]:
            wc = n['white_chi2_contribution']
            lines.append(f"| {n['slice']} | {fmt(n['data'])} | {fmt(n['theory'])} | {fmt(n['residual_after_gate_handling'])} | {fmt(n['diag_pull_sigma'])} | {fmt(wc if wc is not None else float('nan'))} |")
    lines.append("")
    lines.append("## Interpretation")
    lines.extend([f"- {s}" for s in summary["interpretation"]])
    lines.append("")
    lines.append("## v21.1 candidate gate")
    lines.append("A v21.1 candidate remains blocked until the thin-slice residuals above receive a native reproducible resolution or formal attribution without changing physics, thresholds, datasets, or adding parameters.")
    (OUTDIR / "REPORT.md").write_text("\n".join(lines) + "\n", encoding="utf-8")

    print("Wrote:")
    print(OUTDIR / "row20_thin_slice_attribution.json")
    print(OUTDIR / "REPORT.md")
    print("\nTop thin-slice targets:")
    for i, r in enumerate(summary["thin_slice_targets"][:5], 1):
        print(f"{i}. {r['block']} {r['slice']}: delta chi2={r['delta_chi2_removed']:.6f}")


if __name__ == "__main__":
    main()
