#!/usr/bin/env python3
"""
verify_paper_sparc.py — re-derives every number quoted in
sparc_one_scale_halo.tex from the repository's pipeline outputs
(Science/HiddenBranch/output/*.json) and physical constants.
Exit 0 iff all pass.

  S1  Amplitude section            — checks 1-6
  S2  Shape section                — checks 7-12
  S3  Unit system                  — checks 13-16
  S4  Falsifier numbers            — checks 17-18
"""
import sys, os, json, math
import mpmath as mp

mp.mp.dps = 30
PASS = FAIL = 0
def check(n, desc, cond):
    global PASS, FAIL
    ok = bool(cond); PASS += ok; FAIL += (not ok)
    print(f"  [{'PASS' if ok else 'FAIL'}] {n:>2}. {desc}")

HERE = os.path.dirname(os.path.abspath(__file__))
ROOT = os.path.dirname(os.path.dirname(os.path.dirname(HERE)))
OUT = os.path.join(ROOT, "Science", "HiddenBranch", "output")
J = lambda f: json.load(open(os.path.join(OUT, f)))
v5 = J("sparc_universality_v5_results.json")
v6 = J("sparc_calibration_v6_results.json")
v7 = J("sparc_shape_v7_results.json")
v8 = J("sparc_backreaction_v8_results.json")
v9 = J("sparc_masked_v9_results.json")

print("S1  Amplitude (Sec. 4)")
check(1, "universality slope 0.029 +/- 0.076 (n=129)",
      v5["slope_log_aeff_vs_log_Rd"] == 0.029 and v5["slope_se"] == 0.076
      and v5["n_galaxies"] == 129)
check(2, "rival (core tracks baryonic radius) excluded at 13.5 sigma",
      abs(v5["model_R_rc_eq_Rb_over_pi"]["sigma_away"] - 13.49) < 0.02)
check(3, "BTFR scatter 0.285 dex", v5["btfr_scatter_dex"] == 0.285)
g = v6["gap1"]
F = float
check(4, "calibration: gas-dominated (n=51, U=0.5) a0 = 1.172e-10 = "
         "0.976 x MOND",
      F(g["upsilon_0.5"]["gas_dominated"]["a0"]) == 1.172e-10 and
      g["upsilon_0.5"]["gas_dominated"]["n"] == 51 and
      F(g["vs_mond_1.2e-10"]) == 0.976)
check(5, "systematic envelope: full 1.67/1.33e-10 (U=0.5/0.7), "
         "gas U=0.7 1.05e-10",
      abs(F(g["upsilon_0.5"]["full_sample"]["a0"]) - 1.671e-10) < 1e-13 and
      abs(F(g["upsilon_0.7"]["full_sample"]["a0"]) - 1.330e-10) < 1e-13 and
      abs(F(g["upsilon_0.7"]["gas_dominated"]["a0"]) - 1.046e-10) < 1e-13)
cc = v6["gap2"]["candidate_couplings"]
check(6, "coupling table: a^3/pi^2 -> 0.979 kpc; a^3 -> 9.7; "
         "a^(5/2)/pi^2 -> 11.5; a^(7/2)/pi^2 -> 0.08; a^3/pi -> 3.1",
      abs(F(cc["alpha^3/pi^2 (Codex v4)"]["implied_rc_kpc"]) - 0.9789) < 1e-3
      and abs(F(cc["alpha^3 (bare three-layer)"]["implied_rc_kpc"]) - 9.66)
      < 0.01
      and abs(F(cc["alpha^(5/2)/pi^2"]["implied_rc_kpc"]) - 11.46) < 0.01
      and abs(F(cc["alpha^(7/2)/pi^2"]["implied_rc_kpc"]) - 0.0836) < 1e-3
      and abs(F(cc["alpha^3/pi (edge-normalized)"]["implied_rc_kpc"]) - 3.075)
      < 0.01)

print("S2  Shape (Sec. 5)")
check(7, "naive fit: 125 constrained of 138; median 1.26 kpc, "
         "16-84%% [0.56, 3.15], scatter 0.375 dex",
      v7["n_constrained"] == 125 and v7["n_galaxies_fit"] == 138 and
      abs(v7["rc_kpc"]["median"] - 1.255) < 1e-3 and
      abs(v7["rc_kpc"]["p16"] - 0.561) < 1e-3 and
      abs(v7["rc_kpc"]["p84"] - 3.153) < 1e-3 and
      v7["universality_scatter_dex"] == 0.375)
check(8, "naive slope 0.623 +/- 0.082 (7.6 sigma)",
      v7["rc_vs_Rdisk_logslope"]["value"] == 0.623 and
      v7["rc_vs_Rdisk_logslope"]["se"] == 0.082 and
      abs(0.623/0.082 - 7.6) < 0.05)
check(9, "split: clean (60) median 0.889 kpc; baryon-dominated (65) "
         "2.232 kpc",
      v8["n"]["clean"] == 60 and v8["n"]["glarey"] == 65 and
      v8["rc_median_kpc"]["clean"] == 0.889 and
      v8["rc_median_kpc"]["glarey"] == 2.232)
check(10, "pure-glare model rejected within clean half: residual slope "
          "0.426 +/- 0.080 (5.3 sigma)",
      "0.426" in v8["v8b_model_tests"]["glare_only"] and
      "5.3 sigma" in v8["v8b_model_tests"]["glare_only"])
check(11, "masked ladder: 0.618+/-0.080 (120, 1.26) -> 0.363+/-0.092 "
          "(76, 1.00) -> 0.254+/-0.110 (46, 0.89)",
      v9["n"] == 76 and v9["slope"] == [0.363, 0.092] and
      abs(v9["median_rc"] - 0.997) < 0.001)
check(12, "candidate match in shapes: a^3/pi^2 median/pred = 1.282 "
          "(rivals 0.13-15x off)",
      float(v7["candidate_match_med_over_pred"]["alpha^3/pi^2"]) == 1.282 and
      float(v7["candidate_match_med_over_pred"]["alpha^3"]) == 0.13 and
      float(v7["candidate_match_med_over_pred"]["alpha^(7/2)/pi^2"]) == 15.012)

print("S3  Unit system (Sec. 6)")
PI = mp.pi
OM = 4*PI**3 + PI**2 + PI
AL = 1/OM
a0 = mp.mpf("1.172e-10"); c = mp.mpf("2.99792458e8")
Lh = AL**3*c**2/(PI**2*a0)
kpc = mp.mpf("3.0857e19"); yr = mp.mpf(86400)*mp.mpf("365.25")
tau = Lh/c
T = (PI/AL)*tau
check(13, "L_h = alpha^3 c^2/(pi^2 a0) = %s kpc (0.98)" % mp.nstr(Lh/kpc, 4),
      abs(Lh/kpc - mp.mpf("0.9785")) < 0.001)
check(14, "tau_u = L_h/c = %s yr (3.19e3); T = (pi/alpha) tau = %s Myr "
          "(1.374)" % (mp.nstr(tau/yr, 5), mp.nstr(T/yr/1e6, 5)),
      abs(tau/yr - 3191.4) < 1 and abs(T/yr/1e6 - mp.mpf("1.3739")) < 0.001)
check(15, "cT/L_h = pi/alpha = %s (430.51)" % mp.nstr(c*T/Lh, 8),
      abs(c*T/Lh - PI/AL) < mp.mpf(10)**-20)
check(16, "a0 T^2/L_h = alpha exactly: %s vs %s"
         % (mp.nstr(a0*T**2/Lh, 10), mp.nstr(AL, 10)),
      abs(a0*T**2/Lh - AL) < mp.mpf(10)**-25)

print("S4  Falsifiers (Sec. 7)")
check(17, "F1 bound: strict-mask slope 0.254 +/- 0.110 (2.3 sigma)",
      abs(0.254/0.110 - 2.31) < 0.01)
check(18, "F3 window: gas-dominated calibrations 1.046-1.172e-10 inside "
          "the stated (1.05-1.35)e-10 stability window",
      1.04e-10 <= float(g["upsilon_0.7"]["gas_dominated"]["a0"]) <= 1.35e-10 and
      1.05e-10 <= float(g["calibrated_a0"]) <= 1.35e-10)

print(f"\n{'='*60}\nRESULT: {PASS} PASS / {FAIL} FAIL")
sys.exit(0 if FAIL == 0 else 1)