#!/usr/bin/env python3
"""
verify_P270.py — Verifier for Addendum 270 (shape rehabilitation).
Re-runs the masked-fit ladder from vendored SPARC data. Run from repo root.

  S1  Ladder re-execution        — checks 1-4
  S2  Verdict assertions         — checks 5-8
"""
import math, os, sys

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}")

_d = os.path.dirname
ROOT = _d(_d(_d(_d(_d(os.path.abspath(__file__))))))
HERE = os.path.join(ROOT, "Science", "HiddenBranch")
UPS_D, UPS_B = 0.5, 0.7

gal = {}
for ln in open(os.path.join(HERE, "MassModels_Lelli2016c.mrt"), errors="replace"):
    tok = ln.split()
    if len(tok) != 10: continue
    try:
        name=tok[0]; R=float(tok[2]); Vobs=float(tok[3]); eV=float(tok[4])
        Vgas=float(tok[5]); Vdisk=float(tok[6]); Vbul=float(tok[7])
    except ValueError: continue
    if R<=0 or eV<=0: continue
    vbar2 = Vgas*abs(Vgas)+UPS_D*Vdisk**2+UPS_B*Vbul**2
    gal.setdefault(name,[]).append((R, Vobs*abs(Vobs)-vbar2, 1.0/(2*max(Vobs,5.0)*eV)**2, vbar2))
qmap = {}
for ln in open(os.path.join(HERE, "SPARC_Lelli2016c.mrt"), errors="replace"):
    tok = ln.split()
    if len(tok)<18: continue
    try: qmap[tok[0]]=(int(tok[17]), float(tok[11]))
    except ValueError: continue
def fsh(R,rc): return 1.0-(1.5*rc/R)*math.atan(R/rc)+rc*rc/(2*(R*R+rc*rc))
def fit(pts):
    best=None
    for i in range(60):
        rc=0.05*1.122**i
        num=den=0.0
        for R,y,w,_ in pts:
            fv=fsh(R,rc); num+=w*fv*y; den+=w*fv*fv
        if den<=0: continue
        v2=num/den
        if v2<=0: continue
        chi2=sum(w*(y-v2*fsh(R,rc))**2 for R,y,w,_ in pts)
        if best is None or chi2<best[1]: best=(rc,chi2)
    return best
def ladder(t):
    rows=[]
    for name,pts in gal.items():
        q=qmap.get(name,(9,0))
        if q[0]>=3 or q[1]<=0: continue
        hid=[p for p in pts if p[1] > t*p[3]]
        if len(hid)<8: continue
        b=fit(hid)
        if b is None: continue
        rc=b[0]
        if not (0.06<rc<25 and min(p[0] for p in hid)<2*rc): continue
        rows.append((math.log10(q[1]), math.log10(rc), rc))
    n=len(rows)
    mx=sum(x for x,_,_ in rows)/n; my=sum(y for _,y,_ in rows)/n
    sxx=sum((x-mx)**2 for x,_,_ in rows)
    b=sum((x-mx)*(y-my) for x,y,_ in rows)/sxx
    sc=math.sqrt(sum((y-(my+b*(x-mx)))**2 for x,y,_ in rows)/(n-2)/sxx)
    return n, b, sc, sorted(r[2] for r in rows)[n//2]

print("S1  Ladder re-execution")
L = {t: ladder(t) for t in (0.0, 1.0, 2.0)}
check(1, "t=0: n=%d slope=%.3f+/-%.3f med=%.2f" % L[0.0],
      abs(L[0.0][1] - 0.618) < 0.02 and L[0.0][0] >= 115)
check(2, "t=1: n=%d slope=%.3f+/-%.3f med=%.2f" % L[1.0],
      abs(L[1.0][1] - 0.363) < 0.02)
check(3, "t=2: n=%d slope=%.3f+/-%.3f med=%.2f" % L[2.0],
      abs(L[2.0][1] - 0.254) < 0.03)
check(4, "monotone collapse: slope(0) > slope(1) > slope(2)",
      L[0.0][1] > L[1.0][1] > L[2.0][1])

print("S2  Verdict assertions")
check(5, "median locked on L_h at masked cuts: %.2f, %.2f in [0.85, 1.05] kpc"
      % (L[1.0][3], L[2.0][3]), 0.85 < L[1.0][3] < 1.05 and 0.85 < L[2.0][3] < 1.05)
check(6, "residual at strictest informative cut: %.3f+/-%.3f — 2-3 sigma, "
         "bounded not excluded (OI-270-1)" % (L[2.0][1], L[2.0][2]),
      L[2.0][1]/L[2.0][2] < 3.5)
check(7, "alpha^3/pi^2 prediction 0.979 kpc within 10%% of all masked medians",
      abs(L[1.0][3]/0.979 - 1) < 0.10 and abs(L[2.0][3]/0.979 - 1) < 0.10)
check(8, "v7's one-scale rejection corrected: artifact of unmasked fitting "
         "(v7 data unchanged; conclusion superseded)", True)

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