"""
POWER TEST + amplitude reality-check (the referee's 3 corrections).
(1) Is the BOSS kappa_fg actually in PHYSICAL convergence units? Compare its std to the expected
    foreground-convergence RMS (~0.01-0.02). If ours is far smaller, it is a SUPPRESSED proxy, NOT a
    physical convergence -> then only the SIGN test is meaningful and we must stop saying "fully calibrated".
(2) Power test: inject signals through the ACTUAL pipeline (the saved 648-SN kappa_fg dataset) and find the
    slope the test can DETECT at 3 sigma. Is +1.5+/-2.9 underpowered, or does it constrain a given amplitude?
(3) Resolve: calibrated -> constrains large amplitudes; proxy -> inconclusive sign test. Not both.
"""
import numpy as np
d=np.load("boss_data/kappa_fg_dataset.npz")
k=d['kappa_fg']; y=d['dmu']; e=d['err']; mass=d['hostmass']; ra=d['ra']; w=1/e**2; N=len(k)

def gls(X,yy,ww): A=(X.T*ww)@X; return np.linalg.solve(A,(X.T*ww)@yy), np.linalg.inv(A)
def slope_se(yy):
    X=np.column_stack([np.ones_like(k),k,mass]); b=gls(X,yy,w)[0][1]
    pat=(ra//18).astype(int); sl=[gls(X[pat!=p],yy[pat!=p],w[pat!=p])[0][1] for p in np.unique(pat)]
    sl=np.array(sl); se=np.sqrt((len(sl)-1)/len(sl)*np.sum((sl-sl.mean())**2)); return b,se

b0,se0=slope_se(y)
print("="*78); print("POWER TEST + amplitude reality-check — BOSS kappa_fg x Pantheon+ (648 SNe)"); print("="*78)
print(f"\n  measured slope b = {b0:+.2f} +/- {se0:.2f} (jackknife); N={N}")

# (1) is kappa_fg physical?
print(f"\n(1) AMPLITUDE REALITY-CHECK:")
print(f"    kappa_fg std = {k.std():.4f}; expected PHYSICAL foreground-convergence RMS ~ 0.010-0.020")
ratio=k.std()/0.015
print(f"    ratio (ours / physical) ~ {ratio:.2f}  ->  " +
      ("PHYSICAL-ish (amplitude meaningful)" if ratio>0.6 else
       "SUPPRESSED PROXY (~%.0f%% of physical): a galaxy-density proxy, NOT a true convergence map." % (ratio*100)))

# (2) power test: detection threshold + recovery of injected signals
print(f"\n(2) POWER TEST (inject b_inj*kappa_fg through the actual pipeline):")
print(f"    {'b_inj':>8s} {'recovered b':>12s} {'significance':>13s}")
for binj in (-15,-8,-4,-2.17,0,2.17,4,8,15):
    br,ser=slope_se(y+binj*k); print(f"    {binj:8.2f} {br:12.2f} {br/ser:+12.1f}σ")
b3=3*se0
print(f"\n    3-sigma detection threshold (proxy units): |b| > {b3:.1f}. 2-sigma: |b| > {2*se0:.1f}.")
print(f"    measured {b0:+.2f}+/-{se0:.2f} is consistent with ZERO ({abs(b0/se0):.1f}sigma).")

# (3) honest resolution
print(f"\n(3) RESOLUTION (honesty-guard, the referee-correct):")
if ratio<0.6:
    print(f"    kappa_fg is a SUPPRESSED galaxy-density PROXY (~{ratio*100:.0f}% of physical convergence), so the")
    print(f"    proxy->physical amplitude is UNCERTAIN -> the SLOPE cannot be cleanly compared to the physical DFD")
    print(f"    (+7..+15) or LCDM (-2.17) predictions. ROBUST statement = SIGN TEST: no preferred sign (+0.5sigma).")
    print(f"    => 'best available CALIBRATED GALAXY-DENSITY FOREGROUND PROXY', NOT 'fully calibrated convergence'.")
    print(f"    => inconclusive sign test; NO DFD signal; washes out prior hints; NOT a DFD win (and not a clean")
    print(f"       physical constraint either). Decisive amplitude test needs a real weak-lensing kappa map + more SNe.")
else:
    print(f"    kappa_fg is ~physical -> the null DOES constrain: a DFD slope above {b3:.1f} is excluded at 3sigma.")
