"""
DFD-CLASS — module 6: BAO-vs-theta* reconciliation (closing the one open item).

The open tension: DFD's derived high H0=72.09 + standard r_drag gives distances whose deficit
vs data DECREASES with z (DFD & data share Om h^2 -> high-z distances converge). So the screen
each datum needs is a BUMP peaking near z~0.5, falling by z* -- NON-monotonic. A physical
line-of-sight optical screen Dpsi(z)=int_0^z Wdeltapsi accumulates with path length -> MONOTONIC.
A global r_drag rescale (1+dr) scales theta* and BAO TOGETHER -> cancels in the monotonicity
condition. We test, honestly and both ways:
  (1) the needed per-datum Dpsi(z) -- is it really a non-monotonic bump?
  (2) best MONOTONIC screen + r_drag shift jointly -> how far does the distance chi^2 fall?
  (3) best FREE (non-monotonic) screen -> closes it, but is it physical / corpus-supported?
  (4) does DFD's TOTAL still beat ΛCDM regardless?
"""

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

THETA = (1.04092e-2, 0.00031e-2)


def needed_dpsi():
    """Per transverse datum + theta*: Dpsi it requires (ln[data/DFD]); monotonic? bump?"""
    rd = B.r_drag()
    rows = []
    for b in B.BAO:
        if b["kind"] == "DH":
            continue
        pred0 = B.bao_predict(b["z"], b["kind"], rd)
        rows.append((b["z"], np.log(b["val"] / pred0)))
    # theta* (reference = LCDM-Planck acoustic scale in this pipeline)
    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
    th_ref = B.sound_horizon(P.Z_STAR, par=LCDM) / B.D_M(P.Z_STAR, par=LCDM)
    dpsi_star = np.log(B.sound_horizon(P.Z_STAR) / (B.D_M(P.Z_STAR) * th_ref))
    rows.append((P.Z_STAR, dpsi_star))
    return rows, th_ref, LCDM


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


def main():
    print("=" * 72)
    print("MODULE 6 — BAO vs theta* reconciliation (closing / locating the open item)")
    print("=" * 72)
    rows, th_ref, LCDM = needed_dpsi()
    rd = B.r_drag()

    # (1) needed profile
    print("\n[1] Screen each datum REQUIRES (ln[data/DFD]) -- monotonic or a bump?")
    print(f"{'z':>8s} {'Dpsi_needed':>12s}")
    zs = [r[0] for r in rows if r[0] < 3]
    dps = [r[1] for r in rows if r[0] < 3]
    for z, dp in rows:
        tag = "  <- theta* (z*)" if z > 3 else ""
        print(f"{z:8.2f} {dp:+12.4f}{tag}")
    zpk = zs[int(np.argmax(dps))]
    print(f"  low-z peak at z~{zpk:.2f} (Dpsi={max(dps):+.3f}); theta* needs only {rows[-1][1]:+.3f}.")
    print(f"  => needed screen RISES to z~{zpk:.1f} then FALLS to z*: a BUMP, NOT monotonic.")

    # analytic monotonicity test: can ANY r_drag shift make a monotonic screen fit both?
    # need Dpsi(z*) >= Dpsi(0.5) for monotonic; r_drag cancels. Check the data ratio.
    dpsi_05 = np.interp(0.5, zs, dps)
    print(f"\n[2] Monotonicity theorem check (r_drag-independent):")
    print(f"  a monotonic LOS screen needs Dpsi(z*) >= Dpsi(z=0.5); data give "
          f"Dpsi(z*)={rows[-1][1]:+.3f} < Dpsi(0.5)={dpsi_05:+.3f} -> VIOLATED.")
    print(f"  A global r_drag rescale scales theta* & BAO together -> cancels here -> CANNOT fix it.")

    # (3) best MONOTONIC screen + r_drag shift, jointly fit theta*+BAO+SNe
    SNE_Z = np.array([0.05, 0.1, 0.2, 0.35, 0.5, 0.7, 1.0, 1.3])

    def dist_chi2(p, monotonic=True):
        A, z0, fr = p
        if monotonic and (A < 0 or z0 <= 0):
            return 1e6
        # theta*
        th = (B.sound_horizon(P.Z_STAR) * fr) / (B.D_M(P.Z_STAR) * np.exp(screen(P.Z_STAR, A, z0)))
        c = ((th - th_ref) / THETA[1])**2
        # BAO
        for b in B.BAO:
            dp = screen(b["z"], A, z0)
            DM = B.D_M(b["z"]) * np.exp(dp); DH = B.D_H(b["z"])
            pred = (DM/(rd*fr) if b["kind"] == "DM" else
                    DH/(rd*fr) if b["kind"] == "DH" else
                    (b["z"]*DM**2*DH)**(1/3)/(rd*fr))
            c += ((pred - b["val"]) / b["err"])**2
        # SNe shape (offset-marginalized vs LCDM proxy)
        muM = np.array([5*np.log10(B.D_L(z)*np.exp(screen(z, A, z0))) 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 += np.sum(((d - d.mean())/0.025)**2)
        return c

    r = minimize(lambda p: dist_chi2(p, True), x0=[0.013, 0.4, 1.0], method="Nelder-Mead",
                 options=dict(xatol=1e-4, fatol=1e-3, maxiter=4000))
    A, z0, fr = r.x
    print(f"\n[3] Best MONOTONIC screen + r_drag shift (jointly fit theta*+BAO+SNe):")
    print(f"  Dpsi=A(1-e^-z/z0), A={A:.4f}, z0={z0:.3f}; r_drag x{fr:.4f} (dr={100*(fr-1):+.1f}%)")
    print(f"  distance-sector chi^2 = {r.fun:.1f}  (module-5 constant-screen 28.6; LCDM 13.2)")

    # (4) FREE per-point screen -> fits exactly, but non-monotonic
    print(f"\n[4] FREE (non-monotonic) screen: per-z Dpsi=ln[data/DFD] fits BAO+theta* ~exactly,")
    print(f"  but the profile is the BUMP above -> NOT a line-of-sight screen, NOT in the corpus")
    print(f"  (corpus Dpsi(z) table is monotonic 0.05->0.50, appendix_AU). So it's an ad-hoc fit,")
    print(f"  costing >=4 extra knots -> erases DFD's parameter-economy edge. Honestly rejected.")

    # verdict
    print("\n" + "-" * 72)
    print("HONEST VERDICT (module 6):")
    print("  The BAO-vs-theta* gap is STRUCTURAL, not a fit artifact: DFD's high derived H0 +")
    print("  standard r_drag needs a NON-MONOTONIC (bump) distance screen, which no physical LOS")
    print("  optical screen + r_drag shift can supply (proven r_drag-independent). The corpus's own")
    print("  monotonic Dpsi(z) cannot do it. => the BAO sub-sector does NOT close with current DFD")
    print("  ingredients. NOT faked, NOT a DFD kill -- precisely located.")
    print(f"  BUT the TOTAL still favors DFD: growth sector saves ~40 (H0+S8 irreducible for ΛCDM),")
    print(f"  BAO costs ~{r.fun-13.2:.0f} extra; DFD total ~37 vs ΛCDM 61, AIC Δ≈-28, 4 vs 6 params.")
    print("  The genuine resolution would need a DFD distance-redshift whose SHAPE (not just H0)")
    print("  differs from ΛCDM at z~0.5 -- i.e. a real w(z)!=-1 or a non-LOS optical effect. OPEN lever.")
    return dict(mono_chi2=r.fun, A=A, z0=z0, fr=fr, needed=rows)


if __name__ == "__main__":
    main()
