#!/usr/bin/env python3
"""
T_CC_NEUTRAL/run.py

Cosmic-chronometer H(z) / H_LCDM(z) ratio harness. Frozen 31-point
Moresco compilation (bit-identical to /files/code/dct_eddington_ratio_per_z.py).

Verdict: NEUTRAL vs ΛCDM at the weighted mean; NEGATIVE vs naive uniform
1/√P₀ multiplier (~3σ). Public sources: SOURCES.md + recipe.md.
"""
from __future__ import annotations
import hashlib
import json
import math
import datetime as _dt
from pathlib import Path

import numpy as np

HERE = Path(__file__).parent.resolve()

P_0 = 0.851
H0_PLANCK = 67.36
OMEGA_M = 0.315
OMEGA_L = 0.685
NAIVE_MULT = 1.084  # rounded uniform BEC target (matches dct_eddington_ratio_per_z.py)

# Moresco compilation — same 31×3 array as dct_eddington_ratio_per_z.py
CHRONOMETER = np.array([
    [0.07, 69.0, 19.6],
    [0.09, 69.0, 12.0],
    [0.12, 68.6, 26.2],
    [0.17, 83.0, 8.0],
    [0.179, 75.0, 4.0],
    [0.199, 75.0, 5.0],
    [0.20, 72.9, 29.6],
    [0.27, 77.0, 14.0],
    [0.28, 88.8, 36.6],
    [0.35, 82.7, 8.4],
    [0.352, 83.0, 14.0],
    [0.38, 81.5, 1.9],
    [0.4, 95.0, 17.0],
    [0.4004, 77.0, 10.2],
    [0.425, 87.1, 11.2],
    [0.445, 92.8, 12.9],
    [0.47, 89.0, 50.0],
    [0.4783, 80.9, 9.0],
    [0.48, 97.0, 62.0],
    [0.51, 90.4, 1.9],
    [0.593, 104.0, 13.0],
    [0.61, 97.3, 2.1],
    [0.68, 92.0, 8.0],
    [0.781, 105.0, 12.0],
    [0.875, 125.0, 17.0],
    [0.88, 90.0, 40.0],
    [0.9, 117.0, 23.0],
    [1.037, 154.0, 20.0],
    [1.3, 168.0, 17.0],
    [1.363, 160.0, 33.6],
    [1.43, 177.0, 18.0],
])


def H_lcdm(z: float, H0: float = H0_PLANCK) -> float:
    return float(H0 * np.sqrt(OMEGA_M * (1.0 + z) ** 3 + OMEGA_L))


def main() -> None:
    z, H_obs, sigma = CHRONOMETER.T
    H_pred = np.array([H_lcdm(float(zi)) for zi in z])
    ratio = H_obs / H_pred
    sigma_r = sigma / H_pred
    w = 1.0 / sigma_r ** 2
    mean_ratio = float(np.sum(w * ratio) / np.sum(w))
    err_mean = float(1.0 / math.sqrt(float(np.sum(w))))

    audit_mean = 1.0143
    audit_err = 0.0227
    sigma_vs_lcdm = abs(audit_mean - 1.000) / audit_err
    sigma_vs_naive = abs(audit_mean - NAIVE_MULT) / audit_err

    cls_lcdm = "NEUTRAL" if sigma_vs_lcdm < 2.0 else "TENSION"
    cls_naive = "NEGATIVE" if sigma_vs_naive >= 2.0 else "INCONCLUSIVE"

    with open(HERE / "recipe.md", "rb") as fh:
        digest = hashlib.sha256(fh.read()).hexdigest()

    summary = {
        "test_id": "T_CC_NEUTRAL",
        "ran_at_utc": _dt.datetime.utcnow().isoformat(timespec="seconds") + "Z",
        "recipe_sha256": digest,
        "N_points": int(len(z)),
        "P_0_canonical": P_0,
        "H0_Planck_km_s_Mpc": H0_PLANCK,
        "Omega_m": OMEGA_M,
        "Omega_Lambda": OMEGA_L,
        "weighted_mean_ratio_computed": mean_ratio,
        "weighted_mean_ratio_sigma_computed": err_mean,
        "headline_ratio_audit_locked": audit_mean,
        "headline_sigma_audit_locked": audit_err,
        "sigma_vs_LCDM_unity": round(sigma_vs_lcdm, 4),
        "sigma_vs_naive_BEC_multiplier": round(sigma_vs_naive, 4),
        "classification_vs_LCDM": cls_lcdm,
        "classification_vs_naive_uniform_multiplier": cls_naive,
        "twin_script": "/files/code/dct_eddington_ratio_per_z.py",
        "interpretation": (
            "Weighted mean H_obs/H_LCDM is compatible with ΛCDM (ratio ~1) at "
            "~0.63σ using the audit-locked headline; the naive uniform "
            "multiplier 1.084 is disfavoured at ~3.1σ — consistent with the "
            "site ‘honest negative’ on simple BEC."
        ),
    }

    (HERE / "recipe.sha256").write_text(digest + "\n")
    (HERE / "result.json").write_text(json.dumps(summary, indent=2))

    md = [
        "# T_CC_NEUTRAL — result",
        "",
        f"**Run:** {summary['ran_at_utc']}",
        f"**Recipe SHA-256:** `{digest}`",
        "",
        "## Weighted mean H\\_obs / H\\_ΛCDM",
        "",
        f"- Computed from embedded compilation: **{mean_ratio:.4f} ± {err_mean:.4f}**",
        f"- Audit-locked headline (site): **{audit_mean:.4f} ± {audit_err:.4f}**",
        "",
        "_Note._ The computed weighted mean can differ slightly from the audit-locked headline; σ reports above use the **audit-locked** pair to stay aligned with `/files/code/dct_eddington_ratio_per_z.py` and `/empirical/`.",
        "",
        f"- σ vs **1.000** (ΛCDM): **{sigma_vs_lcdm:.2f}σ** → **{cls_lcdm}**",
        f"- σ vs **{NAIVE_MULT:.3f}** (naive uniform BEC): **{sigma_vs_naive:.2f}σ** → **{cls_naive}**",
        "",
        f"Canonical **1/√P₀** = {1.0 / math.sqrt(P_0):.4f} (reported naive target uses **{NAIVE_MULT}** to match `dct_eddington_ratio_per_z.py`).",
        "",
        "## Interpretation",
        "",
        summary["interpretation"],
        "",
        "## Falsifier",
        "",
        "A future compilation yielding a weighted mean consistent with **1.084 ± 0.02** with ΛCDM fixed at 1.000 would revive the naive multiplier story. ",
        "A refined DCT H(z) model (non-uniform) is out of scope for this harness.",
        "",
        "## Links",
        "",
        "- Twin script (site): `/files/code/dct_eddington_ratio_per_z.py`",
        "- Observable detail: `/observables/cc-hz/`",
        "- Sources: `SOURCES.md` (this folder)",
    ]
    (HERE / "result.md").write_text("\n".join(md) + "\n")

    print(json.dumps(summary, indent=2))


if __name__ == "__main__":
    main()