"""
DFD SN-LENSING SIGN DISCRIMINATOR -- forecast + figure (theorem-grade, data-gated).

The test: regress standardized SN Ia Hubble residuals delta_mu on the FOREGROUND weak-lensing
convergence kappa_fg along each line of sight.
  - DFD (optical screen, n=e^psi): a SN seen THROUGH foreground mass looks FARTHER -> POSITIVE slope.
        s_DFD = (5/ln10) * Q * g_geo,  Q = 1/mu(xbar),  mu(x)=x/(1+x),  a0 = 2 sqrt(alpha) c H(z).
  - LCDM (GR magnification): the same SN looks BRIGHTER/CLOSER -> NEGATIVE slope.
        s_LCDM = -5/ln10 = -2.17  (exact, convention-fixed; = d/dk[-2.5 log10(1+2k)] at k->0).
The SIGN is the discriminator (theorem given the DFD optical axiom); the magnitude band only sets power.

Power: sigma_s = sigma_int / (sqrt(N) * sigma_kappa_eff);  Z = (s_DFD - s_LCDM)/sigma_s.
Anchored so 648 SNe + a TRUE convergence map = 6 sigma => pins sigma_kappa_eff(true) = 0.0025.

HONEST STATUS: existing data ~1-1.5 sigma (map-noise limited, NOT theory limited). This is a
pre-registered prediction; the decisive >=3 sigma test awaits LSST/Roman convergence maps. NOT a detection.

Outputs: console table + fig_snlensing_discriminator.png (2 panels). numpy + matplotlib only.
"""
import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

# ----------------------------------------------------------------------------- constants
ALPHA   = 1.0 / 137.035999          # fine-structure constant
LN10    = np.log(10.0)
PREFAC  = 5.0 / LN10                 # = 2.17147, the magnitude-system prefactor (exact)
SIGMA_INT = 0.13                     # SN Ia intrinsic scatter (mag)

def mu(x):
    return x / (1.0 + x)

# ----------------------------------------------------------------------------- DFD slope
def s_dfd(xbar=0.24, g_geo=1.0):
    """DFD slope d(delta_mu)/d(kappa): (5/ln10) * Q * g_geo, Q = 1/mu(xbar). Strictly > 0."""
    Q = 1.0 / mu(xbar)
    return PREFAC * Q * g_geo, Q

# ----------------------------------------------------------------------------- LCDM slope
def s_lcdm():
    """GR magnification slope: -2.5*log10(1+2k) -> d/dk at 0 = -5/ln10. Exact, robustly negative."""
    return -PREFAC

# ----------------------------------------------------------------------------- forecast
def sigma_slope(N, sigma_kappa_eff, sigma_int=SIGMA_INT):
    return sigma_int / (np.sqrt(N) * sigma_kappa_eff)

def Zsep(N, sigma_kappa_eff, sep):
    return sep / sigma_slope(N, sigma_kappa_eff)

def N_for_sigma(sep, sigma_kappa_eff, target=3.0, sigma_int=SIGMA_INT):
    return (target * sigma_int / (sep * sigma_kappa_eff)) ** 2

# ----------------------------------------------------------------------------- main
def main():
    sD_cen, Qcen = s_dfd(0.24, 1.0)      # kernel-matched central
    sD_lo, _     = s_dfd(0.45, 1.0)      # weak boost
    sD_hi, _     = s_dfd(0.18, 2.0)      # strong boost + full gamma=1 doubling
    sL           = s_lcdm()
    bD           = 10.0                  # conservative slope hard-coded in preflight_forecast.py
    sep          = bD - sL               # physical DFD-vs-LCDM separation = 12.17

    print("=" * 78)
    print("DFD SN-LENSING SIGN DISCRIMINATOR -- forecast (data-gated, pre-registered)")
    print("=" * 78)
    print(f"  5/ln10 prefactor (exact)        = {PREFAC:.5f}")
    print(f"  Q = 1/mu(xbar=0.24)             = {Qcen:.3f}")
    print(f"  s_DFD central (kernel-matched)  = +{sD_cen:.2f} mag/kappa  [band +{sD_lo:.1f} .. +{sD_hi:.1f}]")
    print(f"  s_DFD conservative (preflight)  = +{bD:.2f} mag/kappa  (bD hard-coded)")
    print(f"  s_LCDM (GR magnification)       = {sL:+.2f} mag/kappa  (= -5/ln10, exact, NEGATIVE)")
    print(f"  physical separation s_DFD-s_LCDM= {sep:.2f} mag/kappa  (definite SIGN FLIP)")

    # anchor: 648 SNe + true map = 6 sigma -> solve for sigma_kappa_eff(true)
    N_anchor, Z_anchor = 648, 6.0
    skeff_true = Z_anchor * SIGMA_INT / (sep * np.sqrt(N_anchor))
    print(f"\n  ANCHOR: {N_anchor} SNe + TRUE map = {Z_anchor:.0f}sigma  ->  sigma_kappa_eff(true) = {skeff_true:.4f}")

    rows = [
        ("DES-era proxy (current)",        1378, 0.0002),
        ("DES-era proxy (best case)",      4000, 0.0002),
        ("DES-era, IF true map now",       1378, 0.0025),
        ("LSST Y1 (good wl kappa)",       30000, 0.0013),
        ("LSST Y1 (deep kappa)",          30000, 0.0019),
        ("Roman (deep, low-noise kappa)", 20000, 0.0025),
    ]
    print(f"\n  {'survey (kappa map)':32s} {'N':>7s} {'sk_eff':>8s} {'sigma_s':>8s} {'Z':>7s} {'N_3sig':>8s}")
    print("  " + "-" * 74)
    table = []
    for lab, N, sk in rows:
        ss = sigma_slope(N, sk)
        Z  = Zsep(N, sk, sep)
        N3 = N_for_sigma(sep, sk)
        table.append((lab, N, sk, ss, Z, N3))
        print(f"  {lab:32s} {N:7d} {sk:8.4f} {ss:8.2f} {Z:6.2f}s {N3:8.0f}")

    print("\n  READ-OUT: bottleneck is the MAP (Z ~ sigma_kappa_eff, linear), not the SNe (sqrt(N)).")
    print("  With the noisy proxy ~22,000 SNe for 3sigma; with a true map only ~160. LSST/Roman clear 3sigma.")
    print("\n  HONEST STATUS: existing data ~1-1.5 sigma (all DFD-positive sign, map-noise limited, NOT")
    print("  theory limited). This is a PRE-REGISTERED prediction; decisive >=3sigma awaits LSST/Roman.")
    print("  NOT a detection, NOT a falsification.")

    make_figure(sD_cen, sD_lo, sD_hi, bD, sL, table, skeff_true, sep)
    print("\n  Figure written: fig_snlensing_discriminator.png")
    return dict(s_DFD_central=sD_cen, s_DFD_band=(sD_lo, sD_hi), bD=bD, s_LCDM=sL,
                separation=sep, skeff_true=skeff_true, table=table)

# ----------------------------------------------------------------------------- figure
def make_figure(sD_cen, sD_lo, sD_hi, bD, sL, table, skeff_true, sep):
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 5.2))

    # --- Panel 1: opposite-sign slope cartoon -------------------------------
    k = np.linspace(0.0, 0.05, 200)
    dfd_line  = bD * k                          # conservative DFD slope (positive)
    dfd_lo    = sD_lo * k
    dfd_hi    = sD_hi * k
    lcdm_line = -2.5 * np.log10(1.0 + 2.0 * k)  # exact GR magnification residual (negative)

    ax1.fill_between(k, dfd_lo, dfd_hi, color="#2a7fff", alpha=0.18,
                     label=r"DFD band  $s\in[+7,+28]$")
    ax1.plot(k, dfd_line, color="#1457d6", lw=2.6,
             label=r"DFD central  $s\simeq+11$  (looks FARTHER)")
    ax1.plot(k, lcdm_line, color="#d62728", lw=2.6,
             label=r"$\Lambda$CDM  $s=-2.17$  (looks BRIGHTER)")
    ax1.axhline(0, color="0.55", lw=0.9)
    ax1.set_xlabel(r"foreground convergence  $\kappa_{\rm fg}$")
    ax1.set_ylabel(r"SN Hubble residual  $\delta\mu$  (mag)")
    ax1.set_title("Opposite-sign slope: the discriminator")
    ax1.text(0.026, bD * 0.045 * 0.62, "FARTHER\n(DFD)", color="#1457d6",
             fontsize=10, ha="center", va="center", fontweight="bold")
    ax1.text(0.030, lcdm_line[-1] * 0.62, "CLOSER\n($\\Lambda$CDM)", color="#d62728",
             fontsize=10, ha="center", va="center", fontweight="bold")
    ax1.legend(loc="upper left", fontsize=8.5, framealpha=0.9)
    ax1.grid(alpha=0.25)

    # --- Panel 2: achieved sigma vs N, per map quality ----------------------
    Ngrid = np.logspace(2, 4.6, 300)
    surveys = [
        ("proxy (current, sk=0.0002)", 0.0002, "#9467bd", "-"),
        ("DES true map (sk=0.0025)",   0.0025, "#8c564b", "--"),
        ("LSST Y1 (sk=0.0013)",        0.0013, "#2ca02c", "-"),
        ("Roman deep (sk=0.0025)",     0.0025, "#ff7f0e", "-."),
    ]
    seen = set()
    for lab, sk, col, ls in surveys:
        if sk in seen and "Roman" in lab:   # Roman shares sk with DES true; offset visually only in label
            pass
        seen.add(sk)
        Z = Zsep(Ngrid, sk, sep)
        ax2.plot(Ngrid, Z, color=col, ls=ls, lw=2.2, label=lab)
    ax2.axhline(3.0, color="k", lw=1.4, ls=":")
    ax2.text(1.1e2, 3.25, r"$3\sigma$ threshold", fontsize=9)
    # markers for the real anchors
    ax2.scatter([1378], [Zsep(1378, 0.0002, sep)], color="#9467bd", zorder=5, s=45,
                edgecolor="k", label="current proxy (~0.74$\\sigma$)")
    ax2.scatter([30000], [Zsep(30000, 0.0013, sep)], color="#2ca02c", zorder=5, s=45, edgecolor="k")
    ax2.scatter([20000], [Zsep(20000, 0.0025, sep)], color="#ff7f0e", zorder=5, s=45, edgecolor="k")
    ax2.set_xscale("log")
    ax2.set_yscale("log")
    ax2.set_xlabel(r"number of SNe  $N$")
    ax2.set_ylabel(r"achieved separation  $Z$  ($\sigma$)")
    ax2.set_title("Power: the MAP is the bottleneck, not the SNe")
    ax2.legend(loc="upper left", fontsize=8.0, framealpha=0.9)
    ax2.grid(alpha=0.25, which="both")

    fig.suptitle("DFD SN-lensing sign discriminator  (pre-registered, data-gated; existing data ~1-1.5$\\sigma$)",
                 fontsize=12.5, y=1.01)
    fig.tight_layout()
    fig.savefig("fig_snlensing_discriminator.png", dpi=150, bbox_inches="tight")
    plt.close(fig)

if __name__ == "__main__":
    main()
