"""
DFD-CLASS — Joint likelihood / scoreboard (module 5).

Sums every dataset into ONE chi^2 for DFD vs LambdaCDM, with explicit parameter counting.
Two sectors:
  GROWTH sector (screen-INDEPENDENT, DFD wins on derived params):
     H0 (SH0ES), S8 (DES+KiDS weak lensing), fsigma8 (BOSS/eBOSS RSD).
  DISTANCE sector (one shared optical screen Dpsi(z) must satisfy ALL):
     CMB acoustic scale theta*, BAO (12 pts), SNe (Pantheon shape).

The decisive question: can ONE physical screen Dpsi(z) (monotonic, ->0 at z=0) fit CMB theta*
AND BAO AND SNe together, and does DFD's growth-sector win outweigh its distance-sector cost in
the TOTAL chi^2 -- at FEWER free parameters than LCDM's 6? Honest, honesty-guard both ways.
"""

import numpy as np
from scipy.optimize import minimize
import dfd_params as P
import background as B

# ---------- growth-sector results from module 4 (CAMB), screen-independent ----------
GROWTH = dict(
    dfd=dict(H0=72.09, S8=0.7845, fs8_chi2=6.57),   # perturbations.py outputs
    lcdm=dict(H0=67.36, S8=0.8316, fs8_chi2=7.30),
)
SH0ES = (73.04, 1.04)
S8_WL = (0.766, 0.020)
THETA = (1.04092e-2, 0.00031e-2)     # Planck 100*theta_MC/100
N_FS8 = 8

# ---------- LambdaCDM background (Planck) for distance-sector proxy ----------
class LCDM:
    H0 = P.PLANCK["H0"]; OMEGA_M = P.PLANCK["Omega_m"]
    OMEGA_R = P.OMEGA_R_H2 / P.PLANCK["h"]**2
    OMEGA_LAMBDA = 1 - P.PLANCK["Omega_m"] - OMEGA_R
    OMEGA_B_H2 = P.OMEGA_B_H2; OMEGA_G_H2 = P.OMEGA_G_H2


def screen(z, A, z0):
    return A * (1.0 - np.exp(-z / max(z0, 1e-3)))


# ---------- distance-sector chi^2 (DFD with screen) ----------
# theta* reference = the LCDM-Planck value IN THIS SAME analytic pipeline (LCDM fits Planck's
# acoustic scale by construction), so the absolute r_s systematic cancels and the comparison is
# fair: LCDM theta* chi^2 = 0, DFD must use its screen to reach the same acoustic scale.
def theta_chi2(A, z0, ref):
    rs = B.sound_horizon(P.Z_STAR)
    DA = B.D_M(P.Z_STAR) * np.exp(screen(P.Z_STAR, A, z0))
    th = rs / DA
    return ((th - ref) / THETA[1])**2, th


def bao_chi2(A, z0, par=None, rd=None):
    par = par or P
    rd = rd if rd is not None else (B.r_drag() if par is P else P.PLANCK["r_drag"])
    chi2 = 0.0
    for b in B.BAO:
        dp = screen(b["z"], A, z0) if par is P else 0.0
        DM = B.D_M(b["z"], par) * np.exp(dp)
        DH = B.D_H(b["z"], par)
        if b["kind"] == "DM":
            pred = DM / rd
        elif b["kind"] == "DH":
            pred = DH / rd
        else:
            pred = (b["z"] * DM**2 * DH)**(1/3) / rd
        chi2 += ((pred - b["val"]) / b["err"])**2
    return chi2


SNE_Z = np.array([0.05, 0.1, 0.2, 0.35, 0.5, 0.7, 1.0, 1.3])
SNE_ERR = 0.025   # Pantheon+ binned ~mag


def sne_chi2(A, z0, par=None):
    """Offset-marginalized SNe shape chi^2 vs LCDM-Planck distance proxy (Pantheon-quality)."""
    par = par or P
    mu_model = np.array([5*np.log10(B.D_L(z, par) * (np.exp(screen(z, A, z0)) if par is P else 1.0))
                         for z in SNE_Z])
    mu_data = np.array([5*np.log10(B.D_L(z, LCDM)) for z in SNE_Z])  # proxy (LCDM fits Pantheon ~1sig)
    d = mu_model - mu_data
    off = np.sum(d) / len(d)            # analytic marginalization over absolute offset (M_B/H0)
    return np.sum(((d - off) / SNE_ERR)**2)


def distance_chi2(params, ref):
    A, z0 = params
    c_th, _ = theta_chi2(A, z0, ref)
    return c_th + bao_chi2(A, z0) + sne_chi2(A, z0)


def main():
    print("=" * 74)
    print("DFD-CLASS JOINT LIKELIHOOD — full DFD vs LCDM scoreboard")
    print("=" * 74)

    # ---- growth sector (no screen) ----
    g = GROWTH
    h0_d = ((g["dfd"]["H0"] - SH0ES[0]) / SH0ES[1])**2
    h0_l = ((g["lcdm"]["H0"] - SH0ES[0]) / SH0ES[1])**2
    s8_d = ((g["dfd"]["S8"] - S8_WL[0]) / S8_WL[1])**2
    s8_l = ((g["lcdm"]["S8"] - S8_WL[0]) / S8_WL[1])**2
    fs8_d = g["dfd"]["fs8_chi2"]; fs8_l = g["lcdm"]["fs8_chi2"]

    # ---- distance sector: fit DFD screen (A,z0); LCDM defines the theta* reference ----
    TH_REF = LCDM_theta()                          # LCDM fits Planck's acoustic scale
    res = minimize(distance_chi2, x0=[0.013, 0.4], args=(TH_REF,), method="Nelder-Mead",
                   options=dict(xatol=1e-4, fatol=1e-3))
    A, z0 = res.x
    th_d, th_val = theta_chi2(A, z0, TH_REF)
    bao_d = bao_chi2(A, z0)
    sne_d = sne_chi2(A, z0)
    # LCDM distance sector (theta* = 0 by construction; it IS the reference)
    th_l = 0.0
    bao_l = bao_chi2(0, 1, par=LCDM, rd=P.PLANCK["r_drag"])
    sne_l = sne_chi2(0, 1, par=LCDM)

    print(f"\nDFD optical screen fit: Dpsi(z) = {A:.4f}*(1-e^(-z/{z0:.3f}))  "
          f"[Dpsi(z*)={screen(P.Z_STAR,A,z0):.4f}, Dpsi(0.5)={screen(0.5,A,z0):.4f}]")

    rows = [
        ("H0 (SH0ES)",        1, h0_d, h0_l),
        ("S8 (DES+KiDS WL)",  1, s8_d, s8_l),
        ("fsigma8 (RSD)",     N_FS8, fs8_d, fs8_l),
        ("CMB theta* (acoustic)", 1, th_d, th_l),
        ("BAO (12 pts)",      12, bao_d, bao_l),
        ("SNe (Pantheon shape)", len(SNE_Z), sne_d, sne_l),
    ]
    print(f"\n{'dataset':24s} {'N':>3s} {'chi2_DFD':>9s} {'chi2_LCDM':>10s}   winner")
    tot_d = tot_l = tot_n = 0
    for name, N, cd, cl in rows:
        tot_d += cd; tot_l += cl; tot_n += N
        w = "DFD" if cd < cl - 0.5 else ("LCDM" if cl < cd - 0.5 else "tie")
        print(f"{name:24s} {N:3d} {cd:9.2f} {cl:10.2f}   {w}")
    print("-" * 60)
    print(f"{'TOTAL':24s} {tot_n:3d} {tot_d:9.2f} {tot_l:10.2f}")

    # ---- parameter counting ----
    # LCDM: {ombh2, omch2, H0, As, ns, tau} = 6
    # DFD : derived {H0, omch2(theta_i fixed), As, ns} + free {screen A, z0, tau} ... count free:
    p_l = 6
    p_dfd = 2 + 1 + 1   # screen(A,z0) + theta_i(misalignment angle) + tau ; H0/ns/As/Ombh2 derived/BBN
    print(f"\nFree parameters:  DFD ~ {p_dfd}  (screen A,z0 + theta_i + tau; H0/n_s/A_s derived)")
    print(f"                  LCDM = {p_l}  (ombh2,omch2,H0,A_s,n_s,tau all fitted)")
    aic_d = tot_d + 2 * p_dfd
    aic_l = tot_l + 2 * p_l
    print(f"chi^2/dof:  DFD {tot_d/(tot_n-p_dfd):.2f}   LCDM {tot_l/(tot_n-p_l):.2f}")
    print(f"AIC (chi2+2k):  DFD {aic_d:.1f}   LCDM {aic_l:.1f}  -> "
          f"{'DFD preferred' if aic_d < aic_l else 'LCDM preferred'} (delta={aic_d-aic_l:+.1f})")

    print("\n" + "-" * 74)
    print("HONEST VERDICT (module 5 — the full scoreboard):")
    print(f"  GROWTH sector (H0+S8+fsigma8): DFD {h0_d+s8_d+fs8_d:.1f} vs LCDM {h0_l+s8_l+fs8_l:.1f}"
          f"  -> DFD WINS by {h0_l+s8_l+fs8_l-(h0_d+s8_d+fs8_d):.0f} (the combined H0/S8/growth tension).")
    print(f"  DISTANCE sector (theta*+BAO+SNe): DFD {th_d+bao_d+sne_d:.1f} vs LCDM {th_l+bao_l+sne_l:.1f}"
          f"  -> LCDM wins by {th_d+bao_d+sne_d-(th_l+bao_l+sne_l):.0f} (the screen can't fully")
    print(f"     reconcile theta* [{screen(P.Z_STAR,A,z0):.3f}] vs BAO [~0.034]: residual BAO chi2={bao_d:.0f}).")
    if tot_d < tot_l:
        print(f"  TOTAL: DFD {tot_d:.0f} < LCDM {tot_l:.0f} at FEWER params ({p_dfd}<{p_l}) -> DFD ahead overall.")
    else:
        print(f"  TOTAL: DFD {tot_d:.0f} vs LCDM {tot_l:.0f}. DFD's growth win {'covers' if tot_d<tot_l else 'does NOT fully cover'}"
              f" its distance cost; at {p_dfd} vs {p_l} params, AIC {'favors DFD' if aic_d<aic_l else 'favors LCDM'}.")
    print(f"  The ONE unresolved item is the acoustic-scale-vs-BAO screen amplitude (high-H0/std-r_drag).")
    return dict(tot_dfd=tot_d, tot_lcdm=tot_l, aic_dfd=aic_d, aic_lcdm=aic_l, A=A, z0=z0)


def LCDM_theta():
    rs = B.sound_horizon(P.Z_STAR, par=LCDM)
    DA = B.D_M(P.Z_STAR, par=LCDM)
    return rs / DA


if __name__ == "__main__":
    main()
