"""
nR cosmology vs. LambdaCDM: strict low-redshift SN + BAO + optional local H0 prior.
Free rd, full diagnostics, production-ready.
Required files:
Pantheon+SH0ES.dat
Pantheon+SH0ES_STAT+SYS.cov
boss_dr12_cov.npy
"""
from pathlib import Path
import numpy as np
import pandas as pd
from scipy.integrate import cumulative_trapezoid
from scipy.linalg import solve_triangular
from scipy.optimize import differential_evolution, minimize
# ----------------------- Configuration -----------------------
SN_FILE = Path("Pantheon+SH0ES.dat")
SN_COV_FILE = Path("Pantheon+SH0ES_STAT+SYS.cov")
BAO_COV_FILE = Path("boss_dr12_cov.npy")
C_KM_S = 299792.458
TCMB = 2.7255
NEFF = 3.046
# BAO consensus (BOSS DR12)
BAO_Z = np.array([0.38, 0.51, 0.61])
BAO_OBS = np.array([
1512.0, 81.2,
1975.0, 90.9,
2307.0, 99.0
], dtype=float)
RD_FID_MPC = 147.78 # only used for scaling convention
# Analysis switches
USE_LOCAL_H0_PRIOR = True
H0_LOC = 73.04
SIG_H0 = 1.04
FREE_RD = True # if False, rd is fixed to RD_FIXED
RD_FIXED = 147.78
# Grid resolution for distance integrals
N_GRID = 25000
# ----------------------- File validation -----------------------
def load_pantheon_covariance(path):
if not path.exists():
raise FileNotFoundError(f"Missing covariance file: {path}")
with path.open("r", encoding="utf-8") as f:
first_line = f.readline().strip()
try:
n = int(first_line)
except ValueError as exc:
raise ValueError(
"Pantheon covariance first line must contain matrix size N."
) from exc
raw = np.loadtxt(path, skiprows=1, dtype=float).ravel()
if raw.size != n * n:
raise ValueError(
f"Covariance declares {n}x{n}, but contains {raw.size} values."
)
cov = raw.reshape(n, n)
if not np.all(np.isfinite(cov)):
raise ValueError("Pantheon covariance contains NaN or infinity.")
asym = np.max(np.abs(cov - cov.T))
scale = max(float(np.max(np.abs(cov))), 1.0)
if asym > 1e-6 * scale:
raise ValueError(f"Covariance is not symmetric enough: {asym:g}")
return 0.5 * (cov + cov.T)
def validate_covariance(cov, name):
cov = np.asarray(cov, dtype=float)
if cov.ndim != 2 or cov.shape[0] != cov.shape[1]:
raise ValueError(f"{name}: covariance must be square.")
if not np.all(np.isfinite(cov)):
raise ValueError(f"{name}: covariance contains non-finite values.")
scale = max(float(np.max(np.abs(cov))), 1.0)
if np.max(np.abs(cov - cov.T)) > 1e-10 * scale:
raise ValueError(f"{name}: covariance is not symmetric.")
try:
L = np.linalg.cholesky(cov)
except Exception as exc:
raise ValueError(f"{name}: covariance is not positive definite.") from exc
return L
def load_data():
if not SN_FILE.exists():
raise FileNotFoundError(f"Missing Pantheon+ data: {SN_FILE}")
if not BAO_COV_FILE.exists():
raise FileNotFoundError(f"Missing official BAO covariance: {BAO_COV_FILE}")
sn = pd.read_csv(SN_FILE, sep=r"\s+", comment="#")
required = {"zHD", "m_b_corr", "IS_CALIBRATOR"}
missing = required - set(sn.columns)
if missing:
raise ValueError(f"Pantheon+ file is missing columns: {sorted(missing)}")
C_all = load_pantheon_covariance(SN_COV_FILE)
if len(sn) != C_all.shape[0]:
raise ValueError(
f"Pantheon rows ({len(sn)}) do not match covariance size ({C_all.shape[0]})."
)
# Exclude calibrators; absolute scale comes from BAO + optional local H0 prior
is_cal = pd.to_numeric(sn["IS_CALIBRATOR"], errors="raise").to_numpy()
keep = (is_cal == 0)
if not np.any(keep):
raise ValueError("No non-calibrator supernovae selected.")
indices = np.flatnonzero(keep)
sn = sn.iloc[indices].copy().reset_index(drop=True)
C_sn = C_all[np.ix_(indices, indices)]
z_sn = pd.to_numeric(sn["zHD"], errors="raise").to_numpy(float)
m_sn = pd.to_numeric(sn["m_b_corr"], errors="raise").to_numpy(float)
if not np.all(np.isfinite(z_sn)) or not np.all(np.isfinite(m_sn)):
raise ValueError("SN data contain non-finite values.")
if np.any(z_sn <= 0):
raise ValueError("SN redshifts must be strictly positive.")
C_bao = np.load(BAO_COV_FILE, allow_pickle=False)
if C_bao.shape != (6, 6):
raise ValueError(f"BAO covariance must be 6x6; received {C_bao.shape}.")
L_sn = validate_covariance(C_sn, "Pantheon+ covariance")
L_bao = validate_covariance(C_bao, "BOSS DR12 covariance")
return z_sn, m_sn, L_sn, L_bao
# ----------------------- Cosmological models -----------------------
def omega_radiation(H0):
h = H0 / 100.0
omega_gamma = 2.469e-5 * (TCMB / 2.7255) ** 4
return omega_gamma * (1.0 + 0.22710731766 * NEFF) / h**2
def E2_nr(z, H0, Om, xi, alpha, beta):
z = np.asarray(z, dtype=float)
Or = omega_radiation(H0)
Ode = 1.0 - Om - Or
zp1 = 1.0 + z
nR = xi * zp1**alpha * (zp1**beta - 1.0)
return Om * zp1**3 + Or * zp1**4 + Ode + nR
def E2_lcdm(z, H0, Om):
z = np.asarray(z, dtype=float)
Or = omega_radiation(H0)
Ode = 1.0 - Om - Or
zp1 = 1.0 + z
return Om * zp1**3 + Or * zp1**4 + Ode
def distance_grid_generic(z_eval, H0, e2_func, e2_args, n_grid=N_GRID):
z_eval = np.asarray(z_eval, dtype=float)
if np.any(z_eval < 0) or not np.all(np.isfinite(z_eval)):
raise ValueError("Evaluation redshifts must be finite and non-negative.")
zmax = float(np.max(z_eval))
if zmax == 0.0:
return np.zeros_like(z_eval)
zg = np.linspace(0.0, zmax, n_grid)
e2 = e2_func(zg, H0, *e2_args)
if not np.all(np.isfinite(e2)) or np.any(e2 <= 0):
raise ValueError("Non-physical model: E(z)^2 is non-positive or non-finite.")
invE = 1.0 / np.sqrt(e2)
integral = cumulative_trapezoid(invE, zg, initial=0.0)
dc_grid = (C_KM_S / H0) * integral
return np.interp(z_eval, zg, dc_grid)
def predict_bao_generic(H0, e2_func, e2_args, rd_mpc):
dc = distance_grid_generic(BAO_Z, H0, e2_func, e2_args)
scale = RD_FID_MPC / rd_mpc
dm_scaled = dc * scale
e2 = e2_func(BAO_Z, H0, *e2_args)
H = H0 * np.sqrt(e2)
H_scaled = H / scale # H * (rd / rd_fid)
pred = np.empty(6, dtype=float)
pred[0::2] = dm_scaled
pred[1::2] = H_scaled
return pred
# ----------------------- Likelihood core -----------------------
def chi2_sn_generic(z_sn, m_obs, L_sn, H0, e2_func, e2_args):
dc = distance_grid_generic(z_sn, H0, e2_func, e2_args)
dl = (1.0 + z_sn) * dc
if np.any(dl <= 0) or not np.all(np.isfinite(dl)):
raise ValueError("Non-positive or non-finite luminosity distance.")
mu = 5.0 * np.log10(dl) + 25.0
residual0 = m_obs - mu
y = solve_triangular(L_sn, residual0, lower=True)
u = solve_triangular(L_sn, np.ones_like(residual0), lower=True)
denom = float(u @ u)
if not np.isfinite(denom) or denom <= 0:
raise ValueError("Invalid profiled nuisance-parameter denominator.")
M_hat = float(u @ y / denom)
whitened = y - M_hat * u
chi2 = float(whitened @ whitened)
return chi2, M_hat
def chi2_bao_generic(L_bao, H0, e2_func, e2_args, rd_mpc):
pred = predict_bao_generic(H0, e2_func, e2_args, rd_mpc)
residual = BAO_OBS - pred
whitened = solve_triangular(L_bao, residual, lower=True)
return float(whitened @ whitened)
# ----------------------- Objectives & Bounds -----------------------
def get_bounds_nr():
bounds = [
(55.0, 85.0), # H0
(0.15, 0.45), # Om
(-1.5, 1.5), # xi
(-4.0, 4.0), # alpha
(-4.0, 4.0), # beta
]
if FREE_RD:
bounds.append((130.0, 160.0)) # rd
return bounds
def get_bounds_lcdm():
bounds = [
(55.0, 85.0), # H0
(0.15, 0.45), # Om
]
if FREE_RD:
bounds.append((130.0, 160.0)) # rd
return bounds
def objective_nr(theta, z_sn, m_sn, L_sn, L_bao):
theta = np.asarray(theta, dtype=float)
if FREE_RD:
H0, Om, xi, alpha, beta, rd = theta
else:
H0, Om, xi, alpha, beta = theta
rd = RD_FIXED
try:
zcheck = np.linspace(0.0, max(float(np.max(z_sn)), float(np.max(BAO_Z))), 3000)
e2 = E2_nr(zcheck, H0, Om, xi, alpha, beta)
if not np.all(np.isfinite(e2)) or np.any(e2 <= 0):
return 1e100
sn_val, _ = chi2_sn_generic(z_sn, m_sn, L_sn, H0, E2_nr, (Om, xi, alpha, beta))
bao_val = chi2_bao_generic(L_bao, H0, E2_nr, (Om, xi, alpha, beta), rd)
total = sn_val + bao_val
if USE_LOCAL_H0_PRIOR:
total += ((H0 - H0_LOC) / SIG_H0) ** 2
return total if np.isfinite(total) else 1e100
except (ValueError, FloatingPointError, np.linalg.LinAlgError):
return 1e100
def objective_lcdm(theta, z_sn, m_sn, L_sn, L_bao):
theta = np.asarray(theta, dtype=float)
if FREE_RD:
H0, Om, rd = theta
else:
H0, Om = theta
rd = RD_FIXED
try:
zcheck = np.linspace(0.0, max(float(np.max(z_sn)), float(np.max(BAO_Z))), 3000)
e2 = E2_lcdm(zcheck, H0, Om)
if not np.all(np.isfinite(e2)) or np.any(e2 <= 0):
return 1e100
sn_val, _ = chi2_sn_generic(z_sn, m_sn, L_sn, H0, E2_lcdm, (Om,))
bao_val = chi2_bao_generic(L_bao, H0, E2_lcdm, (Om,), rd)
total = sn_val + bao_val
if USE_LOCAL_H0_PRIOR:
total += ((H0 - H0_LOC) / SIG_H0) ** 2
return total if np.isfinite(total) else 1e100
except (ValueError, FloatingPointError, np.linalg.LinAlgError):
return 1e100
# ----------------------- Diagnostics helper -----------------------
def evaluate_model(theta, z_sn, m_sn, L_sn, L_bao, is_nr=True):
"""Return detailed chi2 pieces and M for a given parameter vector."""
theta = np.asarray(theta, dtype=float)
if is_nr:
if FREE_RD:
H0, Om, xi, alpha, beta, rd = theta
else:
H0, Om, xi, alpha, beta = theta
rd = RD_FIXED
e2_func = E2_nr
e2_args = (Om, xi, alpha, beta)
else:
if FREE_RD:
H0, Om, rd = theta
else:
H0, Om = theta
rd = RD_FIXED
e2_func = E2_lcdm
e2_args = (Om,)
sn_val, M_hat = chi2_sn_generic(z_sn, m_sn, L_sn, H0, e2_func, e2_args)
bao_val = chi2_bao_generic(L_bao, H0, e2_func, e2_args, rd)
h0_val = ((H0 - H0_LOC) / SIG_H0) ** 2 if USE_LOCAL_H0_PRIOR else 0.0
total = sn_val + bao_val + h0_val
return {
"H0": H0,
"Om": Om,
"rd": rd,
"chi2_sn": sn_val,
"chi2_bao": bao_val,
"chi2_H0": h0_val,
"chi2_total": total,
"M": M_hat,
"params": theta,
}
# ----------------------- Execution Runner -----------------------
def run_pipeline():
z_sn, m_sn, L_sn, L_bao = load_data()
bounds_nr = get_bounds_nr()
bounds_lcdm = get_bounds_lcdm()
print("Configuration")
print(f" USE_LOCAL_H0_PRIOR = {USE_LOCAL_H0_PRIOR} ({H0_LOC} ± {SIG_H0})")
print(f" FREE_RD = {FREE_RD}")
print(f" SN sample = non-calibrators only (M profiled)")
print()
print("--- Fitting nR Model ---")
de_nr = differential_evolution(
objective_nr, bounds=bounds_nr, args=(z_sn, m_sn, L_sn, L_bao),
seed=20261012, strategy="best1bin", maxiter=900, popsize=18,
tol=1e-7, polish=False, workers=1
)
fit_nr = minimize(
objective_nr, x0=de_nr.x, args=(z_sn, m_sn, L_sn, L_bao),
method="L-BFGS-B", bounds=bounds_nr,
options={"maxiter": 8000, "ftol": 1e-12}
)
print("--- Fitting LambdaCDM Model ---")
de_lcdm = differential_evolution(
objective_lcdm, bounds=bounds_lcdm, args=(z_sn, m_sn, L_sn, L_bao),
seed=20261012, strategy="best1bin", maxiter=900, popsize=18,
tol=1e-7, polish=False, workers=1
)
fit_lcdm = minimize(
objective_lcdm, x0=de_lcdm.x, args=(z_sn, m_sn, L_sn, L_bao),
method="L-BFGS-B", bounds=bounds_lcdm,
options={"maxiter": 8000, "ftol": 1e-12}
)
# Detailed evaluation
res_nr = evaluate_model(fit_nr.x, z_sn, m_sn, L_sn, L_bao, is_nr=True)
res_lcdm = evaluate_model(fit_lcdm.x, z_sn, m_sn, L_sn, L_bao, is_nr=False)
k_nr = len(bounds_nr)
k_lcdm = len(bounds_lcdm)
aic_nr = 2 * k_nr + res_nr["chi2_total"]
aic_lcdm = 2 * k_lcdm + res_lcdm["chi2_total"]
delta_aic = aic_nr - aic_lcdm
print("\n" + "=" * 60)
print("DIRECT MODEL COMPARISON: nR vs. LambdaCDM")
print("=" * 60)
print("\nnR best-fit")
print(f" parameters : {res_nr['params']}")
print(f" H0 : {res_nr['H0']:.3f}")
print(f" Om : {res_nr['Om']:.4f}")
print(f" rd : {res_nr['rd']:.2f} Mpc")
print(f" M (profiled) : {res_nr['M']:.4f}")
print(f" chi2_SN : {res_nr['chi2_sn']:.3f}")
print(f" chi2_BAO : {res_nr['chi2_bao']:.3f}")
print(f" chi2_H0 : {res_nr['chi2_H0']:.3f}")
print(f" chi2_total : {res_nr['chi2_total']:.3f}")
print(f" AIC : {aic_nr:.3f}")
print("\nLambdaCDM best-fit")
print(f" parameters : {res_lcdm['params']}")
print(f" H0 : {res_lcdm['H0']:.3f}")
print(f" Om : {res_lcdm['Om']:.4f}")
print(f" rd : {res_lcdm['rd']:.2f} Mpc")
print(f" M (profiled) : {res_lcdm['M']:.4f}")
print(f" chi2_SN : {res_lcdm['chi2_sn']:.3f}")
print(f" chi2_BAO : {res_lcdm['chi2_bao']:.3f}")
print(f" chi2_H0 : {res_lcdm['chi2_H0']:.3f}")
print(f" chi2_total : {res_lcdm['chi2_total']:.3f}")
print(f" AIC : {aic_lcdm:.3f}")
print("\n" + "-" * 60)
print(f"Delta AIC (nR - LambdaCDM) = {delta_aic:.3f}")
if delta_aic < -2:
print("Conclusion: Data favor nR over LambdaCDM (AIC).")
elif delta_aic > 2:
print("Conclusion: LambdaCDM is favored over nR (AIC).")
else:
print("Conclusion: Statistical tie on AIC.")
print("=" * 60)
return fit_nr, fit_lcdm, res_nr, res_lcdm
if __name__ == "__main__":