"""
DECISIVE SN-LENSING TEST: DES-SN5YR (~1600 SNe inside the DES footprint) x DES Y3 foreground convergence.
DFD predicts Hubble-residual slope b>0 vs foreground kappa; LCDM cosmological constant predicts b~0/<=0.
Honest error = bootstrap. ~4x the Pantheon+ DES overlap -> the strongest existing-data test.
"""
import sys, numpy as np, healpy as hp
MAPF = sys.argv[1] if len(sys.argv)>1 else 'des_y3/wiener_tomo1.fits'
ZCUT = float(sys.argv[2]) if len(sys.argv)>2 else 0.2

# ---- DES-SN5YR ----
lines=[l for l in open('des_sn5yr/DES-Dovekie_Metadata.csv') if l.strip() and not l.startswith('#')]
hdr=lines[0].replace('VARNAMES:','').split(); idx={n:hdr.index(n) for n in hdr}
arr=np.array([r.split() for r in lines[1:] if r.startswith('SN') and len(r.split())>=len(hdr)+1])
g=lambda n: arr[:,idx[n]+1].astype(float); gs=lambda n: arr[:,idx[n]+1]
ra=g('HOST_RA'); dec=g('HOST_DEC'); z=g('zHD'); mures=g('MURES'); muerr=g('MUERR'); ids=gs('IDSURVEY')
sel=(ids=='10')&(ra>-900)&(dec>-900)&(z>ZCUT)&(z<0.9)&np.isfinite(mures)&(muerr>0)&(muerr<1)
ra,dec,z,res,err=ra[sel],dec[sel],z[sel],mures[sel],muerr[sel]
w=1/err**2; res=res-np.average(res,weights=w)
print(f"DES-SN5YR: {len(z)} DES SNe (z>{ZCUT}), <z>={z.mean():.3f}, residual rms={res.std():.3f} mag")

# ---- DES foreground convergence + footprint mask ----
k=hp.smoothing(hp.read_map(MAPF), fwhm=np.radians(0.25))
m=hp.read_map('des_y3/glimpse_mask.fits')
pix=hp.ang2pix(1024, np.radians(90-dec), np.radians(ra%360))
inside=(m[pix]>0.5)
kap=k[pix]
good=inside&np.isfinite(kap)
zz,rr2,ww,kk=z[good],res[good],w[good],kap[good]
print(f"map={MAPF.split('/')[-1]}  {good.sum()} SNe inside DES footprint")
if good.sum()<30: print("too few"); sys.exit()

fit=lambda r,kv,wv:(lambda kb:np.sum(wv*(kv-kb)*r)/np.sum(wv*(kv-kb)**2))(np.average(kv,weights=wv))
b=fit(rr2,kk,ww)
rng=np.random.default_rng(0); bs=np.array([fit(rr2[i],kk[i],ww[i]) for i in (rng.integers(0,len(zz),len(zz)) for _ in range(3000))])
sig=b/bs.std()
print(f"\n  SLOPE  delta-mu = b * kappa_fg :  b = {b:+.2f}   bootstrap sigma_b = {bs.std():.2f}")
print(f"  SIGNIFICANCE = {sig:+.2f} sigma   ->  {'DFD-sign (positive)' if b>0 else 'LCDM-sign (<=0)'}")
# rotation-null cross-check (rotate within footprint)
nl=[]
for s in range(300):
    p2=hp.ang2pix(1024,np.radians(90-dec[good]),np.radians((ra[good]+rng.uniform(30,330))%360))
    ins2=m[p2]>0.5
    if ins2.sum()<30: continue
    nl.append(fit(rr2[ins2],k[p2][ins2],ww[ins2]))
nl=np.array(nl); print(f"  rotation-null: {nl.mean():+.2f} +/- {nl.std():.2f}  -> null-sigma {b/nl.std():+.2f}")
