#!/usr/bin/env python3
"""
Dimensional Coherence Theory (DCT) - Cosmological Tensions Analysis
by Nolan G. Parrott

Analyzes DCT predictions against three major cosmological tensions:
1. Hubble Tension (H0 early vs late)
2. DESI BAO (dark energy equation of state)
3. S8 Tension (structure growth suppression)
"""

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec
from scipy.optimize import minimize
from scipy.interpolate import interp1d

# ============================================================
# Published Measurements
# ============================================================

# H0 measurements: (name, value, sigma, type)
H0_data = {
    'Planck':   (67.4, 0.5,  'early'),
    'SH0ES':    (73.04, 1.04, 'late'),
    'TRGB':     (69.8, 1.7,  'late'),
    'TDCOSMO':  (74.2, 1.6,  'late'),
    'H0LiCOW':  (73.3, 1.8,  'late'),
    'MCP':      (73.9, 3.0,  'late'),
}

# S8 measurements: (name, value, sigma, type)
S8_data = {
    'Planck':  (0.832, 0.013, 'early'),
    'DES_Y3':  (0.776, 0.017, 'late'),
    'KiDS':    (0.766, 0.017, 'late'),
    'HSC':     (0.769, 0.033, 'late'),
}

# DESI BAO w0-wa
DESI_w0 = -0.727
DESI_w0_err = 0.067
DESI_wa = -1.05
DESI_wa_err = 0.31

# BBN constraint
BBN_DG_limit = 0.05  # |DeltaG/G| < 0.05

# Cosmological constants
H0_planck = 67.4   # km/s/Mpc
H0_late_avg = 73.0  # approximate late-universe average
t_universe = 13.8   # Gyr
H0_to_invGyr = 1.0 / 978.0  # km/s/Mpc -> 1/Gyr

# ============================================================
# DCT Model Functions
# ============================================================

def P_func(t, alpha, t0):
    """DCT coherence parameter P(t) = 1 - alpha * exp(-t/t0)"""
    return 1.0 - alpha * np.exp(-t / t0)

def Pdot_func(t, alpha, t0):
    """Time derivative of P(t)"""
    return (alpha / t0) * np.exp(-t / t0)

def H_tilde(H_bare, t, alpha, t0):
    """DCT physical Hubble parameter: H_tilde = (1/sqrt(P)) * (H + Pdot/(2P))"""
    P = P_func(t, alpha, t0)
    Pdot = Pdot_func(t, alpha, t0)
    if P <= 0:
        return np.inf
    return (1.0 / np.sqrt(P)) * (H_bare + Pdot / (2.0 * P))

def H_tilde_array(H_bare, t_arr, alpha, t0):
    """Vectorized H_tilde"""
    P = P_func(t_arr, alpha, t0)
    Pdot = Pdot_func(t_arr, alpha, t0)
    mask = P > 0
    result = np.full_like(t_arr, np.inf)
    result[mask] = (1.0 / np.sqrt(P[mask])) * (H_bare + Pdot[mask] / (2.0 * P[mask]))
    return result

def z_to_t(z):
    """Approximate lookback time mapping: t(z) ~ t_universe / (1+z)^(3/2) (matter-dominated approx)"""
    return t_universe / (1.0 + z)**1.5

def dct_weff(z, alpha, t0):
    """DCT effective equation of state: w_eff(z) = -2/3 + 1/(3*H*t0)
    where H is the Hubble rate at redshift z."""
    t = z_to_t(z)
    # H(z) in 1/Gyr units, approximate as H0*(1+z)^(3/2) for matter era
    H_z = H0_planck * H0_to_invGyr * (1.0 + z)**1.5
    correction = 1.0 / (3.0 * H_z * t0)
    return -2.0/3.0 + correction

def desi_w(z):
    """DESI CPL parameterization: w(z) = w0 + wa * z/(1+z)"""
    return DESI_w0 + DESI_wa * z / (1.0 + z)

def dct_s8_suppression(alpha, t0):
    """DCT structure growth suppression factor.

    The coherence drift term Pdot/(2P) acts as an effective friction
    on structure growth, suppressing sigma8 relative to Planck prediction.

    Suppression ~ sqrt(P_late / P_early) * exp(-integral of Pdot/(2P*H) dt)
    For small alpha: suppression ~ 1 - alpha/2 * (1 - exp(-t_uni/t0))
    """
    P_late = P_func(t_universe, alpha, t0)
    P_early = P_func(0.38, alpha, t0)  # t ~ 0.38 Gyr at z~1100 (recombination)

    # Numerical integration of the suppression integral
    t_arr = np.linspace(0.38, t_universe, 1000)
    dt = t_arr[1] - t_arr[0]
    P_arr = P_func(t_arr, alpha, t0)
    Pdot_arr = Pdot_func(t_arr, alpha, t0)
    # H(t) ~ H0 * (t_universe/t)^(2/3) rough approximation
    H_arr = H0_planck * H0_to_invGyr * (t_universe / t_arr)**(2.0/3.0)

    integrand = Pdot_arr / (2.0 * P_arr * H_arr)
    integral = np.trapz(integrand, t_arr)

    suppression = np.sqrt(P_late) * np.exp(-integral)
    return suppression

def bbn_constraint(alpha, t0):
    """BBN constraint: |DeltaG/G| = |1/P - 1| at BBN epoch (t ~ 0.01 Gyr, T ~ 1 MeV)"""
    t_bbn = 0.01  # Gyr (approximate)
    P_bbn = P_func(t_bbn, alpha, t0)
    if P_bbn <= 0:
        return np.inf
    return abs(1.0 / P_bbn - 1.0)

# ============================================================
# Chi-squared Functions
# ============================================================

def chi2_hubble(alpha, t0):
    """Chi-squared for Hubble tension data under DCT."""
    chi2 = 0.0

    # At early universe (CMB, t ~ 0.38 Gyr), H_tilde should match Planck
    t_early = 0.38
    H_tilde_early = H_tilde(H0_planck, t_early, alpha, t0)
    # Planck measurement
    chi2 += ((H_tilde_early - H0_data['Planck'][0]) / H0_data['Planck'][1])**2

    # At late universe (t ~ t_universe), H_tilde should match late measurements
    t_late = t_universe
    H_tilde_late = H_tilde(H0_planck, t_late, alpha, t0)

    for name, (val, sig, typ) in H0_data.items():
        if typ == 'late':
            chi2 += ((H_tilde_late - val) / sig)**2

    return chi2

def chi2_hubble_lcdm():
    """Chi-squared for Hubble data under LCDM (H0 = Planck value everywhere)."""
    chi2 = 0.0
    for name, (val, sig, typ) in H0_data.items():
        chi2 += ((H0_planck - val) / sig)**2
    return chi2

def chi2_desi(alpha, t0):
    """Chi-squared for DESI BAO data under DCT."""
    # Compare DCT w_eff at z=0 and its derivative to DESI w0, wa
    w_dct_0 = dct_weff(0.001, alpha, t0)  # near z=0

    # Compute effective wa from DCT
    z_arr = np.array([0.001, 0.3, 0.5, 0.7, 1.0])
    w_dct_arr = np.array([dct_weff(z, alpha, t0) for z in z_arr])
    w_desi_arr = np.array([desi_w(z) for z in z_arr])

    # Fit DCT w(z) to CPL form to extract effective w0, wa
    # w_CPL(z) = w0 + wa * z/(1+z)
    # Use least squares
    A = np.column_stack([np.ones(len(z_arr)), z_arr / (1.0 + z_arr)])
    result = np.linalg.lstsq(A, w_dct_arr, rcond=None)
    w0_dct, wa_dct = result[0]

    chi2 = ((w0_dct - DESI_w0) / DESI_w0_err)**2
    chi2 += ((wa_dct - DESI_wa) / DESI_wa_err)**2
    return chi2

def chi2_desi_lcdm():
    """Chi-squared for DESI data under LCDM (w=-1, wa=0)."""
    chi2 = ((-1.0 - DESI_w0) / DESI_w0_err)**2
    chi2 += ((0.0 - DESI_wa) / DESI_wa_err)**2
    return chi2

def chi2_s8(alpha, t0):
    """Chi-squared for S8 data under DCT."""
    suppression = dct_s8_suppression(alpha, t0)
    s8_dct = S8_data['Planck'][0] * suppression  # DCT prediction

    chi2 = 0.0
    for name, (val, sig, typ) in S8_data.items():
        if typ == 'late':
            chi2 += ((s8_dct - val) / sig)**2
    # Also penalize deviation from Planck at early times (should be close)
    chi2 += ((S8_data['Planck'][0] - S8_data['Planck'][0]) / S8_data['Planck'][1])**2  # trivially 0
    return chi2

def chi2_s8_lcdm():
    """Chi-squared for S8 data under LCDM (no suppression)."""
    chi2 = 0.0
    s8_lcdm = S8_data['Planck'][0]
    for name, (val, sig, typ) in S8_data.items():
        chi2 += ((s8_lcdm - val) / sig)**2
    return chi2

def chi2_bbn(alpha, t0):
    """BBN constraint penalty."""
    dg = bbn_constraint(alpha, t0)
    if dg > BBN_DG_limit:
        return ((dg - BBN_DG_limit) / 0.01)**2  # steep penalty
    return 0.0

def chi2_total(params):
    """Total joint chi-squared for DCT."""
    alpha, t0 = params
    if alpha < 0 or alpha > 1 or t0 < 0.1:
        return 1e10
    return chi2_hubble(alpha, t0) + chi2_desi(alpha, t0) + chi2_s8(alpha, t0) + chi2_bbn(alpha, t0)

# ============================================================
# Parameter Scan
# ============================================================

print("=" * 70)
print("DCT Cosmological Tensions Analysis")
print("Dimensional Coherence Theory by Nolan G. Parrott")
print("=" * 70)

# Scan alpha and t0
N_alpha = 200
N_t0 = 200
alpha_arr = np.linspace(0.05, 0.4, N_alpha)
t0_arr = np.linspace(5.0, 200.0, N_t0)

chi2_hubble_grid = np.zeros((N_alpha, N_t0))
chi2_desi_grid = np.zeros((N_alpha, N_t0))
chi2_s8_grid = np.zeros((N_alpha, N_t0))
chi2_bbn_grid = np.zeros((N_alpha, N_t0))
chi2_total_grid = np.zeros((N_alpha, N_t0))

print("\nScanning parameter space: alpha in [0.05, 0.4], t0 in [5, 200] Gyr...")
for i, alpha in enumerate(alpha_arr):
    for j, t0 in enumerate(t0_arr):
        ch = chi2_hubble(alpha, t0)
        cd = chi2_desi(alpha, t0)
        cs = chi2_s8(alpha, t0)
        cb = chi2_bbn(alpha, t0)
        chi2_hubble_grid[i, j] = ch
        chi2_desi_grid[i, j] = cd
        chi2_s8_grid[i, j] = cs
        chi2_bbn_grid[i, j] = cb
        chi2_total_grid[i, j] = ch + cd + cs + cb

# Find best-fit
idx = np.unravel_index(np.argmin(chi2_total_grid), chi2_total_grid.shape)
alpha_best = alpha_arr[idx[0]]
t0_best = t0_arr[idx[1]]
chi2_min = chi2_total_grid[idx]

# Refine with optimizer
from scipy.optimize import minimize
res = minimize(chi2_total, [alpha_best, t0_best], method='Nelder-Mead',
               options={'xatol': 1e-6, 'fatol': 1e-6})
alpha_best, t0_best = res.x
chi2_min = res.fun

print(f"\nBest-fit DCT parameters:")
print(f"  alpha = {alpha_best:.4f}")
print(f"  t0    = {t0_best:.2f} Gyr")
print(f"  Joint chi2 (DCT)  = {chi2_min:.2f}")

# Compute individual contributions at best fit
ch_best = chi2_hubble(alpha_best, t0_best)
cd_best = chi2_desi(alpha_best, t0_best)
cs_best = chi2_s8(alpha_best, t0_best)
cb_best = chi2_bbn(alpha_best, t0_best)

print(f"\n  chi2 breakdown:")
print(f"    Hubble: {ch_best:.2f}")
print(f"    DESI:   {cd_best:.2f}")
print(f"    S8:     {cs_best:.2f}")
print(f"    BBN:    {cb_best:.2f}")

# LCDM comparison
chi2_lcdm_h = chi2_hubble_lcdm()
chi2_lcdm_d = chi2_desi_lcdm()
chi2_lcdm_s = chi2_s8_lcdm()
chi2_lcdm_total = chi2_lcdm_h + chi2_lcdm_d + chi2_lcdm_s

print(f"\nLCDM chi2:")
print(f"    Hubble: {chi2_lcdm_h:.2f}")
print(f"    DESI:   {chi2_lcdm_d:.2f}")
print(f"    S8:     {chi2_lcdm_s:.2f}")
print(f"    Total:  {chi2_lcdm_total:.2f}")

delta_chi2 = chi2_lcdm_total - chi2_min
print(f"\nDelta chi2 (LCDM - DCT) = {delta_chi2:.2f}")
print(f"DCT improvement: {delta_chi2:.1f} units of chi2 (2 extra parameters)")

# Derived quantities at best fit
P_now = P_func(t_universe, alpha_best, t0_best)
H_tilde_now = H_tilde(H0_planck, t_universe, alpha_best, t0_best)
ratio = H_tilde_now / H0_planck
supp = dct_s8_suppression(alpha_best, t0_best)
s8_predicted = 0.832 * supp
w_eff_0 = dct_weff(0.01, alpha_best, t0_best)
dg_bbn = bbn_constraint(alpha_best, t0_best)

print(f"\nDerived quantities at best fit:")
print(f"  P(t_now)        = {P_now:.4f}")
print(f"  H_tilde(now)    = {H_tilde_now:.2f} km/s/Mpc")
print(f"  H_tilde/H_bare  = {ratio:.4f}")
print(f"  S8 suppression  = {supp:.4f}")
print(f"  S8 predicted    = {s8_predicted:.4f}")
print(f"  w_eff(z~0)      = {w_eff_0:.4f}")
print(f"  |DG/G| at BBN   = {dg_bbn:.5f} (limit: {BBN_DG_limit})")

# ============================================================
# Multi-Panel Figure
# ============================================================

fig = plt.figure(figsize=(18, 14))
fig.suptitle('Dimensional Coherence Theory: Cosmological Tensions Analysis\nNolan G. Parrott',
             fontsize=16, fontweight='bold', y=0.98)

gs = GridSpec(3, 3, figure=fig, hspace=0.35, wspace=0.35)

# --- Panel 1: Chi-squared surface (Hubble) ---
ax1 = fig.add_subplot(gs[0, 0])
AA, TT = np.meshgrid(alpha_arr, t0_arr, indexing='ij')
levels_h = np.linspace(chi2_hubble_grid.min(), chi2_hubble_grid.min() + 30, 15)
cp1 = ax1.contourf(AA, TT, chi2_hubble_grid, levels=levels_h, cmap='viridis')
ax1.plot(alpha_best, t0_best, 'r*', markersize=15, label='Best fit')
ax1.set_xlabel(r'$\alpha$', fontsize=12)
ax1.set_ylabel(r'$t_0$ [Gyr]', fontsize=12)
ax1.set_title(r'$\chi^2$ Surface: Hubble Tension', fontsize=11)
plt.colorbar(cp1, ax=ax1, label=r'$\chi^2$')
ax1.legend(fontsize=9)

# --- Panel 2: Chi-squared surface (Joint) ---
ax2 = fig.add_subplot(gs[0, 1])
levels_j = np.linspace(chi2_total_grid.min(), chi2_total_grid.min() + 50, 20)
cp2 = ax2.contourf(AA, TT, chi2_total_grid, levels=levels_j, cmap='magma')
ax2.plot(alpha_best, t0_best, 'c*', markersize=15, label='Best fit')
# Draw 1-sigma, 2-sigma contours (Delta chi2 = 2.30, 6.18 for 2 params)
ax2.contour(AA, TT, chi2_total_grid, levels=[chi2_min + 2.30, chi2_min + 6.18],
            colors=['cyan', 'white'], linewidths=[1.5, 1.0], linestyles=['solid', 'dashed'])
ax2.set_xlabel(r'$\alpha$', fontsize=12)
ax2.set_ylabel(r'$t_0$ [Gyr]', fontsize=12)
ax2.set_title(r'Joint $\chi^2$ Surface (all tensions)', fontsize=11)
plt.colorbar(cp2, ax=ax2, label=r'$\chi^2$')
ax2.legend(fontsize=9)

# --- Panel 3: H0 measurements vs DCT prediction ---
ax3 = fig.add_subplot(gs[0, 2])
colors_h0 = {'early': '#2196F3', 'late': '#F44336'}
y_pos = 0
labels = []
for name, (val, sig, typ) in H0_data.items():
    ax3.errorbar(val, y_pos, xerr=sig, fmt='o', color=colors_h0[typ],
                 markersize=8, capsize=4, capthick=2, linewidth=2)
    labels.append(name)
    y_pos += 1

# DCT predictions
ax3.axvline(H0_planck, color='#2196F3', linestyle='--', alpha=0.5, label=r'$\Lambda$CDM ($H_0=67.4$)')
ax3.axvline(H_tilde_now, color='#4CAF50', linestyle='-', linewidth=2.5,
            label=f'DCT $\\tilde{{H}}_0$={H_tilde_now:.1f}')
ax3.axvspan(H_tilde_now - 1.0, H_tilde_now + 1.0, alpha=0.15, color='green')

ax3.set_yticks(range(len(labels)))
ax3.set_yticklabels(labels, fontsize=10)
ax3.set_xlabel(r'$H_0$ [km/s/Mpc]', fontsize=12)
ax3.set_title('Hubble Constant Measurements', fontsize=11)
ax3.legend(fontsize=8, loc='lower right')
ax3.set_xlim(64, 80)

# --- Panel 4: w(z) comparison ---
ax4 = fig.add_subplot(gs[1, 0])
z_plot = np.linspace(0.01, 2.0, 200)
w_dct_plot = np.array([dct_weff(z, alpha_best, t0_best) for z in z_plot])
w_desi_plot = np.array([desi_w(z) for z in z_plot])

ax4.plot(z_plot, w_dct_plot, 'b-', linewidth=2.5, label='DCT $w_{eff}(z)$')
ax4.plot(z_plot, w_desi_plot, 'r--', linewidth=2, label='DESI CPL fit')
ax4.axhline(-1.0, color='gray', linestyle=':', alpha=0.5, label=r'$\Lambda$CDM ($w=-1$)')
ax4.axhline(-2.0/3.0, color='blue', linestyle=':', alpha=0.3, label='DCT asymptote $w=-2/3$')

# DESI error band at z=0
ax4.errorbar(0.0, DESI_w0, yerr=DESI_w0_err, fmt='rs', markersize=8, capsize=4,
             label=f'DESI $w_0$={DESI_w0}', zorder=5)

ax4.set_xlabel('Redshift $z$', fontsize=12)
ax4.set_ylabel('$w(z)$', fontsize=12)
ax4.set_title('Dark Energy Equation of State', fontsize=11)
ax4.legend(fontsize=8, loc='lower left')
ax4.set_ylim(-1.8, -0.3)
ax4.set_xlim(0, 2.0)

# --- Panel 5: S8 measurements ---
ax5 = fig.add_subplot(gs[1, 1])
y_pos = 0
labels_s8 = []
colors_s8 = {'early': '#2196F3', 'late': '#FF9800'}
for name, (val, sig, typ) in S8_data.items():
    ax5.errorbar(val, y_pos, xerr=sig, fmt='s', color=colors_s8[typ],
                 markersize=8, capsize=4, capthick=2, linewidth=2)
    labels_s8.append(name)
    y_pos += 1

ax5.axvline(0.832, color='#2196F3', linestyle='--', alpha=0.5, label=r'Planck $\Lambda$CDM')
ax5.axvline(s8_predicted, color='#4CAF50', linestyle='-', linewidth=2.5,
            label=f'DCT $S_8$={s8_predicted:.3f}')
ax5.axvspan(s8_predicted - 0.015, s8_predicted + 0.015, alpha=0.15, color='green')

ax5.set_yticks(range(len(labels_s8)))
ax5.set_yticklabels(labels_s8, fontsize=10)
ax5.set_xlabel(r'$S_8$', fontsize=12)
ax5.set_title(r'$S_8$ Measurements vs DCT', fontsize=11)
ax5.legend(fontsize=8, loc='lower right')
ax5.set_xlim(0.7, 0.88)

# --- Panel 6: P(t) evolution ---
ax6 = fig.add_subplot(gs[1, 2])
t_plot = np.linspace(0.01, t_universe, 500)
P_plot = P_func(t_plot, alpha_best, t0_best)
Pdot_plot = Pdot_func(t_plot, alpha_best, t0_best)

ax6.plot(t_plot, P_plot, 'b-', linewidth=2.5, label=r'$P(t)$')
ax6.axhline(1.0, color='gray', linestyle=':', alpha=0.5)
ax6.axvline(t_universe, color='red', linestyle='--', alpha=0.3, label='$t_{now}$')

# Mark key epochs
ax6.axvline(0.38, color='orange', linestyle=':', alpha=0.5, label='CMB ($z\\sim1100$)')
ax6.axvline(0.01, color='purple', linestyle=':', alpha=0.5, label='BBN')

ax6_twin = ax6.twinx()
ax6_twin.plot(t_plot, Pdot_plot, 'r--', linewidth=1.5, alpha=0.7, label=r'$\dot{P}(t)$')
ax6_twin.set_ylabel(r'$\dot{P}(t)$ [Gyr$^{-1}$]', fontsize=11, color='red')
ax6_twin.tick_params(axis='y', labelcolor='red')

ax6.set_xlabel('Cosmic time $t$ [Gyr]', fontsize=12)
ax6.set_ylabel(r'$P(t)$', fontsize=12)
ax6.set_title('DCT Coherence Parameter Evolution', fontsize=11)
ax6.legend(fontsize=8, loc='center right')
ax6.set_xlim(0, t_universe + 0.5)

# --- Panel 7: H_tilde evolution with z ---
ax7 = fig.add_subplot(gs[2, 0])
z_h_plot = np.linspace(0.01, 3.0, 300)
t_h_plot = np.array([z_to_t(z) for z in z_h_plot])
H_bare_z = H0_planck * (1.0 + z_h_plot)**1.5  # matter-dominated approx
H_tilde_z = np.array([H_tilde(H0_planck * (1+z)**1.5, z_to_t(z), alpha_best, t0_best)
                       for z in z_h_plot])
ratio_z = H_tilde_z / H_bare_z

ax7.plot(z_h_plot, ratio_z, 'b-', linewidth=2.5, label=r'$\tilde{H}/H_{bare}$')
ax7.axhline(1.0, color='gray', linestyle=':', alpha=0.5)
ax7.axhline(73.0/67.4, color='red', linestyle='--', alpha=0.5,
            label=f'SH0ES/Planck = {73.0/67.4:.3f}')
ax7.set_xlabel('Redshift $z$', fontsize=12)
ax7.set_ylabel(r'$\tilde{H}/H_{bare}$', fontsize=12)
ax7.set_title('DCT Hubble Enhancement vs Redshift', fontsize=11)
ax7.legend(fontsize=9)
ax7.set_xlim(0, 3)

# --- Panel 8: Chi-squared comparison bar chart ---
ax8 = fig.add_subplot(gs[2, 1])
categories = ['Hubble', 'DESI', 'S8', 'TOTAL']
chi2_dct_vals = [ch_best, cd_best, cs_best, chi2_min]
chi2_lcdm_vals = [chi2_lcdm_h, chi2_lcdm_d, chi2_lcdm_s, chi2_lcdm_total]

x = np.arange(len(categories))
width = 0.35
bars1 = ax8.bar(x - width/2, chi2_lcdm_vals, width, label=r'$\Lambda$CDM',
                color='#F44336', alpha=0.8, edgecolor='black')
bars2 = ax8.bar(x + width/2, chi2_dct_vals, width, label='DCT',
                color='#4CAF50', alpha=0.8, edgecolor='black')

# Add value labels
for bar in bars1:
    ax8.text(bar.get_x() + bar.get_width()/2., bar.get_height() + 0.3,
             f'{bar.get_height():.1f}', ha='center', va='bottom', fontsize=8)
for bar in bars2:
    ax8.text(bar.get_x() + bar.get_width()/2., bar.get_height() + 0.3,
             f'{bar.get_height():.1f}', ha='center', va='bottom', fontsize=8)

ax8.set_xticks(x)
ax8.set_xticklabels(categories, fontsize=11)
ax8.set_ylabel(r'$\chi^2$', fontsize=12)
ax8.set_title(r'$\chi^2$ Comparison: $\Lambda$CDM vs DCT', fontsize=11)
ax8.legend(fontsize=10)

# --- Panel 9: Summary text ---
ax9 = fig.add_subplot(gs[2, 2])
ax9.axis('off')
# (summary text is rendered below using plain text to avoid LaTeX issues)
# Use a simpler text approach to avoid LaTeX issues
ax9.text(0.05, 0.95, f"DCT Best-Fit Parameters\n{'='*30}",
         transform=ax9.transAxes, fontsize=12, verticalalignment='top',
         fontfamily='monospace', fontweight='bold')
summary_text = (
    f"\nalpha     = {alpha_best:.4f}\n"
    f"t0        = {t0_best:.1f} Gyr\n"
    f"P(now)    = {P_now:.4f}\n"
    f"H_tilde   = {H_tilde_now:.1f} km/s/Mpc\n"
    f"S8_DCT    = {s8_predicted:.3f}\n"
    f"w_eff(0)  = {w_eff_0:.3f}\n"
    f"|dG/G|    = {dg_bbn:.4f}\n\n"
    f"Chi-squared Comparison\n{'='*30}\n"
    f"LCDM total: {chi2_lcdm_total:.1f}\n"
    f"DCT total:  {chi2_min:.1f}\n"
    f"Delta chi2: {delta_chi2:.1f}\n"
    f"(2 extra params)"
)
ax9.text(0.05, 0.82, summary_text,
         transform=ax9.transAxes, fontsize=10, verticalalignment='top',
         fontfamily='monospace',
         bbox=dict(boxstyle='round', facecolor='lightyellow', alpha=0.8))

import os as _os_outdir
_OUTDIR = _os_outdir.path.join(_os_outdir.path.dirname(_os_outdir.path.abspath(__file__)), 'outputs')
_os_outdir.makedirs(_OUTDIR, exist_ok=True)
plt.savefig(_os_outdir.path.join(_OUTDIR, 'tensions_results.png'),
            dpi=150, bbox_inches='tight', facecolor='white')
print(f"\nFigure saved to tensions_results.png")
print("=" * 70)