#!/usr/bin/env python3 """ T_H0_ENVELOPE/run.py Two-method H_0 envelope test: SH0ES vs Planck ratio against DCT canonical 1/sqrt(P_0). Plus a 17-method scatter overview. The single quantitative test is the SH0ES/Planck ratio agreement with DCT's 1.0830 prediction. """ from __future__ import annotations import json, hashlib, datetime as _dt, math from pathlib import Path import numpy as np import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt HERE = Path(__file__).parent.resolve() P0 = 0.851 RATIO_DCT = 1.0 / math.sqrt(P0) # 1.0830 # Two-method H_PLANCK, S_PLANCK = 67.40, 0.50 # Aghanim+ 2020 H_SHOES, S_SHOES = 73.04, 1.04 # Riess+ 2022 # 17-method compilation (mirrors the exploratory script) COMPILATION = [ ("Planck PR3 CMB", 67.40, 0.50, 0.0, "Aghanim 2020"), ("BAO + BBN inverse", 67.70, 0.80, 0.5, "Macri 2018"), ("DESI BAO + BBN", 67.60, 1.10, 0.5, "DESI 2024"), ("ACT DR4 CMB", 67.90, 1.50, 0.0, "Aiola 2020"), ("TRGB", 69.80, 1.70, 1.5, "Freedman 2021"), ("LMC eclipsing binary", 72.00, 1.30, 2.0, "Pietrzynski 2019"), ("Mira variables", 72.70, 4.60, 2.5, "Huang 2020"), ("Surface brightness fluct.", 73.30, 2.50, 3.0, "Khetan 2021"), ("Type II SCM", 75.80, 5.00, 4.0, "de Jaeger 2020"), ("Tully-Fisher", 75.10, 2.80, 4.5, "Kourkchi 2020"), ("GW170817 standard siren", 70.00, 12.0, 5.0, "LVC 2017"), ("SH0ES Cepheid + SNe", 73.04, 1.04, 6.0, "Riess 2022"), ("JWST Cepheid", 72.60, 2.00, 6.0, "Riess 2024"), ("H0LiCOW lensing", 73.30, 1.80, 7.0, "Wong 2020"), ("TDCOSMO", 67.40, 4.10, 7.0, "Birrer 2020"), ("Pesce megamaser NGC4258", 73.90, 3.00, 8.5, "Pesce 2020"), ("Pesce megamaser UGC3789", 75.40, 3.40, 8.5, "Pesce 2020"), ] def main() -> None: # Two-method ratio R_obs = H_SHOES / H_PLANCK sR = R_obs * math.sqrt((S_SHOES/H_SHOES)**2 + (S_PLANCK/H_PLANCK)**2) sigma_dist = (R_obs - RATIO_DCT) / sR if abs(sigma_dist) <= 1: cls = "PASS" elif abs(sigma_dist) <= 3: cls = "SOFT_PASS" elif abs(sigma_dist) <= 5: cls = "RESIDUAL_TENSION" else: cls = "MISS" # 17-method weighted least squares slope vs rank (post-hoc rank caveat) H = np.array([row[1] for row in COMPILATION]) sH = np.array([row[2] for row in COMPILATION]) rank = np.array([row[3] for row in COMPILATION]) w = 1.0 / sH**2 # weighted linear fit A = np.vstack([rank, np.ones_like(rank)]).T W = np.diag(w) cov = np.linalg.inv(A.T @ W @ A) pars = cov @ A.T @ W @ H slope, intercept = pars s_slope = math.sqrt(cov[0, 0]) s_intercept = math.sqrt(cov[1, 1]) slope_sigma = abs(slope) / s_slope # constant-H null H_const = np.sum(H * w) / np.sum(w) chi2_null = float(np.sum(w * (H - H_const)**2)) chi2_lin = float(np.sum(w * (H - (slope * rank + intercept))**2)) delta_chi2 = chi2_null - chi2_lin # H_E and H_phys envelope from fit H_at_rank0 = intercept H_at_rank8p5 = intercept + 8.5 * slope R_fit = H_at_rank8p5 / H_at_rank0 # Jackknife on slope slopes = [] for i in range(len(H)): keep = np.arange(len(H)) != i Ai = np.vstack([rank[keep], np.ones(keep.sum())]).T Wi = np.diag(w[keep]) ci = np.linalg.inv(Ai.T @ Wi @ Ai) pi_ = ci @ Ai.T @ Wi @ H[keep] slopes.append(pi_[0]) slopes = np.array(slopes) out = { "test_id": "T_H0_ENVELOPE", "ran_at_utc": _dt.datetime.utcnow().isoformat(timespec="seconds") + "Z", "P_0": P0, "DCT_ratio_prediction": RATIO_DCT, "two_method_test": { "H_Planck": H_PLANCK, "sigma_Planck": S_PLANCK, "H_SH0ES": H_SHOES, "sigma_SH0ES": S_SHOES, "ratio_observed": R_obs, "ratio_observed_sigma": sR, "sigma_distance_from_DCT": sigma_dist, "classification": cls, }, "seventeen_method_test": { "n": len(H), "weighted_constant_H_km_s_Mpc": H_const, "chi2_constant": chi2_null, "chi2_linear": chi2_lin, "delta_chi2_constant_minus_linear": delta_chi2, "slope_km_s_Mpc_per_rank": slope, "slope_sigma_from_zero": slope_sigma, "H_at_rank_0_fit": H_at_rank0, "H_at_rank_8p5_fit": H_at_rank8p5, "envelope_ratio_fit": R_fit, "jackknife_slope_min_max": [float(slopes.min()), float(slopes.max())], "jackknife_slope_std": float(slopes.std(ddof=1)), "post_hoc_rank_caveat": ( "Information-density rank assignment is post-hoc; pre-" "registration of the rank is required to promote this from " "post-dictive to pre-dictive. Result is NEUTRAL until that " "happens." ), }, "BAO_coupling_caveat": ( "Homogeneous conformal P(t) cancels from radial null chi(z); legacy " "DESI Delta-chi^2 from rescaling D_M is a retracted map, not a " "homogeneous prediction. Operational H_phys targets are separate " "from photon BAO rulers (DCT-BAO-01; /observables/bao/)." ), } with open(HERE / "recipe.md", "rb") as fh: digest = hashlib.sha256(fh.read()).hexdigest() out["recipe_sha256"] = digest (HERE / "recipe.sha256").write_text(digest + "\n") (HERE / "result.json").write_text(json.dumps(out, indent=2)) md = [ "# T_H0_ENVELOPE — result", "", f"**Run:** {out['ran_at_utc']}", f"**Recipe SHA-256:** `{digest}`", "", "## Two-method ratio test (the operative DCT prediction)", "", f"- DCT canonical ratio `1/√P₀ = 1/√{P0} = {RATIO_DCT:.4f}`", f"- Observed ratio `H_SH0ES / H_Planck = {R_obs:.4f} ± {sR:.4f}`", f"- σ-distance DCT from observed: `{sigma_dist:+.2f}σ`", f"- **Classification: {cls}**", "", "## 17-method scatter test (post-hoc rank — see caveat)", "", f"- Constant-H best fit: `{H_const:.2f}` km/s/Mpc; χ² = {chi2_null:.2f}", f"- Linear-in-rank fit: slope = `{slope:+.4f} ± {s_slope:.4f}` km/s/Mpc/rank, intercept `{intercept:.2f}`", f"- χ²(linear) = {chi2_lin:.2f}; Δχ² = {delta_chi2:.2f}", f"- Slope σ from zero: `{slope_sigma:.2f}σ`", f"- Fit endpoints: H@rank0 = {H_at_rank0:.2f}, H@rank8.5 = {H_at_rank8p5:.2f}, ratio = {R_fit:.4f}", f"- Jackknife slope range: [{slopes.min():.4f}, {slopes.max():.4f}], std {slopes.std(ddof=1):.4f}", "", "**Post-hoc rank caveat:**", out["seventeen_method_test"]["post_hoc_rank_caveat"], "", "## BAO coupling caveat", "", out["BAO_coupling_caveat"], "", "## Falsifier signature", "", f"A measured ratio outside `1.0830 ± 0.005` at >5σ falsifies DCT's " "conformal-frame H₀ mechanism as applied to this two-method ratio. " "Today's two-method ratio matches DCT at the 0.02σ level. This ratio " "test is logically separate from photon null BAO χ(z); see BAO caveat " "above.", ] (HERE / "result.md").write_text("\n".join(md) + "\n") # figure fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(8, 6.5)) ax1.errorbar([row[3] for row in COMPILATION], H, yerr=sH, fmt="o", color="C0", ecolor="gray", capsize=3, alpha=0.85) ax1.axhline(H_const, color="C2", ls="--", label=f"const-H fit = {H_const:.2f}") xx = np.linspace(0, 9, 100) ax1.plot(xx, intercept + slope*xx, color="C3", label=f"linear fit slope={slope:+.3f}") ax1.axhline(H_PLANCK, color="C0", ls=":", alpha=0.5) ax1.axhline(H_SHOES, color="C3", ls=":", alpha=0.5) ax1.set_xlabel("Information-density rank (post-hoc)") ax1.set_ylabel("H₀ (km/s/Mpc)") ax1.set_title("17-method H₀ vs information-density rank (caveat: post-hoc)") ax1.legend(fontsize=8) ax2.axvline(RATIO_DCT, color="C0", lw=2, label=f"DCT 1/√P₀ = {RATIO_DCT:.4f}") ax2.axvline(R_obs, color="C3", lw=2, label=f"H_SH0ES/H_Planck = {R_obs:.4f}") ax2.axvspan(R_obs - sR, R_obs + sR, alpha=0.2, color="C3", label="±1σ") ax2.set_xlim(1.05, 1.13) ax2.set_yticks([]) ax2.set_title("Two-method ratio test") ax2.set_xlabel("H_phys / H_E") ax2.legend() fig.tight_layout() fig.savefig(HERE / "figure.png", dpi=140) plt.close(fig) print(f"Two-method: ratio_obs={R_obs:.4f}+/-{sR:.4f}, DCT={RATIO_DCT:.4f}, " f"sigma_dist={sigma_dist:+.2f} -> {cls}") print(f"17-method: constant chi2={chi2_null:.2f}, linear chi2={chi2_lin:.2f}, " f"slope={slope:+.4f}+/-{s_slope:.4f} -> {slope_sigma:.2f}σ") if __name__ == "__main__": main()