"""
DFD-CLASS — module 7a: does an evolving w(z) close the BAO gap? (feasibility, NOT a DFD fit)

Module 6 proved: BAO closes only if DFD's distance SHAPE (not just H0) differs from LCDM at
z~0.5 -> a derived w(z)!=-1. Here we test the LEVER's viability: with DFD's DERIVED H0=72.09 and
Om=0.275 held FIXED, and r_drag STANDARD (DE negligible at z>1000, so r_drag unchanged), scan CPL
w(a)=w0+wa(1-a) and ask:
  - is there a (w0,wa) that fits CMB theta* + BAO + SNe TOGETHER (closes the gap)?
  - is it DESI-2024-aligned (w0>-1, wa<0)?
This is a DIAGNOSTIC of whether the w(z) lever exists. It does NOT claim DFD predicts this w(z) --
that is decided by the corpus (module-7 analysis). Derived-only rule respected: the answer here
just tells us the TARGET; the corpus tells us DFD's ACTUAL w(z).
"""

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

THETA = (1.04092e-2, 0.00031e-2)
SNE_Z = np.array([0.05, 0.1, 0.2, 0.35, 0.5, 0.7, 1.0, 1.3])


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 rho_de(a, w0, wa):
    return a**(-3*(1 + w0 + wa)) * np.exp(-3 * wa * (1 - a))


def E_cpl(z, w0, wa):
    a = 1.0 / (1 + z)
    return np.sqrt(P.OMEGA_M * (1 + z)**3 + P.OMEGA_R * (1 + z)**4
                   + P.OMEGA_LAMBDA * rho_de(a, w0, wa))


def D_M_cpl(z, w0, wa):
    integ, _ = quad(lambda zz: 1.0 / E_cpl(zz, w0, wa), 0.0, z, limit=200)
    return (P.C_KMS / P.H0) * integ


def D_H_cpl(z, w0, wa):
    return P.C_KMS / (P.H0 * E_cpl(z, w0, wa))


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


TH_REF = theta_ref()
RD = B.r_drag()           # standard sound horizon, unchanged by late-time w(z)


def dist_chi2(w0, wa, parts=False):
    # theta*: r_s(z*) standard / D_M(z*) with CPL late-time expansion
    # (sound horizon integral itself is z>1000 where DE is negligible -> use standard RD-scaled DA)
    DA_star = D_M_cpl(P.Z_STAR, w0, wa)
    th = B.sound_horizon(P.Z_STAR) / DA_star
    c_th = ((th - TH_REF) / THETA[1])**2
    c_bao = 0.0
    for b in B.BAO:
        DM = D_M_cpl(b["z"], w0, wa); DH = D_H_cpl(b["z"], w0, wa)
        pred = (DM/RD if b["kind"] == "DM" else DH/RD if b["kind"] == "DH"
                else (b["z"]*DM**2*DH)**(1/3)/RD)
        c_bao += ((pred - b["val"]) / b["err"])**2
    muM = np.array([5*np.log10(D_M_cpl(z, w0, wa)*(1+z)) for z in SNE_Z])
    muD = np.array([5*np.log10(B.D_L(z, LCDM)) for z in SNE_Z])
    d = muM - muD
    c_sne = np.sum(((d - d.mean()) / 0.025)**2)
    if parts:
        return c_th, c_bao, c_sne
    return c_th + c_bao + c_sne


def main():
    print("=" * 72)
    print("MODULE 7a — does an evolving w(z) close BAO? (feasibility, fixed derived H0/Om)")
    print("=" * 72)

    # baseline w=-1 (the current DFD)
    c0 = dist_chi2(-1.0, 0.0, parts=True)
    print(f"\n[w=-1, current DFD]  theta*={c0[0]:.1f}  BAO={c0[1]:.1f}  SNe={c0[2]:.1f}  "
          f"total={sum(c0):.1f}  (LCDM distance-sector ~13.2)")

    # minimize over (w0, wa)
    res = minimize(lambda p: dist_chi2(*p), x0=[-0.9, -0.3], method="Nelder-Mead",
                   options=dict(xatol=1e-3, fatol=1e-2, maxiter=2000))
    w0, wa = res.x
    cb = dist_chi2(w0, wa, parts=True)
    print(f"\n[best-fit CPL w(z)]  w0={w0:+.3f}  wa={wa:+.3f}")
    print(f"   theta*={cb[0]:.1f}  BAO={cb[1]:.1f}  SNe={cb[2]:.1f}  total={sum(cb):.1f}")
    print(f"   w(z=0)={w0:+.3f}, w(z=0.5)={w0+wa*(1-1/1.5):+.3f}, w(early)={w0+wa:+.3f}")

    # DESI 2024 DR1 best-fit (w0wa CDM, DESI+CMB+PantheonPlus): w0=-0.827, wa=-0.75
    DESI = (-0.827, -0.75)
    cd = dist_chi2(*DESI, parts=True)
    print(f"\n[DESI-2024 w0wa = ({DESI[0]}, {DESI[1]})]  theta*={cd[0]:.1f}  BAO={cd[1]:.1f}  "
          f"SNe={cd[2]:.1f}  total={sum(cd):.1f}")

    print("\n" + "-" * 72)
    print("HONEST VERDICT (module 7a — is the w(z) lever viable?):")
    closes = sum(cb) < 16
    near_desi = abs(w0 - DESI[0]) < 0.25 and (wa < 0)
    print(f"  w=-1 distance chi^2 {sum(c0):.0f} -> best evolving-w {sum(cb):.0f} "
          f"(LCDM ~13). The lever {'CLOSES the gap' if closes else 'helps but does not fully close'}.")
    print(f"  Required: w0={w0:+.2f}, wa={wa:+.2f} (quintessence-like, w>-1 today). DESI-2024 prefers "
          f"({DESI[0]},{DESI[1]}) -> {'SAME DIRECTION/aligned' if near_desi else 'different'}.")
    print(f"  => IF DFD's psi-DE derives a w(z) of THIS form, BAO closes AND it matches DESI = full win.")
    print(f"     The corpus (module-7 analysis) decides whether DFD's ACTUAL w(z) is this or exactly -1.")
    print(f"     NOT a DFD fit -- this is the TARGET the derivation must hit. Honesty-guard respected.")
    return dict(w0=w0, wa=wa, chi2_w1=sum(c0), chi2_best=sum(cb), chi2_desi=sum(cd), desi=DESI)


if __name__ == "__main__":
    main()
