"""
DFD-CLASS — module 10: the mu-fight + bookkeeping (the referee's brutal test).

Two questions, resolved separately and honestly:

(A) THE mu-FIGHT. The real screen integrand is NOT density x boost; it is
        s(delta) = (1+delta) * Q_mu(g/a0),   Q_mu(x) = (1+x)/x,   x = g/a0.
    More mass raises (1+delta) but ALSO raises g -> raises x -> LOWERS Q_mu. They fight.
    Halos: high delta but x>>1 so Q_mu->1 (boost dies). Voids: low delta but x<<1 so Q_mu~10
    (boost lives). The referee's claim: the screen is carried by void/filament sightlines, NOT host
    halos. We integrate s(delta) over the real lognormal density PDF with x(delta), galaxy-weighted,
    and report whether the NET is enhanced and FROM WHERE.

(B) THE BOOKKEEPING. Report three numbers so the comparison to Dpsi_need=0.022 is unambiguous:
        Dpsi_std   = screen with Q_mu=1 (ordinary lensing/magnification, ΛCDM-shared)
        Dpsi_total = screen with the full DFD Q_mu (coherent optical dilation, corpus eq:psitomo)
        Dpsi_excess= Dpsi_total - Dpsi_std  (the DFD-SPECIFIC part beyond standard lensing)
    Closure question (vs bare DFD): Dpsi_need is TOTAL relative to no-optical-effects bare DFD, and
    the corpus screen D_M^DFD=e^{Dpsi}D_bare is a COHERENT DILATION (shifts the BAO D_M, unlike pure
    magnification) -> the closure comparison is Dpsi_total vs 0.022. The EXCESS is the ΛCDM-discriminator.
"""

import numpy as np
from scipy.integrate import quad
import dfd_params as P

DPSI_LIN = 0.007
DPSI_NEED = 0.022
SIGMA = 0.9        # density-field rms at the screen smoothing scale (~few Mpc, z~0.5)
B_GAL = 1.0        # galaxy selection bias exponent (galaxies in overdensities)


def lognormal_pdf(opd, sigma=SIGMA):
    """PDF of (1+delta)=opd for a lognormal field with <delta>=0."""
    if opd <= 0:
        return 0.0
    g = np.log(opd) + 0.5*sigma**2
    return np.exp(-g**2/(2*sigma**2)) / (opd * sigma*np.sqrt(2*np.pi))


def x_of_opd(opd, x_norm, p=0.5):
    """local g/a0 vs overdensity. Acceleration scales sub-linearly with enclosed density: x ~ (1+d)^p.
    x_norm sets the cosmic-mean acceleration environment of the screen LOS."""
    return x_norm * opd**p


def Qmu(x):
    return (1+x)/x


def moments(x_norm, with_Q=True, gal=False):
    """<s> over the PDF; s=(1+d)*Q_mu (or (1+d) if with_Q=False). gal=True weights by (1+d)^B_GAL."""
    def w(opd):
        base = lognormal_pdf(opd)
        if gal:
            base *= opd**B_GAL
        return base
    def num(opd):
        s = opd * (Qmu(x_of_opd(opd, x_norm)) if with_Q else 1.0)
        return w(opd) * s
    N = quad(num, 1e-3, 50, limit=120)[0]
    D = quad(w, 1e-3, 50, limit=120)[0]
    return N/D


def env_decomposition(x_norm):
    """Contribution to the galaxy-weighted screen from voids / filaments / halos."""
    bins = [("void", 1e-3, 0.5), ("filament", 0.5, 3.0), ("halo", 3.0, 50.0)]
    out = {}
    tot = 0.0
    for name, lo, hi in bins:
        def num(opd):
            return lognormal_pdf(opd)*opd**B_GAL * opd*Qmu(x_of_opd(opd, x_norm))
        c = quad(num, lo, hi, limit=80)[0]
        out[name] = c; tot += c
    for k in out:
        out[k] /= tot
    return out


def main():
    print("=" * 74)
    print("MODULE 10 — the mu-fight + bookkeeping (the referee's brutal test)")
    print("=" * 74)

    # (A) the mu-fight: scan x_norm, find what galaxy-weighted enhancement it gives
    print("\n[A] mu-FIGHT: galaxy-weighted screen enhancement B_NL vs cosmic-mean x_norm of the LOS:")
    print(f"{'x_norm':>8s} {'B_NL':>7s} {'Dpsi_tot':>9s} {'void%':>6s} {'fil%':>6s} {'halo%':>6s}")
    target_xn = None
    for xn in [0.1, 0.2, 0.3, 0.5, 1.0, 2.0, 5.0]:
        s_gal = moments(xn, with_Q=True, gal=True)
        s_std_lin = moments(xn, with_Q=False, gal=True)   # standard, no Q
        # B_NL relative to the LINEAR (no-Q, no-nonlinear) baseline that gave 0.007
        B = s_gal / s_std_lin * Qmu(xn)    # the Q enhancement over the std galaxy screen
        dpsi_tot = DPSI_LIN * B
        env = env_decomposition(xn)
        print(f"{xn:8.2f} {B:7.2f} {dpsi_tot:9.4f} {100*env['void']:6.0f} "
              f"{100*env['filament']:6.0f} {100*env['halo']:6.0f}")
    # closure x_norm = where dpsi_tot is CLOSEST to the need (steep function -> unique crossing near x~1)
    xn_grid = np.linspace(0.3, 3.0, 80)
    dps = [DPSI_LIN * (moments(x, True, True)/moments(x, False, True)) * Qmu(x) for x in xn_grid]
    target_xn = xn_grid[int(np.argmin(np.abs(np.array(dps) - DPSI_NEED)))]
    print(f"\n  -> Dpsi_total = need ({DPSI_NEED}) at x_norm ~ {target_xn:.2f}  (NOT 0.2-0.3 -- that overshoots).")
    print(f"  -> EXPONENTIALLY SENSITIVE to x_norm: 0.54 at x=0.1 down to 0.009 at x=5. The closure")
    print(f"     magnitude is NOT robustly predicted; it hinges on the real LOS acceleration environment.")

    # (B) the bookkeeping at the x_norm that closes (or best)
    xn = target_xn
    s_gal_Q = moments(xn, with_Q=True, gal=True)
    s_gal_std = moments(xn, with_Q=False, gal=True)
    dpsi_std = DPSI_LIN                                  # ordinary lensing/mag (ΛCDM-shared) ~0.006-0.007
    dpsi_total = DPSI_LIN * (s_gal_Q/s_gal_std) * Qmu(xn)
    dpsi_excess = dpsi_total - dpsi_std
    print(f"\n[B] BOOKKEEPING (at x_norm={xn}):")
    print(f"   Dpsi_std    (Q_mu=1, ΛCDM-shared lensing)      = {dpsi_std:.4f}")
    print(f"   Dpsi_total  (full DFD optical dilation)        = {dpsi_total:.4f}")
    print(f"   Dpsi_excess (DFD-specific, beyond std lensing) = {dpsi_excess:.4f}")
    print(f"   need (vs BARE DFD) = {DPSI_NEED}")

    print("\n[C] sigma8 consistency: the sigma8 scale (8 Mpc/h, denser) sits at x~5.85 -> Q_mu=1.17.")
    print(f"   The screen needs x_norm~{xn} (void-dominated LOS) << 5.85 -> SCALE-SEPARATED. One mu-law")
    print(f"   gives a big screen (x~{xn}) AND modest growth (x~5.85). The referee's x_BAO~0.2-0.3 target.")

    print("\n" + "-" * 74)
    print("HONEST VERDICT (module 10):")
    env = env_decomposition(xn)
    carrier = max(env, key=env.get)
    print(f"  mu-FIGHT outcome: NET ENHANCED, carried by the {carrier.upper()} ({100*env[carrier]:.0f}% of the")
    print(f"  screen) -- {'voids/filaments as the referee predicted (host halos contribute '+str(round(100*env['halo']))+'%, mu-boost dies there)' if carrier!='halo' else 'HALOS (unexpected -- check)'}.")
    print(f"  CLOSURE (Dpsi_total vs need): {dpsi_total:.3f} vs {DPSI_NEED} at x_norm~{xn:.2f} -> "
          f"{'CLOSES' if dpsi_total>=0.9*DPSI_NEED else 'short'}, IF the screen is a coherent dilation")
    print(f"  (corpus eq:psitomo D_M=e^Dpsi D_bare, shifts the BAO scale) -- which is the corpus claim.")
    print(f"  DISCRIMINATOR (excess vs need): {dpsi_excess:.3f} vs {DPSI_NEED} -> if 'need' means EXTRA")
    print(f"  beyond ΛCDM lensing, still {'closes' if dpsi_excess>=0.9*DPSI_NEED else 'SHORT by '+str(round((DPSI_NEED-dpsi_excess)*1000))+'e-3'} -- the bookkeeping decides.")
    print(f"  ⇒ DECISIVE TEST (real maps): is x_bar_BAO/SN ~ 0.2-0.3 (closes) or ~1 (fails) on real")
    print(f"  galaxy-selected void/filament sightlines? + does the BAO/SN residual carry the DFD mu-law")
    print(f"  scale-dependent correlation with foreground density AFTER standard lensing is removed?")
    return dict(target_xn=target_xn, dpsi_std=dpsi_std, dpsi_total=dpsi_total,
                dpsi_excess=dpsi_excess, carrier=carrier, env=env)


if __name__ == "__main__":
    main()
