from dataclasses import dataclass, field import math from typing import Optional from dataclasses import dataclass import math @dataclass class Parameter: label: str value: float min: float max: float step: float @dataclass class Toggle: label: str value: bool @dataclass class SimulationParams: # --- Cosmology --- H0: Parameter = field(default_factory=lambda: Parameter("H0 (base)", 67.4, 63.0, 70.0, 0.1 )) Om: Parameter = field(default_factory=lambda: Parameter("Omega_m", 0.315, 0.28, 0.38, 0.001 )) # --- Crust enhancement --- dsw: Toggle = field(default_factory=lambda: Toggle("DSW Profile C", True )) B: Parameter = field(default_factory=lambda: Parameter("B (amplitude)", 0.25, 0.10, 0.40, 0.005 )) gamma: Parameter = field(default_factory=lambda: Parameter("gamma (DE weight)",3.64, 0.0, 5.0, 0.01 )) z_s: Parameter = field(default_factory=lambda: Parameter("z_s (scale)", 0.65, 0.1, 2.5, 0.01 )) z_harm: Parameter = field(default_factory=lambda: Parameter("z_harm", -0.25, -0.50, 0.10, 0.01 )) alpha: Parameter = field(default_factory=lambda: Parameter("alpha (tail)", 0.50, 0.30, 1.20, 0.01 )) # --- G_eff suppression: TWO Weinberg-anchored channels --------------- # The canonical S8 = 0.7923 comes from BOTH G_eff dips below (NOT a single # eta=0.181 channel — that alone gives ~0.826, barely below LCDM). Both # efficiencies are fixed by the Weinberg angle (alpha_mf = 0.3008): # * downstream channel (here): eta = 2*alpha^2*(1-sin2thW) = 0.139, # a narrow dip at z_dip=0.52 (subcritical exit). Note 1/(1+alpha)=1-sin2thW. # * crossing channel (eta_mix below): eta = 2*alpha^2 = 0.181, a BROAD dip # at Z_CRIT=1.588 (the M=1 transcritical crossing) — the DOMINANT term. # See desi-dark-energy-crust.qmd "The crust suppresses structure growth". eta_crust:Parameter = field(default_factory=lambda: Parameter("eta_crust", 0.139, 0.0, 0.25, 0.0005)) z_dip: Parameter = field(default_factory=lambda: Parameter("z_dip", 0.52, 0.1, 1.5, 0.01 )) sigma_sup:Parameter = field(default_factory=lambda: Parameter("sigma_sup", 0.30, 0.05, 1.0, 0.01 )) # --- Vortex mixing (the DOMINANT G_eff channel, at the M=1 crossing) --- mixing: Toggle = field(default_factory=lambda: Toggle("Include vortex mixing", True)) eta_mix: Parameter = field(default_factory=lambda: Parameter("eta_mix", 0.181, 0.0, 0.25, 0.001 )) sigma_mix:Parameter = field(default_factory=lambda: Parameter("sigma_mix", 0.56, 0.15, 1.2, 0.01 )) # --- Bulk deficit --- C: Parameter = field(default_factory=lambda: Parameter("C (depth)", 1.00, 0.0, 1.5, 0.02 )) z_b: Parameter = field(default_factory=lambda: Parameter("z_b (onset)", 2.20, 0.3, 3.0, 0.05 )) def _make_params(b, zs, gam, c, zb, eta, zd, sig, mo, em, sm, dsw, h0, om, zh, al) -> SimulationParams: """Map Pr(...) positional arguments onto a SimulationParams instance.""" p = SimulationParams() # Cosmology p.H0.value = h0 p.Om.value = om # Crust enhancement p.dsw.value = dsw p.B.value = b p.z_s.value = zs p.gamma.value = gam p.z_harm.value = zh p.alpha.value = al # Crust suppression p.eta_crust.value = eta p.z_dip.value = zd p.sigma_sup.value = sig # Vortex mixing p.mixing.value = mo p.eta_mix.value = em p.sigma_mix.value = sm # Bulk deficit p.C.value = c p.z_b.value = zb return p # ── Physical / survey constants ──────────────────────────────── Or = 9.1e-5 # Radiation density parameter RD_PLANCK = 147.09 # Sound horizon [Mpc], Planck 2018 H0_SHOES = 73.04 # SH0ES local H0 [km/s/Mpc] Z_CRIT = 1.588 # Transcritical crossing redshift Z_B_DSW = 2.20 # DSW soliton edge redshift # CPL / DESI best-fit reference values DW0, DWA = -0.42, -1.75 S0, SA = 0.21, 0.58 RH = -0.85 # Standard ΛCDM reference (Planck 2018, fixed — used for Jia H0 inference) OM_LCDM = 0.315 OL_LCDM = 0.685 # Jia et al. observed H0(z) data points JIA_Z = [0.1, 0.3, 0.5, 0.7, 2.5 ] JIA_H0 = [72.20, 71.62, 69.78, 68.13, 67.23] JIA_ERR = [0.19, 0.36, 0.51, 0.67, 0.84 ] # DESI Y3 BAO measurements — each entry: (z, type, value, sigma, label) # type: 'DV' | 'DM' | 'DH' DESI_BAO = [ dict(z=0.295, type='DV', val=7.9417, sig=0.0764, label='BGS'), dict(z=0.510, type='DM', val=13.5876, sig=0.1685, label='LRG1'), dict(z=0.510, type='DH', val=21.8629, sig=0.4296, label='LRG1'), dict(z=0.706, type='DM', val=17.3507, sig=0.1793, label='LRG2'), dict(z=0.706, type='DH', val=19.4553, sig=0.3338, label='LRG2'), dict(z=0.934, type='DM', val=21.5756, sig=0.1600, label='LRG3+ELG1'), dict(z=0.934, type='DH', val=17.6415, sig=0.2009, label='LRG3+ELG1'), dict(z=1.321, type='DM', val=27.6009, sig=0.3226, label='ELG2'), dict(z=1.321, type='DH', val=14.1760, sig=0.2253, label='ELG2'), dict(z=1.484, type='DM', val=30.5119, sig=0.7629, label='QSO'), dict(z=1.484, type='DH', val=12.8170, sig=0.5145, label='QSO'), dict(z=2.330, type='DH', val=8.6315, sig=0.1018, label='Lyα'), dict(z=2.330, type='DM', val=38.9890, sig=0.5318, label='Lyα'), ] Z_BAO = [0.30, 0.51, 0.71, 0.93, 1.32, 2.33] # Moraine spline knots (15-knot signed freeform fit, χ²=9.68) # # This is the ORIGINAL freeform fit: 15 free knots, natural-cubic interpolation # below. It is kept as-is for continuity — it is the "f_spline" curve the script # has always produced. SPLINE_Z = [0.00, 0.07, 0.15, 0.23, 0.30, 0.38, 0.45, 0.60, 0.80, 1.00, 1.30, 1.60, 2.00, 2.30, 2.50] SPLINE_A = [+0.1639, +0.2509, +0.0551, +0.7043, -0.3105, -0.3559, +0.2788, -0.1312, +0.7205, -0.1114, +0.1975, -0.0887, -0.3650, +0.4887, +0.2884] # ── Bore-skeleton spline (14-knot, PCHIP) ───────────────────── # The benchmark fit that the paper figure (arxiv/two-tensions.svg, via # scripts/tensions_figure.py) now uses for the "Moraine Crust Spline" curve. # # Produced by data/desi/bao_data/fit_bore_spline.py in --bore mode: # python fit_bore_spline.py --desi --jia --bore \ # --desi-mean … --desi-cov … --h0-local-target 72.0 \ # --n-starts 8 --max-iter 40 # best-fit record: results/bore-desi-jia-72h-{meta,fit-starts}.csv (rank 1), # χ²_total = 10.2138 (χ²_BAO 7.06 + χ²_Jia 3.15, full DESI covariance), # 14 knots, H₀-local target 72.0 km/s/Mpc. # # Knot z-grid is BORE_SKELETON_KNOTS from fit_bore_spline.py (z=1.588 is the # transcritical crossing). fit_bore_spline interpolates the knots with scipy's # PchipInterpolator (monotone cubic Hermite), NOT a natural cubic spline, so the # bore curve is reproduced here with a matching pure-Python PCHIP (bore_interp). BORE_Z = [0.00, 0.07, 0.15, 0.23, 0.33, 0.45, 0.60, 0.80, 1.00, 1.30, 1.588, 2.00, 2.30, 2.50] BORE_A = [+0.197956, +0.243262, +0.000000, +0.676863, -0.593041, +0.380286, -0.209047, +0.826247, -0.148742, +0.182902, -0.031837, -0.356463, +0.485128, +0.180061] # ── Numerical integration (Simpson's rule) ──────────────────── def simp_int(z: float, efn, N: int) -> float: """Integrate 1/√E(z') from 0 to z using Simpson's rule with N steps.""" if z <= 0: return 0.0 h = z / N s = 1.0 / math.sqrt(efn(0)) for i in range(1, N): w = 4 if (i % 2) else 2 s += w / math.sqrt(efn(i * h)) s += 1.0 / math.sqrt(efn(z)) return s * h / 3.0 # ── E² — Hubble rate squared (normalized) ───────────────────── def E2_lcdm(z: float, Om: float, OL: float) -> float: return Om * (1+z)**3 + Or * (1+z)**4 + OL def E2_lcdm_std(z: float) -> float: """Fixed Planck 2018 ΛCDM — used as the observer's assumed model in Jia.""" return OM_LCDM * (1+z)**3 + OL_LCDM def E2_cpl(z: float, w0: float, wa: float, Om: float, OL: float) -> float: """Chevallier-Polarski-Linder dark energy equation of state.""" return (Om * (1+z)**3 + Or * (1+z)**4 + OL * (1+z)**(3*(1+w0+wa)) * math.exp(-3*wa*z/(1+z))) def E2_sub(z: float, B: float, zs: float, C: float, zb: float, gamma: float, Om: float, OL: float, Z_HARM: float, ALPHA_TAIL: float) -> float: """Substrate model E²: ΛCDM + crust enhancement − bulk deficit.""" return (Om * (1+z)**3 + Or * (1+z)**4 + OL * f_mod(z, B, zs, C, zb, gamma, Om, OL, Z_HARM, ALPHA_TAIL, dsw=True)) # ── Comoving distance integrals ─────────────────────────────── def D_M_sub(z: float, B: float, zs: float, C: float, zb: float, gamma: float, Om: float, OL: float, Z_HARM: float, ALPHA_TAIL: float) -> float: return simp_int(z, lambda zz: E2_sub(zz, B, zs, C, zb, gamma, Om, OL, Z_HARM, ALPHA_TAIL), 80) def D_M_cpl(z: float, w0: float, wa: float, Om: float, OL: float) -> float: return simp_int(z, lambda zz: E2_cpl(zz, w0, wa, Om, OL), 60) def D_M_lcdm_std(z: float) -> float: return simp_int(z, E2_lcdm_std, 80) # ── Dark energy density fraction ────────────────────────────── def OmegaL(z: float, Om: float, OL: float) -> float: return OL / (Om * (1+z)**3 + OL) def OmegaL_crit(Om: float, OL: float) -> float: return OmegaL(Z_CRIT, Om, OL) # ── DSW (Dispersive Shock Wave) envelope ────────────────────── def e_DSW(z: float, Z_HARM: float, ALPHA_TAIL: float) -> float: """ DSW envelope function with three regimes: - z > Z_B_DSW : soliton edge (sech² tail) - Z_CRIT ≤ z ≤ Z_B : upstream body (√ profile, Grimshaw-Smyth) - Z_HARM ≤ z < Z_CRIT : downstream tail (power law, exponent α) """ if z > Z_B_DSW: return 1.0 / math.cosh((z - Z_B_DSW) / 0.15)**2 if Z_CRIT <= z <= Z_B_DSW: d = Z_B_DSW - Z_CRIT return math.sqrt((Z_B_DSW - z) / d) if d > 0 else 0.0 if Z_HARM <= z < Z_CRIT: d2 = Z_CRIT - Z_HARM return ((z - Z_HARM) / d2) ** ALPHA_TAIL if d2 > 0 else 0.0 return 0.0 def dsw_product_peak(gamma: float, Z_HARM: float, ALPHA_TAIL: float, Om: float, OL: float) -> float: """Max of eDSW(z) × (ΩΛ(z)/ΩΛ_crit)^γ — used to normalise fCrustDSW.""" mx = 0.0 olc = OmegaL_crit(Om, OL) z = max(0.0, Z_HARM) while z <= Z_B_DSW + 0.3: v = e_DSW(z, Z_HARM, ALPHA_TAIL) * (OmegaL(z, Om, OL) / olc) ** gamma if v > mx: mx = v z += 0.004 return mx if mx > 0 else 1.0 # ── Crust enhancement functions ─────────────────────────────── def f_crust_DSW(z: float, B: float, gamma: float, Z_HARM: float, ALPHA_TAIL: float, Om: float, OL: float) -> float: env = e_DSW(z, Z_HARM, ALPHA_TAIL) if env < 1e-10: return 0.0 peak = dsw_product_peak(gamma, Z_HARM, ALPHA_TAIL, Om, OL) return B * env * (OmegaL(z, Om, OL) / OmegaL_crit(Om, OL))**gamma / peak def df_crust_DSW(z: float, B: float, gamma: float, Z_HARM: float, ALPHA_TAIL: float, Om: float, OL: float) -> float: dz = 0.001 return (f_crust_DSW(z+dz, B, gamma, Z_HARM, ALPHA_TAIL, Om, OL) - f_crust_DSW(max(Z_HARM, z-dz), B, gamma, Z_HARM, ALPHA_TAIL, Om, OL)) / (2*dz) def f_crust_param(z: float, B: float, zs: float) -> float: """Parametric (non-DSW) crust profile.""" return B * z * math.exp(-z / zs) def df_crust_param(z: float, B: float, zs: float) -> float: return B * math.exp(-z / zs) * (1 - z / zs) def f_crust(z: float, B: float, zs: float, gamma: float, Z_HARM: float, ALPHA_TAIL: float, Om: float, OL: float, dsw: bool) -> float: if dsw: return f_crust_DSW(z, B, gamma, Z_HARM, ALPHA_TAIL, Om, OL) return f_crust_param(z, B, zs) def df_crust(z: float, B: float, zs: float, gamma: float, Z_HARM: float, ALPHA_TAIL: float, Om: float, OL: float, dsw: bool) -> float: if dsw: return df_crust_DSW(z, B, gamma, Z_HARM, ALPHA_TAIL, Om, OL) return df_crust_param(z, B, zs) # ── Bulk deficit (Volovik self-tuning) ──────────────────────── def g_bulk(z: float, zb: float) -> float: return z**2 / (z**2 + zb**2) def dg_bulk(z: float, zb: float) -> float: d = z**2 + zb**2 return 2 * z * zb**2 / d**2 # ── Combined modifier fMod = 1 + crust − bulk ───────────────── def f_mod(z: float, B: float, zs: float, C: float, zb: float, gamma: float, Om: float, OL: float, Z_HARM: float, ALPHA_TAIL: float, dsw: bool = True) -> float: return 1 + f_crust(z, B, zs, gamma, Z_HARM, ALPHA_TAIL, Om, OL, dsw) \ - C * g_bulk(z, zb) def df_mod(z: float, B: float, zs: float, C: float, zb: float, gamma: float, Om: float, OL: float, Z_HARM: float, ALPHA_TAIL: float, dsw: bool = True) -> float: return (df_crust(z, B, zs, gamma, Z_HARM, ALPHA_TAIL, Om, OL, dsw) - C * dg_bulk(z, zb)) def f_mod_spline(z: float, C: float, zb: float) -> float: return 1 + spline_interp(z) - C * g_bulk(z, zb) def E2_spline(z: float, C: float, zb: float, Om: float, OL: float) -> float: return Om*(1+z)**3 + Or*(1+z)**4 + OL*f_mod_spline(z, C, zb) def f_mod_bore(z: float, C: float, zb: float) -> float: """f_mod for the 14-knot bore-skeleton spline (PCHIP crust + bulk deficit). Same structure as f_mod_spline — a bulk deficit (1 − C·g_bulk) plus a crust overlay — but the crust is the bore-skeleton PCHIP curve (bore_interp), which mirrors fit_bore_spline.py's f_baseline + crust_spline exactly. """ return 1 + bore_interp(z) - C * g_bulk(z, zb) def E2_bore(z: float, C: float, zb: float, Om: float, OL: float) -> float: return Om*(1+z)**3 + Or*(1+z)**4 + OL*f_mod_bore(z, C, zb) # ── Full Friedmann self-consistency (WIP-18) ────────────────── # The bootstrap E2_sub holds Ω_Λ = 1−Ω_m−Ω_r FIXED and plugs in f(z). That # leaves three ΛCDM-reference leaks (open-problems.qmd#wip-18), all closed here: # (L1) FLATNESS — E²(0) = Ω_m + Ω_r + Ω_Λ·f(0) must equal 1. With a # present-epoch crust crest f(0)=1.25 the bootstrap runs at E²(0)=1.17, # a 9% violation of H(0)≡H0. Fix: Ω_Λ,base = (1−Ω_m−Ω_r)/f(0) — the # Volovik *base* density sits a factor f(0) below today's observed Ω_Λ. # (L2) DE-WEIGHT — the crust's [Ω_Λ(z)/Ω_Λ(z_crit)]^γ factor must use the # model's OWN Ω_Λ(z) = Ω_Λ,base·f(z)/E²(z), not ΛCDM's. # (L3) z_crit — re-solved from M(z) = H(z)·d_proper(z)/c = 1 under H(z). # This is a pure-Python port of friedmann_self_consistency.py's fixed-point, # specialised to the DSW Profile-C crust so the S8 background route (and the # distance integrals) can be evaluated on the self-consistent expansion law. def e_DSW_zc(z: float, Z_HARM: float, ALPHA_TAIL: float, z_crit: float, z_b: float) -> float: """DSW envelope with explicit (self-consistent) z_crit and z_b edges. Identical in form to e_DSW, but the upstream/downstream split is pinned at the passed z_crit (which L3 may move off the global Z_CRIT) and the soliton edge at the passed z_b. """ if z > z_b: return 1.0 / math.cosh((z - z_b) / 0.15)**2 if z_crit <= z <= z_b: d = z_b - z_crit return math.sqrt((z_b - z) / d) if d > 0 else 0.0 if Z_HARM <= z < z_crit: d2 = z_crit - Z_HARM return ((z - Z_HARM) / d2) ** ALPHA_TAIL if d2 > 0 else 0.0 return 0.0 def _solve_zcrit_grid(zg: list, fg: list, Om: float, OL_base: float, H0: float = None) -> float: """Re-solve the transcritical crossing M(z)=1 (L3) on the f(z) grid. M(z) = H(z)·d_proper(z)/c with d_proper = D_C/(1+z), D_C = (c/H0)∫dz'/H. The c and H0 cancel, leaving M(z) = E(z)·∫₀ᶻ dz'/E / (1+z). Returns the first crossing by linear interpolation (reproduces 1.588 under ΛCDM). """ n = len(zg) E = [math.sqrt(max(Om*(1+zg[i])**3 + Or*(1+zg[i])**4 + OL_base*fg[i], 1e-30)) for i in range(n)] inv = [1.0 / e for e in E] Dc = 0.0 s_prev = None for i in range(n): if i > 0: Dc += 0.5 * (inv[i] + inv[i-1]) * (zg[i] - zg[i-1]) s = E[i] * Dc / (1 + zg[i]) - 1.0 if s_prev is not None and s_prev <= 0 < s: z0, z1 = zg[i-1], zg[i] return z0 - s_prev * (z1 - z0) / (s - s_prev) s_prev = s return Z_CRIT def solve_self_consistent(B: float, gamma: float, C: float, zb: float, Om: float, OL: float, Z_HARM: float, ALPHA_TAIL: float, H0: float = H0_SHOES, n_iter: int = 80, tol: float = 1e-11) -> dict: """Fixed-point solve for (Ω_Λ,base, z_crit, Ω_Λ(z)) closing L1/L2/L3. Only meaningful for the DSW crust profile (dsw=True). Returns a dict with the converged scalars and self-consistent closures f(z), E2(z) bound to the converged state — drop-in replacements for f_mod / E2_sub on the flatness-normalised expansion law. """ OLtot = 1.0 - Om - Or NG = 600 zmax = max(zb, 3.0) + 0.5 zg = [zmax * i / NG for i in range(NG + 1)] z0b = max(0.0, Z_HARM) # seed: ΛCDM DE fraction, z_crit at the Planck value olz = [OL / (Om * (1+z)**3 + OL) for z in zg] zc = Z_CRIT OL_base = OLtot def _peak(olz_tab, zc_local, ol_crit): """Peak of envelope·DE-weight over the active band [max(0,z_harm), z_b].""" NB, pk = 400, 0.0 for j in range(NB + 1): zz = z0b + (zb - z0b) * j / NB e = e_DSW_zc(zz, Z_HARM, ALPHA_TAIL, zc_local, zb) ov = (_lerp_curve(zz, zg, olz_tab) / ol_crit) ** gamma v = e * ov if v > pk: pk = v return pk if pk > 1e-12 else 1.0 it = 0 for it in range(n_iter): ol_crit = _lerp_curve(zc, zg, olz) peak = _peak(olz, zc, ol_crit) # f(0) → flatness rescale of the base density (L1) env0 = e_DSW_zc(0.0, Z_HARM, ALPHA_TAIL, zc, zb) crust0 = B * env0 * (olz[0] / ol_crit) ** gamma / peak if env0 > 1e-12 else 0.0 f0 = 1.0 + crust0 - C * g_bulk(0.0, zb) OL_base = OLtot / f0 # f(z), E²(z), and the model's own Ω_Λ(z) on the grid (L2) fg, olz_new = [], [] for i, z in enumerate(zg): env = e_DSW_zc(z, Z_HARM, ALPHA_TAIL, zc, zb) crust = B * env * (olz[i] / ol_crit) ** gamma / peak if env > 1e-12 else 0.0 f = 1.0 + crust - C * g_bulk(z, zb) e2 = max(Om*(1+z)**3 + Or*(1+z)**4 + OL_base*f, 1e-30) fg.append(f) olz_new.append(OL_base * f / e2) zc_new = _solve_zcrit_grid(zg, fg, Om, OL_base, H0) # L3 change = max(abs(olz_new[i] - olz[i]) for i in range(len(zg))) + abs(zc_new - zc) olz, zc = olz_new, zc_new if change < tol: break # bind converged state into the closures olz_f, zc_f, olb_f = list(olz), zc, OL_base ol_crit_f = _lerp_curve(zc_f, zg, olz_f) peak_f = _peak(olz_f, zc_f, ol_crit_f) def f_sc(z: float) -> float: env = e_DSW_zc(z, Z_HARM, ALPHA_TAIL, zc_f, zb) crust = (B * env * (_lerp_curve(z, zg, olz_f) / ol_crit_f) ** gamma / peak_f if env > 1e-12 else 0.0) return 1.0 + crust - C * g_bulk(z, zb) def E2_sc(z: float) -> float: return Om*(1+z)**3 + Or*(1+z)**4 + olb_f * f_sc(z) return dict(OL_base=olb_f, z_crit=zc_f, f0=f_sc(0.0), iters=it + 1, E2_0=E2_sc(0.0), f=f_sc, E2=E2_sc, olz_z=zg, olz=olz_f) # ── Geff suppression (crust + vortex mixing) ────────────────── def g_eff(z: float, eta: float, zd: float, sig: float) -> float: """Gaussian suppression of G_eff around z_dip.""" return 1 - eta * math.exp(-0.5 * ((z - zd) / sig)**2) def dg_eff(z: float, eta: float, zd: float, sig: float) -> float: return eta * (z - zd) / sig**2 * math.exp(-0.5 * ((z - zd) / sig)**2) def eta_mix(z: float, em: float, sm: float) -> float: """Vortex mixing suppression centred on Z_CRIT.""" return em * math.exp(-0.5 * ((z - Z_CRIT) / sm)**2) def g_eff_total(z: float, eta: float, zd: float, sig: float, mixing_on: bool, em: float, sm: float) -> float: g = g_eff(z, eta, zd, sig) if mixing_on: g -= eta_mix(z, em, sm) return max(g, 0.5) # floor at 0.5 — physical lower bound def dg_eff_total(z: float, eta: float, zd: float, sig: float, mixing_on: bool, em: float, sm: float) -> float: dg = dg_eff(z, eta, zd, sig) if mixing_on: dg += em * (z - Z_CRIT) / sm**2 * math.exp(-0.5 * ((z - Z_CRIT) / sm)**2) return dg # ── Effective dark energy EOS and its derivative ────────────── def f_eff(z: float, B: float, zs: float, C: float, zb: float, gamma: float, eta: float, zd: float, sig: float, mixing_on: bool, em: float, sm: float, Om: float, OL: float, Z_HARM: float, ALPHA_TAIL: float) -> float: g = g_eff_total(z, eta, zd, sig, mixing_on, em, sm) f = f_mod(z, B, zs, C, zb, gamma, Om, OL, Z_HARM, ALPHA_TAIL) a = Om * (1+z)**3 / OL return (g - 1) * a + g * f def df_eff(z: float, B: float, zs: float, C: float, zb: float, gamma: float, eta: float, zd: float, sig: float, mixing_on: bool, em: float, sm: float, Om: float, OL: float, Z_HARM: float, ALPHA_TAIL: float) -> float: g = g_eff_total(z, eta, zd, sig, mixing_on, em, sm) dg = dg_eff_total(z, eta, zd, sig, mixing_on, em, sm) f = f_mod(z, B, zs, C, zb, gamma, Om, OL, Z_HARM, ALPHA_TAIL) df = df_mod(z, B, zs, C, zb, gamma, Om, OL, Z_HARM, ALPHA_TAIL) a = Om * (1+z)**3 / OL da = 3 * Om * (1+z)**2 / OL return dg * (a + f) + (g - 1) * da + g * df def w_eff(z: float, B: float, zs: float, C: float, zb: float, gamma: float, eta: float, zd: float, sig: float, mixing_on: bool, em: float, sm: float, Om: float, OL: float, Z_HARM: float, ALPHA_TAIL: float) -> float: """Effective dark energy equation-of-state parameter w(z).""" fe = f_eff(z, B, zs, C, zb, gamma, eta, zd, sig, mixing_on, em, sm, Om, OL, Z_HARM, ALPHA_TAIL) if abs(fe) < 0.01: return -1.0 dfe = df_eff(z, B, zs, C, zb, gamma, eta, zd, sig, mixing_on, em, sm, Om, OL, Z_HARM, ALPHA_TAIL) return -1 + (1+z) * dfe / (3 * fe) # ── Inferred H₀(z) — Jia et al. observer framework ─────────── def H0_inferred(z: float, B: float, zs: float, C: float, zb: float, gamma: float, Om: float, OL: float, Z_HARM: float, ALPHA_TAIL: float, H0_BASE: float) -> float: """ An observer assuming standard ΛCDM measures luminosity distances from the substrate model and infers H₀(z) = H0_BASE × I_ΛCDM / I_sub. Reproduces the Jia et al. running-H₀ diagnostic. """ if z < 0.005: z = 0.005 I_sub = simp_int(z, lambda zz: E2_sub(zz, B, zs, C, zb, gamma, Om, OL, Z_HARM, ALPHA_TAIL), 80) I_lcdm = simp_int(z, E2_lcdm_std, 80) if I_sub < 1e-10: return H0_BASE return H0_BASE * I_lcdm / I_sub def chi_sq_jia(B: float, zs: float, C: float, zb: float, gamma: float, Om: float, OL: float, Z_HARM: float, ALPHA_TAIL: float, H0_BASE: float) -> float: """χ² against Jia et al. running-H₀ observations.""" c2 = 0.0 for i in range(len(JIA_Z)): hm = H0_inferred(JIA_Z[i], B, zs, C, zb, gamma, Om, OL, Z_HARM, ALPHA_TAIL, H0_BASE) r = (hm - JIA_H0[i]) / JIA_ERR[i] c2 += r * r return c2 def chi_sq_jia_spline(C: float, zb: float, Om: float, OL: float, H0_BASE: float) -> float: """χ² against Jia et al. running-H₀ observations — spline model.""" c2 = 0.0 for i in range(len(JIA_Z)): hm = H0_inferred_spline(JIA_Z[i], C, zb, Om, OL, H0_BASE) r = (hm - JIA_H0[i]) / JIA_ERR[i] c2 += r * r return c2 def chi_sq_jia_bore(C: float, zb: float, Om: float, OL: float, H0_BASE: float) -> float: """χ² against Jia et al. running-H₀ observations — bore-skeleton model.""" c2 = 0.0 for i in range(len(JIA_Z)): hm = H0_inferred_bore(JIA_Z[i], C, zb, Om, OL, H0_BASE) r = (hm - JIA_H0[i]) / JIA_ERR[i] c2 += r * r return c2 def chi_sq_jia_lcdm(Om: float, OL: float, H0_BASE: float) -> float: """χ² against Jia et al. running-H₀ observations — flat ΛCDM. An ΛCDM observer fitting ΛCDM data infers the same H₀ at every z, so the only residual comes from H0_BASE ≠ JIA_H0[i]. """ c2 = 0.0 for i in range(len(JIA_Z)): z = JIA_Z[i] if z < 0.005: z = 0.005 I_model = simp_int(z, lambda zz: E2_lcdm(zz, Om, OL), 80) I_lcdm = simp_int(z, E2_lcdm_std, 80) hm = H0_BASE * I_lcdm / I_model if I_model > 1e-10 else H0_BASE r = (hm - JIA_H0[i]) / JIA_ERR[i] c2 += r * r return c2 # ── BAO predictions ─────────────────────────────────────────── def _bao_predict_single(om: float, ol: float, h0: float, e2_fn, n_int: int = 300) -> list[float]: """ Core BAO predictor — takes any E²(z) callable. Returns predicted DM/DH/DV (in units of r_d) for each DESI_BAO entry. """ cf = 299792.458 / (h0 * RD_PLANCK) preds = [] for d in DESI_BAO: z = d['z'] e2 = max(e2_fn(z), 1e-10) dh_rd = cf / math.sqrt(e2) dm_rd = cf * simp_int(z, lambda zz: max(e2_fn(zz), 1e-10), n_int) dv_rd = (z * dh_rd * dm_rd**2) ** (1/3) if d['type'] == 'DM': preds.append(dm_rd) elif d['type'] == 'DH': preds.append(dh_rd) else: preds.append(dv_rd) return preds def bao_predict_dsw(B: float, zs: float, C: float, zb: float, gamma: float, Om: float, OL: float, Z_HARM: float, ALPHA_TAIL: float, H0: float) -> list[float]: """BAO predictions from the full substrate/DSW model.""" return _bao_predict_single( Om, OL, H0, lambda z: E2_sub(z, B, zs, C, zb, gamma, Om, OL, Z_HARM, ALPHA_TAIL)) def bao_predict_lcdm(Om: float, OL: float, H0: float) -> list[float]: """BAO predictions from flat ΛCDM (f_mod = 1 everywhere).""" return _bao_predict_single( Om, OL, H0, lambda z: E2_lcdm_std(z)) def bao_predict_spline(C: float, zb: float, Om: float, OL: float, H0: float) -> list[float]: """BAO predictions from the spline moraine model.""" return _bao_predict_single( Om, OL, H0, lambda z: E2_spline(z, C, zb, Om, OL)) def bao_predict_bore(C: float, zb: float, Om: float, OL: float, H0: float) -> list[float]: """BAO predictions from the 14-knot bore-skeleton moraine model.""" return _bao_predict_single( Om, OL, H0, lambda z: E2_bore(z, C, zb, Om, OL)) # ── BAO χ² (diagonal — no off-diagonal covariance) ────────── def bao_chi2_diagonal(preds: list[float]) -> float: """Diagonal BAO χ² — sum of squared pulls against DESI_BAO. Note: fit_freeform_spline.py uses a full covariance matrix loaded from external files, which accounts for off-diagonal correlations between observations at the same effective redshift (e.g. DM and DH at z=0.510). This diagonal version ignores those correlations because DESI_BAO in this file stores only the marginalised 1-sigma uncertainties. The diagonal result will differ from the full- covariance result by a modest amount. """ return sum(((preds[k] - d['val']) / d['sig'])**2 for k, d in enumerate(DESI_BAO)) # ── Combined total χ² = BAO + Jia (matches fit_freeform_spline.py) ── def total_chi2_dsw(B: float, zs: float, C: float, zb: float, gamma: float, Om: float, OL: float, Z_HARM: float, ALPHA_TAIL: float, H0: float) -> dict: """Total χ² for the full substrate / DSW model: BAO + Jia.""" preds = bao_predict_dsw(B, zs, C, zb, gamma, Om, OL, Z_HARM, ALPHA_TAIL, H0) c2_bao = bao_chi2_diagonal(preds) c2_jia = chi_sq_jia(B, zs, C, zb, gamma, Om, OL, Z_HARM, ALPHA_TAIL, H0) return {'chi2_bao': c2_bao, 'chi2_jia': c2_jia, 'chi2_total': c2_bao + c2_jia} def total_chi2_spline(C: float, zb: float, Om: float, OL: float, H0: float) -> dict: """Total χ² for the spline moraine model: BAO + Jia.""" preds = bao_predict_spline(C, zb, Om, OL, H0) c2_bao = bao_chi2_diagonal(preds) c2_jia = chi_sq_jia_spline(C, zb, Om, OL, H0) return {'chi2_bao': c2_bao, 'chi2_jia': c2_jia, 'chi2_total': c2_bao + c2_jia} def total_chi2_bore(C: float, zb: float, Om: float, OL: float, H0: float) -> dict: """Total χ² for the 14-knot bore-skeleton moraine model: BAO + Jia. NB: as with total_chi2_spline, this uses the DIAGONAL BAO covariance (bao_chi2_diagonal), so the χ²_BAO here differs from the full-covariance χ²_BAO=7.06 quoted in results/bore-desi-jia-72h-meta.csv (which was produced by fit_bore_spline.py with the off-diagonal DESI covariance). """ preds = bao_predict_bore(C, zb, Om, OL, H0) c2_bao = bao_chi2_diagonal(preds) c2_jia = chi_sq_jia_bore(C, zb, Om, OL, H0) return {'chi2_bao': c2_bao, 'chi2_jia': c2_jia, 'chi2_total': c2_bao + c2_jia} def total_chi2_lcdm(Om: float, OL: float, H0: float) -> dict: """Total χ² for flat ΛCDM: BAO + Jia.""" preds = bao_predict_lcdm(Om, OL, H0) c2_bao = bao_chi2_diagonal(preds) c2_jia = chi_sq_jia_lcdm(Om, OL, H0) return {'chi2_bao': c2_bao, 'chi2_jia': c2_jia, 'chi2_total': c2_bao + c2_jia} # ── CPL best-fit (grid search over w0, wa) ──────────────────── def fit_CPL(B: float, zs: float, C: float, zb: float, gamma: float, Om: float, OL: float, Z_HARM: float, ALPHA_TAIL: float) -> dict: """ Find best-fit CPL (w0, wa) by minimising distance residuals against the substrate comoving distances at Z_BAO redshifts. Four-pass grid refinement. """ Dt = [D_M_sub(z, B, zs, C, zb, gamma, Om, OL, Z_HARM, ALPHA_TAIL) for z in Z_BAO] def chi(w0, wa): s = 0.0 for k, z in enumerate(Z_BAO): dc = D_M_cpl(z, w0, wa, Om, OL) r = (Dt[k] - dc) / Dt[k] s += r * r return s # Pass 1 — coarse bw, ba, bc = -0.5, -1.5, 1e10 for i in range(31): for j in range(41): w0 = -2.0 + i * 2.5/30 wa = -5.0 + j * 7.0/40 c2 = chi(w0, wa) if c2 < bc: bc, bw, ba = c2, w0, wa # Pass 2 — medium for dw, da, steps in [(0.12, 0.25, 24), (0.02, 0.04, 24), (0.003, 0.006, 20)]: for i in range(steps+1): for j in range(steps+1): w0 = bw - dw + i * 2*dw/steps wa = ba - da + j * 2*da/steps c2 = chi(w0, wa) if c2 < bc: bc, bw, ba = c2, w0, wa return {'w0': round(bw, 4), 'wa': round(ba, 3)} # ── χ² against DESI CPL contours ───────────────────────────── def chi_sq_cpl(cpl: dict) -> float: """Tension with DESI Y3 CPL best fit, using covariance (S0, SA, RH).""" d0 = S0**2 * SA**2 * (1 - RH**2) C00 = SA**2 / d0 C11 = S0**2 / d0 C01 = -RH * S0 * SA / d0 dw0 = cpl['w0'] - DW0 dwa = cpl['wa'] - DWA return C00*dw0**2 + 2*C01*dw0*dwa + C11*dwa**2 # ── S8 growth parameter ─────────────────────────────────────── def compute_S8(B: float, zs: float, C: float, zb: float, gamma: float, eta: float, zd: float, sig: float, mixing_on: bool, em: float, sm: float, Om: float, OL: float, Z_HARM: float, ALPHA_TAIL: float, e2_sub_fn=None) -> float: """ Integrate the linear growth ODE in log-a to get S8 = 0.832 × D/D_ΛCDM. Uses 500 steps (JS used 500 for the standard model, 2000 for spline variant). The crust enters growth through TWO channels (see desi-dark-energy-crust.qmd "The crust suppresses structure growth"): • the BACKGROUND route — the substrate E²(z) sets the Hubble friction (2+dlnH) and the Ω_m(1+z)³/E² source; this is the ONLY place the crust's f(z) touches growth, and it nearly cancels between the low-z crest and the high-z deficit. • the COUPLING route — G_eff(z) multiplies the source only (force coupling, NOT energy content), and does essentially all the suppression. `e2_sub_fn` overrides the substrate E²(z): pass a self-consistent (flatness- normalised) closure from solve_self_consistent to run the background route on the full-Friedmann expansion law (WIP-18). G_eff is untouched — it does not belong in Friedmann. Default (None) uses the bootstrap E2_sub. """ sub_e2 = (e2_sub_fn if e2_sub_fn is not None else (lambda zz: E2_sub(zz, B, zs, C, zb, gamma, Om, OL, Z_HARM, ALPHA_TAIL))) N = 500 la0 = -math.log(31) # a_start = 1/31 (z~30) dl = (0 - la0) / N D, Dp = 1e-5, 1e-5 # substrate growth and its derivative Dl, Dpl = 1e-5, 1e-5 # ΛCDM growth and its derivative for i in range(N): la = la0 + i * dl z = max(1/math.exp(la) - 1, 0) e2 = max(sub_e2(z), 1e-10) e2l = max(E2_lcdm(z, Om, OL), 1e-10) g = g_eff_total(z, eta, zd, sig, mixing_on, em, sm) # Numerical derivative of ln E² (centred, half-step offset). # NB: dlH is d ln H/d ln a, NOT d ln H²/d ln a. The samples are dl apart # (la ± dl/2), so the centred difference divides by dl; the *extra* factor # of 2 in (2*dl) is the E²→H Jacobian (d ln H = ½ d ln E²). The two are # unrelated and must not be "simplified" away — together they give exactly # the d ln H/d ln a that the friction term (2 + dlH) requires. zp = max(1/math.exp(la + dl*0.5) - 1, 0) zm = max(1/math.exp(la - dl*0.5) - 1, 0) e2p = max(sub_e2(zp), 1e-10) e2m = max(sub_e2(zm), 1e-10) dlH = (math.log(e2p) - math.log(e2m)) / (2 * dl) e2pl = max(E2_lcdm(zp, Om, OL), 1e-10) e2ml = max(E2_lcdm(zm, Om, OL), 1e-10) dlHl = (math.log(e2pl) - math.log(e2ml)) / (2 * dl) # Substrate growth ODE src = 1.5 * g * Om * (1+z)**3 / e2 fr = 2 + dlH Dpp = src * D - fr * Dp D += Dp * dl Dp += Dpp * dl # ΛCDM growth ODE (reference) srcl = 1.5 * Om * (1+z)**3 / e2l frl = 2 + dlHl Dppl = srcl * Dl - frl * Dpl Dl += Dpl * dl Dpl += Dppl * dl sigma8 = 0.811 * D / Dl return sigma8 * math.sqrt(Om / 0.3) # ── Moraine spline (natural cubic, 15 knots) ────────────────── def _build_spline_coeffs() -> dict: """Compute natural cubic spline coefficients for the moraine model.""" n = len(SPLINE_Z) - 1 h = [SPLINE_Z[i+1] - SPLINE_Z[i] for i in range(n)] al = [0.0] * (n+1) for i in range(1, n): al[i] = (3/h[i] * (SPLINE_A[i+1] - SPLINE_A[i]) - 3/h[i-1] * (SPLINE_A[i] - SPLINE_A[i-1])) l, mu, zc = [1.0]+[0.0]*n, [0.0]*(n+1), [0.0]*(n+1) for i in range(1, n): l[i] = 2*(SPLINE_Z[i+1] - SPLINE_Z[i-1]) - h[i-1]*mu[i-1] mu[i] = h[i] / l[i] zc[i] = (al[i] - h[i-1]*zc[i-1]) / l[i] l[n] = 1.0 a = list(SPLINE_A) c = [0.0] * (n+1) b = [0.0] * n d = [0.0] * n for j in range(n-1, -1, -1): c[j] = zc[j] - mu[j]*c[j+1] b[j] = (a[j+1]-a[j])/h[j] - h[j]*(c[j+1] + 2*c[j])/3 d[j] = (c[j+1] - c[j]) / (3*h[j]) return dict(a=a, b=b, c=c, d=d, n=n) _SP_COEFFS = None # lazily built on first use def spline_interp(z: float) -> float: global _SP_COEFFS if _SP_COEFFS is None: _SP_COEFFS = _build_spline_coeffs() sp = _SP_COEFFS if z <= SPLINE_Z[0]: return SPLINE_A[0] if z >= SPLINE_Z[sp['n']]: # Exponential tail-off beyond last knot return SPLINE_A[sp['n']] * math.exp(-(z - SPLINE_Z[sp['n']]) * 8) i = next(k for k in range(sp['n']-1, -1, -1) if z >= SPLINE_Z[k]) dz = z - SPLINE_Z[i] return sp['a'][i] + sp['b'][i]*dz + sp['c'][i]*dz**2 + sp['d'][i]*dz**3 # ── Bore-skeleton spline (PCHIP — monotone cubic Hermite) ───── # fit_bore_spline.py interpolates its knots with scipy.interpolate.PchipInterpolator, # so the bore curve must be reproduced with the SAME scheme, not the natural cubic # above. This is a pure-Python port of PCHIP: Fritsch-Carlson interior slopes # (harmonic-weighted secants, zeroed at local extrema) with scipy's one-sided # three-point endpoint rule, then piecewise Hermite evaluation. Verified to # reproduce results/bore-desi-jia-72h-fit-curve.csv to < 1e-6. def _sign(x: float) -> float: return (x > 0) - (x < 0) def _build_pchip_slopes(x: list, y: list) -> list: n = len(x) h = [x[i+1] - x[i] for i in range(n-1)] m = [(y[i+1] - y[i]) / h[i] for i in range(n-1)] # secant slopes d = [0.0] * n for k in range(1, n-1): if m[k-1] * m[k] <= 0: d[k] = 0.0 # local extremum → flat else: w1 = 2*h[k] + h[k-1] w2 = h[k] + 2*h[k-1] d[k] = (w1 + w2) / (w1/m[k-1] + w2/m[k]) # weighted harmonic mean def _edge(h0, h1, m0, m1): d0 = ((2*h0 + h1)*m0 - h0*m1) / (h0 + h1) if _sign(d0) != _sign(m0): return 0.0 if _sign(m0) != _sign(m1) and abs(d0) > 3*abs(m0): return 3*m0 return d0 d[0] = _edge(h[0], h[1], m[0], m[1]) d[n-1] = _edge(h[n-2], h[n-3], m[n-2], m[n-3]) return d _BORE_SLOPES = None # lazily built on first use def bore_interp(z: float) -> float: """PCHIP crust overlay for the 14-knot bore skeleton (BORE_Z, BORE_A). extrapolate=False (matches fit_bore_spline): outside [BORE_Z[0], BORE_Z[-1]] the crust is 0 (the baseline bulk deficit still applies in f_mod_bore). """ global _BORE_SLOPES if _BORE_SLOPES is None: _BORE_SLOPES = _build_pchip_slopes(BORE_Z, BORE_A) d = _BORE_SLOPES if z < BORE_Z[0] or z > BORE_Z[-1]: return 0.0 if z >= BORE_Z[-1]: return BORE_A[-1] i = next(k for k in range(len(BORE_Z)-1) if BORE_Z[k] <= z < BORE_Z[k+1]) h = BORE_Z[i+1] - BORE_Z[i] t = (z - BORE_Z[i]) / h t2, t3 = t*t, t*t*t h00 = 2*t3 - 3*t2 + 1 h10 = t3 - 2*t2 + t h01 = -2*t3 + 3*t2 h11 = t3 - t2 return (h00*BORE_A[i] + h10*h*d[i] + h01*BORE_A[i+1] + h11*h*d[i+1]) # ── Phantom crossing finder ─────────────────────────────────── def find_phantom_crossing(B: float, zs: float, C: float, zb: float, gamma: float, eta: float, zd: float, sig: float, mixing_on: bool, em: float, sm: float, Om: float, OL: float, Z_HARM: float, ALPHA_TAIL: float) -> Optional[str]: """Return redshift (as string) where w(z) crosses −1, or None.""" z = 0.01 w_prev = w_eff(z-0.005, B, zs, C, zb, gamma, eta, zd, sig, mixing_on, em, sm, Om, OL, Z_HARM, ALPHA_TAIL) while z < 3.0: w_cur = w_eff(z, B, zs, C, zb, gamma, eta, zd, sig, mixing_on, em, sm, Om, OL, Z_HARM, ALPHA_TAIL) if (w_prev + 1) * (w_cur + 1) < 0: t = abs(w_prev+1) / (abs(w_prev+1) + abs(w_cur+1)) return f"{z - 0.005 + t*0.005:.2f}" w_prev = w_cur z += 0.005 return None def createLCDM() -> SimulationParams: """ΛCDM preset — no crust, no mixing, flat substrate.""" return _make_params( b=0.0, zs=0.65, gam=1.90, c=0.0, zb=1.0, eta=0.0, zd=0.52, sig=0.30, mo=False, em=0.181, sm=0.56, dsw=True, h0=67.4, om=0.315, zh=0.0, al=0.50 ) def createSubstrate() -> SimulationParams: """★ Substrate preset — full DSW Profile C with vortex mixing.""" return _make_params( b=0.25, zs=0.65, gam=3.64, c=1.00, zb=2.20, eta=0.139, zd=0.52, sig=0.30, mo=True, em=0.181, sm=0.56, dsw=True, h0=67.4, om=0.315, zh=-0.25, al=0.50 ) @dataclass class EnergyDensityCurve: """ Output of SubstrateModel.compute_energy_density(). All arrays are parallel and indexed by redshift step. """ # Redshift axis z: list[float] # Primary curve — f_mod(z) or 1.0 (ΛCDM flat reference) f_mod: list[float] # Spline moraine overlay — f_mod_spline(z), always computed f_spline: list[float] # Spline knot markers for overlay rendering knot_z: list[float] knot_f: list[float] knot_amplitude: list[float] # raw SPLINE_A value — sign encodes ridge vs void # Scalar annotations f_at_zero: float # f_mod(z=0) — marks the crust peak z0: float # left edge of redshift axis (may be negative) z1: float # right edge class SubstrateModel: """ Pure physics model — no drawing, no UI. Wraps a SimulationParams instance and exposes compute methods that return plain data structures ready for plotting or analysis. Usage ----- params = createSubstrate() model = SubstrateModel(params) curve = model.compute_energy_density(steps=200) print(curve.f_at_zero) print(curve.f_mod) """ # Physical range constants (match the JS drawF window) F_MIN = 0.0 F_MAX = 2.2 Z_END = 3.0 def __init__(self, params): self.p = params # ── Convenience accessors ────────────────────────────────── @property def Om(self) -> float: return self.p.Om.value @property def OL(self) -> float: return 1.0 - self.p.Om.value - Or @property def H0(self) -> float: return self.p.H0.value @property def Z_HARM(self) -> float: return self.p.z_harm.value @property def ALPHA_TAIL(self) -> float: return self.p.alpha.value def _params_tuple(self): """Unpack the physics kwargs used by most compute functions.""" p = self.p return dict( B = p.B.value, zs = p.z_s.value, C = p.C.value, zb = p.z_b.value, gamma = p.gamma.value, Om = self.Om, OL = self.OL, Z_HARM = self.Z_HARM, ALPHA_TAIL = self.ALPHA_TAIL, dsw = self.p.dsw.value, ) # ── Primary compute methods ──────────────────────────────── def compute_energy_density(self, steps: int = 200) -> EnergyDensityCurve: """ Compute f_mod(z) and the spline moraine overlay across [z0, Z_END]. Parameters ---------- steps : int Number of evenly-spaced redshift samples. 200–250 matches the original dashboard resolution; increase for higher-precision output. Returns ------- EnergyDensityCurve Parallel arrays of z, f_mod, f_spline, plus knot markers. """ pk = self._params_tuple() z0 = min(0.0, self.Z_HARM - 0.05) z1 = self.Z_END # ── Main f_mod curve ────────────────────────────────── z_axis = [z0 + (z1 - z0) * i / steps for i in range(steps + 1)] f_main = [ f_mod(z, pk['B'], pk['zs'], pk['C'], pk['zb'], pk['gamma'], pk['Om'], pk['OL'], pk['Z_HARM'], pk['ALPHA_TAIL'], pk['dsw']) for z in z_axis ] # ── f(z=0) annotation ───────────────────────────────── f_at_zero = f_mod(0.0, pk['B'], pk['zs'], pk['C'], pk['zb'], pk['gamma'], pk['Om'], pk['OL'], pk['Z_HARM'], pk['ALPHA_TAIL'], pk['dsw']) # ── Spline moraine overlay ──────────────────────────── # Uses a clipped z-axis starting at 0 (spline undefined for z<0) z0_spl = max(0.0, z0) z_spl = [z0_spl + (z1 - z0_spl) * i / 250 for i in range(251)] f_spl_raw = [f_mod_spline(z, pk['C'], pk['zb']) for z in z_spl] # Clip to the visible f window (mirrors the JS bounds check) f_spline = [ f if self.F_MIN <= f <= self.F_MAX else None for f in f_spl_raw ] # ── Spline knot markers ─────────────────────────────── knot_z, knot_f, knot_amp = [], [], [] for k, zk in enumerate(SPLINE_Z): if z0 <= zk <= z1: fk = f_mod_spline(zk, pk['C'], pk['zb']) if self.F_MIN <= fk <= self.F_MAX: knot_z.append(zk) knot_f.append(fk) knot_amp.append(SPLINE_A[k]) return EnergyDensityCurve( z = z_axis, f_mod = f_main, f_spline = f_spline, knot_z = knot_z, knot_f = knot_f, knot_amplitude = knot_amp, f_at_zero = f_at_zero, z0 = z0, z1 = z1, ) def compute_H0_running(self, steps: int = 200) -> dict: """ Inferred H₀(z) across the Jia redshift range [0.005, 2.5]. Returns {'z': [...], 'H0': [...], 'chi_sq': float} """ pk = self._params_tuple() z_ax = [max(0.005, 2.5 * i / steps) for i in range(1, steps + 1)] h0_vals = [ H0_inferred(z, pk['B'], pk['zs'], pk['C'], pk['zb'], pk['gamma'], pk['Om'], pk['OL'], pk['Z_HARM'], pk['ALPHA_TAIL'], self.H0) for z in z_ax ] chi2 = chi_sq_jia(pk['B'], pk['zs'], pk['C'], pk['zb'], pk['gamma'], pk['Om'], pk['OL'], pk['Z_HARM'], pk['ALPHA_TAIL'], self.H0) return {'z': z_ax, 'H0': h0_vals, 'chi_sq': chi2} def compute_w_eff(self, steps: int = 200) -> dict: """ Effective dark energy EOS w(z) and phantom crossing redshift. Returns {'z': [...], 'w': [...], 'phantom_crossing': str | None} """ pk = self._params_tuple() p = self.p z_ax = [self.Z_END * i / steps for i in range(steps + 1)] w_vals = [ w_eff(z, pk['B'], pk['zs'], pk['C'], pk['zb'], pk['gamma'], p.eta_crust.value, p.z_dip.value, p.sigma_sup.value, p.mixing.value, p.eta_mix.value, p.sigma_mix.value, pk['Om'], pk['OL'], pk['Z_HARM'], pk['ALPHA_TAIL']) for z in z_ax ] crossing = find_phantom_crossing( pk['B'], pk['zs'], pk['C'], pk['zb'], pk['gamma'], p.eta_crust.value, p.z_dip.value, p.sigma_sup.value, p.mixing.value, p.eta_mix.value, p.sigma_mix.value, pk['Om'], pk['OL'], pk['Z_HARM'], pk['ALPHA_TAIL'] ) return {'z': z_ax, 'w': w_vals, 'phantom_crossing': crossing} def compute_S8(self) -> float: """S8 growth suppression scalar (bootstrap background E²).""" pk = self._params_tuple() p = self.p return compute_S8( pk['B'], pk['zs'], pk['C'], pk['zb'], pk['gamma'], p.eta_crust.value, p.z_dip.value, p.sigma_sup.value, p.mixing.value, p.eta_mix.value, p.sigma_mix.value, pk['Om'], pk['OL'], pk['Z_HARM'], pk['ALPHA_TAIL'] ) def solve_self_consistency(self) -> dict: """Converged full-Friedmann state (L1/L2/L3) for the DSW crust. Returns the dict from solve_self_consistent (Ω_Λ,base, z_crit, f(0), E²(0)=1, plus f/E2 closures). Only meaningful with dsw=True. """ pk = self._params_tuple() return solve_self_consistent( pk['B'], pk['gamma'], pk['C'], pk['zb'], pk['Om'], pk['OL'], pk['Z_HARM'], pk['ALPHA_TAIL'], self.H0) def compute_S8_sc(self, sc: dict = None) -> float: """S8 with the BACKGROUND route on the self-consistent expansion law. Closes the flatness/DE-weight/z_crit leaks in the substrate E²(z) that feeds the growth ODE (WIP-18); G_eff (the coupling route) is unchanged. """ pk = self._params_tuple() p = self.p if sc is None: sc = self.solve_self_consistency() return compute_S8( pk['B'], pk['zs'], pk['C'], pk['zb'], pk['gamma'], p.eta_crust.value, p.z_dip.value, p.sigma_sup.value, p.mixing.value, p.eta_mix.value, p.sigma_mix.value, pk['Om'], pk['OL'], pk['Z_HARM'], pk['ALPHA_TAIL'], e2_sub_fn=sc['E2'] ) def compute_S8_spline(self) -> float: """S8 for the freeform 15-knot *spline* moraine model (vs the DSW envelope). The DSW Profile-C crust and the spline are two *different fits* to the same DESI BAO + Jia H₀(z) data, so they imply slightly different S8 values. The ONLY difference S8 sees is the BACKGROUND route — the substrate E²(z) that sets the Hubble friction in the growth ODE. The COUPLING route (G_eff at z_crit / z_dip, anchored to the Weinberg angle) is profile-independent and is applied identically to both. The DSW headline value uses the full self-consistent solver (`compute_S8_sc`, closing L1+L2+L3 → 0.8156). That solver is built on the DSW soliton/harmonic-edge envelope (`e_DSW_zc`) and has no spline analogue, so the spline is run on its own *flatness-normalised* (L1) background — the dominant leak by far (L1 ≈ +0.03; the DSW's extra L2+L3 add only ≈ +0.001). The spline's steep crest/trough forcing nearly cancels at first order but leaves a small *net* suppression, so it lands ≈ 0.008–0.011 *below* the DSW: DSW (self-consistent, L1+L2+L3) : S8 = 0.8156 (the headline) spline (flatness-normalised, L1) : S8 = 0.8053 Both sit just *above* the 0.76–0.79 weak-lensing band. Quote them as two genuine, nearly-equal values — NOT the obsolete bootstrap 0.792/0.793 pair, which predates the WIP-18 self-consistency fix and is the *same* model's old leaked-background number, not a real DSW-vs-spline split. """ pk = self._params_tuple() p = self.p Om, OL = pk['Om'], pk['OL'] # spline crust, flatness-normalised base density (L1): Ω_Λ,base = Ω_Λ,tot / f(0) f0 = f_mod_spline(0.0, pk['C'], pk['zb']) OLb = (1.0 - Om - Or) / f0 e2_spline_L1 = lambda z: (Om*(1+z)**3 + Or*(1+z)**4 + OLb * f_mod_spline(z, pk['C'], pk['zb'])) return compute_S8( pk['B'], pk['zs'], pk['C'], pk['zb'], pk['gamma'], p.eta_crust.value, p.z_dip.value, p.sigma_sup.value, p.mixing.value, p.eta_mix.value, p.sigma_mix.value, Om, OL, pk['Z_HARM'], pk['ALPHA_TAIL'], e2_sub_fn=e2_spline_L1) def compute_S8_bore(self) -> float: """S8 for the 14-knot bore-skeleton spline (flatness-normalised, L1). Same treatment as compute_S8_spline (the two are alternative freeform fits to the same data): the crust's f(z) enters growth ONLY through the BACKGROUND route, run on the L1 flatness-normalised expansion law; the G_eff COUPLING route is profile-independent and applied identically. """ pk = self._params_tuple() p = self.p Om, OL = pk['Om'], pk['OL'] # bore crust, flatness-normalised base density (L1): Ω_Λ,base = Ω_Λ,tot / f(0) f0 = f_mod_bore(0.0, pk['C'], pk['zb']) OLb = (1.0 - Om - Or) / f0 e2_bore_L1 = lambda z: (Om*(1+z)**3 + Or*(1+z)**4 + OLb * f_mod_bore(z, pk['C'], pk['zb'])) return compute_S8( pk['B'], pk['zs'], pk['C'], pk['zb'], pk['gamma'], p.eta_crust.value, p.z_dip.value, p.sigma_sup.value, p.mixing.value, p.eta_mix.value, p.sigma_mix.value, Om, OL, pk['Z_HARM'], pk['ALPHA_TAIL'], e2_sub_fn=e2_bore_L1) def compute_S8_decomposition(self) -> dict: """Decompose S8 into the background and coupling routes (WIP-18/WIP-20). Separates the two channels and the bootstrap→self-consistent shift, so the chapter's intermediate numbers (coupling-only 0.783, L1-only 0.814, self-consistent net 0.816) are reproducible. Returns a dict of S8s. """ pk = self._params_tuple() p = self.p Om, OL = pk['Om'], pk['OL'] def s8(e2_fn, geff_on): return compute_S8( pk['B'], pk['zs'], pk['C'], pk['zb'], pk['gamma'], (p.eta_crust.value if geff_on else 0.0), p.z_dip.value, p.sigma_sup.value, (p.mixing.value and geff_on), (p.eta_mix.value if geff_on else 0.0), p.sigma_mix.value, Om, OL, pk['Z_HARM'], pk['ALPHA_TAIL'], e2_sub_fn=e2_fn) e2_flat = lambda z: E2_lcdm(z, Om, OL) # f=1 background e2_boot = lambda z: E2_sub(z, pk['B'], pk['zs'], pk['C'], pk['zb'], pk['gamma'], Om, OL, pk['Z_HARM'], pk['ALPHA_TAIL']) # bootstrap (L1 open) sc = self.solve_self_consistency() f0 = f_mod(0.0, pk['B'], pk['zs'], pk['C'], pk['zb'], pk['gamma'], Om, OL, pk['Z_HARM'], pk['ALPHA_TAIL'], True) OLb = (1.0 - Om - Or) / f0 e2_L1 = lambda z: (Om*(1+z)**3 + Or*(1+z)**4 + OLb*f_mod( z, pk['B'], pk['zs'], pk['C'], pk['zb'], pk['gamma'], Om, OL, pk['Z_HARM'], pk['ALPHA_TAIL'], True)) # flatness only return dict( lcdm = s8(e2_flat, False), # 0.831 coupling_flat = s8(e2_flat, True), # 0.783 background_sc = s8(sc['E2'], False), # background lift bootstrap = s8(e2_boot, True), # 0.792 (leak open) flatness_only = s8(e2_L1, True), # 0.814 (L1 only) self_consistent = s8(sc['E2'], True), # 0.816 (L1+L2+L3) OL_base = OLb, z_crit = sc['z_crit']) def compute_CPL_fit(self) -> dict: """Best-fit CPL (w0, wa) and DESI tension χ².""" pk = self._params_tuple() cpl = fit_CPL(pk['B'], pk['zs'], pk['C'], pk['zb'], pk['gamma'], pk['Om'], pk['OL'], pk['Z_HARM'], pk['ALPHA_TAIL']) return {**cpl, 'chi_sq_desi': chi_sq_cpl(cpl)} def compute_total_chi2(self) -> dict: """Total χ² breakdown (BAO + Jia) for DSW, spline, bore, and ΛCDM. Returns a dict with keys 'dsw', 'spline', 'bore', 'lcdm', each containing {'chi2_bao', 'chi2_jia', 'chi2_total'}. """ pk = self._params_tuple() dsw = total_chi2_dsw( pk['B'], pk['zs'], pk['C'], pk['zb'], pk['gamma'], pk['Om'], pk['OL'], pk['Z_HARM'], pk['ALPHA_TAIL'], self.H0) spline = total_chi2_spline( pk['C'], pk['zb'], pk['Om'], pk['OL'], self.H0) bore = total_chi2_bore( pk['C'], pk['zb'], pk['Om'], pk['OL'], self.H0) lcdm = total_chi2_lcdm(OM_LCDM, OL_LCDM, self.H0) return {'dsw': dsw, 'spline': spline, 'bore': bore, 'lcdm': lcdm} def summary(self, steps: int = 200) -> dict: """Run all diagnostics and return a single summary dict.""" return { 'preset' : 'substrate' if self.p.dsw.value else 'lcdm', 'H0' : self.H0, 'Om' : self.Om, 'S8' : self.compute_S8(), 'CPL' : self.compute_CPL_fit(), 'H0_running' : self.compute_H0_running(steps), 'w_eff' : self.compute_w_eff(steps), 'energy_density' : self.compute_energy_density(steps), 'total_chi2' : self.compute_total_chi2(), } # ── Entry point ─────────────────────────────────────────────── # ── Spline-variant H0 inference (mirrors H0inferredSpline in JS) ── def H0_inferred_spline(z: float, C: float, zb: float, Om: float, OL: float, H0_BASE: float) -> float: """H0(z) inferred by an ΛCDM observer against the spline moraine model.""" if z < 0.005: z = 0.005 I_sub = simp_int(z, lambda zz: E2_spline(zz, C, zb, Om, OL), 300) I_lcdm = simp_int(z, E2_lcdm_std, 80) if I_sub < 1e-10: return H0_BASE return H0_BASE * I_lcdm / I_sub def H0_inferred_bore(z: float, C: float, zb: float, Om: float, OL: float, H0_BASE: float) -> float: """H0(z) inferred by an ΛCDM observer against the bore-skeleton model.""" if z < 0.005: z = 0.005 I_sub = simp_int(z, lambda zz: E2_bore(zz, C, zb, Om, OL), 300) I_lcdm = simp_int(z, E2_lcdm_std, 80) if I_sub < 1e-10: return H0_BASE return H0_BASE * I_lcdm / I_sub # ── Curve-file import (direct-sample mode) ─────────────────── def _load_curve_csv(filepath: str) -> dict: """ Read a bore-desi-jia fit-curve CSV and return sorted arrays for linear interpolation. Expected columns: z, f_substrate_base, f_crust_spline, f_substrate_spline, f_crust_dsw, f_substrate_dsw Returns ------- dict with keys 'z', 'f_dsw', 'f_spline' — each a list[float]. """ import csv zs, f_dsw, f_spl = [], [], [] with open(filepath, newline='') as fh: reader = csv.DictReader(fh) for row in reader: zs.append(float(row['z'])) f_dsw.append(float(row['f_substrate_dsw'])) f_spl.append(float(row['f_substrate_spline'])) return {'z': zs, 'f_dsw': f_dsw, 'f_spline': f_spl} def _lerp_curve(z: float, zs: list, fs: list) -> float: """Linearly interpolate f(z) from sorted parallel arrays zs, fs.""" if z <= zs[0]: return fs[0] if z >= zs[-1]: return fs[-1] # Binary search for bracket lo, hi = 0, len(zs) - 1 while hi - lo > 1: mid = (lo + hi) // 2 if zs[mid] <= z: lo = mid else: hi = mid t = (z - zs[lo]) / (zs[hi] - zs[lo]) if zs[hi] != zs[lo] else 0.0 return fs[lo] + t * (fs[hi] - fs[lo]) def E2_curve(z: float, Om: float, OL: float, zs: list, fs: list) -> float: """E²(z) using a pre-computed f(z) curve (linearly interpolated).""" return Om * (1+z)**3 + Or * (1+z)**4 + OL * _lerp_curve(z, zs, fs) def H0_inferred_curve(z: float, Om: float, OL: float, H0_BASE: float, zs: list, fs: list) -> float: """ H₀(z) inferred by an ΛCDM observer against a model whose f(z) is provided as a pre-computed curve (linearly interpolated). """ if z < 0.005: z = 0.005 I_model = simp_int(z, lambda zz: E2_curve(zz, Om, OL, zs, fs), 300) I_lcdm = simp_int(z, E2_lcdm_std, 80) if I_model < 1e-10: return H0_BASE return H0_BASE * I_lcdm / I_model def export_from_curve(curve_filepath: str, filepath: str = 'substrate_curve_data.csv', H0_BASE: float = 67.4, Om: float = 0.315) -> None: """ Read a bore-desi-jia fit-curve CSV and compute H₀(z) for each row. For every redshift z in the input CSV, this computes: H0_dsw — inferred H₀ using f_substrate_dsw curve H0_spline — inferred H₀ using f_substrate_spline curve H0_lcdm — inferred H₀ assuming pure ΛCDM (f=1) Output CSV columns: z, H0_dsw, H0_spline, H0_lcdm, f_dsw, f_spline, f_lcdm """ OL = 1.0 - Om - Or curve = _load_curve_csv(curve_filepath) zs_arr = curve['z'] fs_dsw = curve['f_dsw'] fs_spline = curve['f_spline'] header = ['z', 'H0_dsw', 'H0_spline', 'H0_lcdm', 'f_dsw', 'f_spline', 'f_lcdm'] rows = [] for i, z in enumerate(zs_arr): if z < 0.005: # Skip z=0 — H0_inferred is undefined there continue h0_dsw = H0_inferred_curve(z, Om, OL, H0_BASE, zs_arr, fs_dsw) h0_spl = H0_inferred_curve(z, Om, OL, H0_BASE, zs_arr, fs_spline) h0_lcdm = H0_inferred_curve(z, OM_LCDM, OL_LCDM, H0_BASE, zs_arr, [1.0]*len(zs_arr)) rows.append([z, h0_dsw, h0_spl, h0_lcdm, fs_dsw[i], fs_spline[i], 1.0]) import csv with open(filepath, 'w', newline='') as fh: w = csv.writer(fh) w.writerow(header) for row in rows: w.writerow([f'{v:.6f}' for v in row]) print(f"Wrote {len(rows)} rows × {len(header)} columns → {filepath}") print(f" z range : {rows[0][0]:.4f} – {rows[-1][0]:.4f}") print(f" source : {curve_filepath}") print(f" H0_BASE : {H0_BASE}") print(f" Om : {Om}") # ── CSV export ──────────────────────────────────────────────── def export_data(steps: int, filepath: str = 'substrate_data.csv') -> None: """ Compute and export a multi-column CSV over z ∈ [0.005, 3.0]. Columns ------- z redshift H0_dsw H0 inferred by ΛCDM observer vs substrate/DSW model H0_spline H0 inferred by ΛCDM observer vs 15-knot spline moraine model H0_bore H0 inferred by ΛCDM observer vs 14-knot bore-skeleton model H0_lcdm H0 inferred vs pure ΛCDM (flat reference, ≈ H0_BASE) f_dsw energy density modifier f_mod(z) — substrate model f_spline energy density modifier f_mod_spline(z) — 15-knot moraine f_bore energy density modifier f_mod_bore(z) — 14-knot bore skeleton f_lcdm energy density modifier — ΛCDM (always 1.0) """ sub = createSubstrate() lcdm = createLCDM() sp = sub # spline shares substrate C, zb, Om, OL # Unpack the parameters we need B = sub.B.value zs = sub.z_s.value C = sub.C.value zb = sub.z_b.value gamma = sub.gamma.value Om = sub.Om.value OL = 1.0 - Om - Or Z_HARM = sub.z_harm.value ALPHA = sub.alpha.value H0 = sub.H0.value z_axis = [max(0.005, 3.0 * i / steps) for i in range(1, steps + 1)] header = ['z', 'H0_dsw', 'H0_spline', 'H0_bore', 'H0_lcdm', 'f_dsw', 'f_spline', 'f_bore', 'f_lcdm'] rows = [] for z in z_axis: h0_dsw = H0_inferred(z, B, zs, C, zb, gamma, Om, OL, Z_HARM, ALPHA, H0) h0_spl = H0_inferred_spline(z, C, zb, Om, OL, H0) h0_bore = H0_inferred_bore(z, C, zb, Om, OL, H0) h0_lcdm = H0_inferred(z, 0.0, zs, 0.0, zb, gamma, OM_LCDM, OL_LCDM, 0.0, ALPHA, H0) fd_dsw = f_mod(z, B, zs, C, zb, gamma, Om, OL, Z_HARM, ALPHA, dsw=True) fd_spl = f_mod_spline(z, C, zb) fd_bore = f_mod_bore(z, C, zb) fd_lcdm = 1.0 rows.append([z, h0_dsw, h0_spl, h0_bore, h0_lcdm, fd_dsw, fd_spl, fd_bore, fd_lcdm]) import csv with open(filepath, 'w', newline='') as fh: w = csv.writer(fh) w.writerow(header) for row in rows: w.writerow([f'{v:.6f}' for v in row]) print(f"Wrote {len(rows)} rows × {len(header)} columns → {filepath}") print(f" z range : {z_axis[0]:.4f} – {z_axis[-1]:.4f}") print(f" steps : {steps}") def export_bao(filepath: str) -> None: """ Export BAO residuals for each DESI_BAO observation point. Columns ------- z observation redshift type DM | DH | DV label survey label (BGS, LRG1, …) observed measured value in units of r_d sigma 1-sigma uncertainty pred_dsw substrate/DSW prediction pred_spline spline moraine prediction pred_lcdm flat ΛCDM prediction err_dsw (pred_dsw - observed) / sigma [pull] err_spline (pred_spline - observed) / sigma err_lcdm (pred_lcdm - observed) / sigma """ sub = createSubstrate() B = sub.B.value zs = sub.z_s.value C = sub.C.value zb = sub.z_b.value gamma = sub.gamma.value Om = sub.Om.value OL = 1.0 - Om - Or Z_HARM = sub.z_harm.value ALPHA = sub.alpha.value H0 = sub.H0.value preds_dsw = bao_predict_dsw(B, zs, C, zb, gamma, Om, OL, Z_HARM, ALPHA, H0) preds_spline = bao_predict_spline(C, zb, Om, OL, H0) preds_bore = bao_predict_bore(C, zb, Om, OL, H0) preds_lcdm = bao_predict_lcdm(OM_LCDM, OL_LCDM, H0) header = ['z', 'type', 'label', 'observed', 'sigma', 'pred_dsw', 'pred_spline', 'pred_bore', 'pred_lcdm', 'err_dsw', 'err_spline', 'err_bore', 'err_lcdm'] import csv with open(filepath, 'w', newline='') as fh: w = csv.writer(fh) w.writerow(header) for k, d in enumerate(DESI_BAO): obs = d['val'] sig = d['sig'] pd = preds_dsw[k] ps = preds_spline[k] pb = preds_bore[k] pl = preds_lcdm[k] w.writerow([ f"{d['z']:.3f}", d['type'], d['label'], f"{obs:.4f}", f"{sig:.4f}", f"{pd:.4f}", f"{ps:.4f}", f"{pb:.4f}", f"{pl:.4f}", f"{(pd - obs) / sig:+.4f}", f"{(ps - obs) / sig:+.4f}", f"{(pb - obs) / sig:+.4f}", f"{(pl - obs) / sig:+.4f}", ]) bao_chi2_dsw_val = sum(((preds_dsw[k] - d['val']) / d['sig'])**2 for k, d in enumerate(DESI_BAO)) bao_chi2_spline_val = sum(((preds_spline[k] - d['val']) / d['sig'])**2 for k, d in enumerate(DESI_BAO)) bao_chi2_bore_val = sum(((preds_bore[k] - d['val']) / d['sig'])**2 for k, d in enumerate(DESI_BAO)) bao_chi2_lcdm_val = sum(((preds_lcdm[k] - d['val']) / d['sig'])**2 for k, d in enumerate(DESI_BAO)) print(f"Wrote {len(DESI_BAO)} BAO points → {filepath}") print(f" BAO χ² DSW : {bao_chi2_dsw_val:.3f}") print(f" BAO χ² spline : {bao_chi2_spline_val:.3f}") print(f" BAO χ² bore : {bao_chi2_bore_val:.3f}") print(f" BAO χ² ΛCDM : {bao_chi2_lcdm_val:.3f}") # ── CLI entry point ─────────────────────────────────────────── if __name__ == '__main__': import argparse parser = argparse.ArgumentParser( description='Substrate Galactic — dark energy / H0 tension model') parser.add_argument( '--output-data', metavar='STEPS', type=int, default=None, help='Export CSV with H0 and energy density curves at STEPS redshift points') parser.add_argument( '--output-file', metavar='FILE', type=str, default='substrate_data.csv', help='Destination CSV filename (default: substrate_data.csv)') parser.add_argument( '--bao-file', metavar='FILE', type=str, default=None, help='Export BAO residuals (DSW / spline / ΛCDM pulls) to FILE') parser.add_argument( '--curve-file', metavar='FILE', type=str, default=None, help='Input a fit-curve CSV (z, f_substrate_spline, f_substrate_dsw) ' 'and compute H0 values from the direct samples') parser.add_argument( '--curve-output', metavar='FILE', type=str, default='substrate_curve_data.csv', help='Destination CSV for --curve-file output (default: substrate_curve_data.csv)') parser.add_argument( '--s8-decompose', action='store_true', help='Print the S8 background/coupling route decomposition (WIP-18/WIP-20)') args = parser.parse_args() if args.curve_file is not None: export_from_curve(curve_filepath=args.curve_file, filepath=args.curve_output) if args.output_data is not None: export_data(steps=args.output_data, filepath=args.output_file) if args.bao_file is not None: export_bao(filepath=args.bao_file) if not any([args.curve_file, args.output_data, args.bao_file]): # Default: print the summary diagnostics as before model = SubstrateModel(createLCDM()) curve = model.compute_energy_density(steps=200) print(f"\n{'='*48}") print(" ΛCDM") print(f"{'='*48}") #print(f" f(z=0) : {curve.f_at_zero:.4f}") #print(f" f_mod range : [{min(curve.f_mod):.4f}, {max(curve.f_mod):.4f}]") #print(f" spline knots : {len(curve.knot_z)}") #print(f" S8 (bootstrap) : {model.compute_S8():.4f}") h0r = model.compute_H0_running(steps=50) print(f" H0 running χ² : {h0r['chi_sq']:.3f}") cpl = model.compute_CPL_fit() print(f" CPL (w0, wa) : ({cpl['w0']}, {cpl['wa']})") print(f" DESI tension χ² : {cpl['chi_sq_desi']:.3f}") crossing = model.compute_w_eff(steps=200)['phantom_crossing'] print(f" Phantom crossing : z = {crossing}") # Total χ² breakdown (BAO + Jia) — matches # fit_freeform_spline.py methodology chi2 = model.compute_total_chi2() print(f"\n {'─'*44}") print(f" Total χ² = BAO + Jia (diagonal covariance)") print(f" {'─'*44}") print(f" {'Model':<12s} {'BAO':>8s} {'Jia':>8s} {'Total':>8s}") print(f" {'─'*12} {'─'*8} {'─'*8} {'─'*8}") d = chi2['lcdm'] name = 'ΛCDM' print(f" {name:<12s} {d['chi2_bao']:8.3f}" f" {d['chi2_jia']:8.3f}" f" {d['chi2_total']:8.3f}") params = createSubstrate() model = SubstrateModel(params) curve = model.compute_energy_density(steps=200) print(f"\n{'='*48}") print(f" Substrate") print(f"{'='*48}") print(f" f(z=0) : {curve.f_at_zero:.4f}") print(f" f_mod range : [{min(curve.f_mod):.4f}, {max(curve.f_mod):.4f}]") print(f" spline knots : {len(curve.knot_z)}") #print(f" S8 (bootstrap) : {model.compute_S8():.4f}") if params.dsw.value and params.B.value > 0: # Full Friedmann self-consistency (WIP-18): close L1/L2/L3 in the # background E²(z) that drives growth. G_eff (coupling route) # stays out of Friedmann, so it is untouched here. sc = model.solve_self_consistency() print(f" S8 (self-consist) : {model.compute_S8_sc(sc):.4f}") print(f" self-con f(0) : {sc['f0']:.4f} Ω_Λ,base = {sc['OL_base']:.4f}" f" z_crit = {sc['z_crit']:.4f} E²(0) = {sc['E2_0']:.4f}") h0r = model.compute_H0_running(steps=50) print(f" H0 running χ² : {h0r['chi_sq']:.3f}") cpl = model.compute_CPL_fit() print(f" CPL (w0, wa) : ({cpl['w0']}, {cpl['wa']})") print(f" DESI tension χ² : {cpl['chi_sq_desi']:.3f}") crossing = model.compute_w_eff(steps=200)['phantom_crossing'] print(f" Phantom crossing : z = {crossing}") # Total χ² breakdown (BAO + Jia) — matches # fit_freeform_spline.py methodology chi2 = model.compute_total_chi2() print(f"\n {'─'*44}") print(f" Total χ² = BAO + Jia (diagonal covariance)") print(f" {'─'*44}") print(f" {'Model':<12s} {'BAO':>8s} {'Jia':>8s} {'Total':>8s}") print(f" {'─'*12} {'─'*8} {'─'*8} {'─'*8}") for tag in ('dsw', 'spline', 'bore', 'lcdm'): d = chi2[tag] name = {'dsw': 'DSW', 'spline': 'Spline (15k)', 'bore': 'Bore (14k)', 'lcdm': 'ΛCDM'}[tag] print(f" {name:<12s} {d['chi2_bao']:8.3f}" f" {d['chi2_jia']:8.3f}" f" {d['chi2_total']:8.3f}") m = SubstrateModel(createSubstrate()) d = m.compute_S8_decomposition() print(f"\n{'='*56}") print(" S8 decomposition — background × coupling routes") print(f" (full Friedmann self-consistency, WIP-18/WIP-20)") print(f"{'='*56}") print(f" ΛCDM (no crust, no G_eff) S8 = {d['lcdm']:.4f}") #print(f" coupling route only (G_eff, flat bg) S8 = {d['coupling_flat']:.4f}") #print(f" background route only (SC E², no G_eff) S8 = {d['background_sc']:.4f}") #print(f" {'-'*52}") #print(f" bootstrap (G_eff + leaked bg) S8 = {d['bootstrap']:.4f}") #print(f" flatness only (G_eff + L1 bg) S8 = {d['flatness_only']:.4f}") #print(f" self-consist (G_eff + L1+L2+L3 bg) S8 = {d['self_consistent']:.4f}") print(f" {'-'*52}") print(f" The substrate models:") print(f" DSW envelope (self-consistent) S8 = {m.compute_S8_sc():.4f}") print(f" spline (15-knot) (flatness-norm L1) S8 = {m.compute_S8_spline():.4f}") print(f" bore (14-knot) (flatness-norm L1) S8 = {m.compute_S8_bore():.4f}") print(f" {'-'*52}") print(f" Ω_Λ,base = {d['OL_base']:.4f} z_crit = {d['z_crit']:.4f}")