substrate_galactic.py
DESI and Hubble cosmology
The dark-energy crust: BAO residuals against DESI, the H_0 ladder, and the S_8 growth suppression from the two Weinberg-anchored G_\mathrm{eff} channels.
Download substrate_galactic.py 1,829 lines · standard library only · run with python3 scripts/substrate_galactic.py
Back to the source-code overview.
scripts/substrate_galactic.py
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}")