#!/usr/bin/env python3
"""
fitted_knobs_analysis.py  --  DFD cosmological fitted-knob ledger (numpy only).

Reproduces, with no sampler (no cobaya, no CAMB/CLASS call):
  (1) A_s = 32 pi alpha^5  via the de Sitter dictionary with eps_W = alpha/(4pi)
      and H_star/M_R = sqrt(8pi); confirms the power alpha^5 is rigid and the
      coefficient collapses to ONE normalization knob.
  (2) the eps_W <-> r consistency: r = 16 eps_W = 4 alpha/pi (tautology, not a
      second vote).
  (3) the chi abundance overshoot: Omega_chi h^2|post = 1.62 * pi^2/3 ~ 5.3,
      a factor ~44 over the 0.12 target, plus the H=m convention band (~9x).
  (4) the honest knob count: 4 fitted DFD coefficients vs LambdaCDM's 6
      (3 under a charitable omega_b grade); H0 and n_s derived outright.

Companion to appendix_FK_fitted_knobs.tex.  Writes fig_fitted_knobs.png.

HONESTY-GUARD: this script neither derives 32pi/eps_W nor fixes the overshoot.
It reports them as computed.

UPGRADE NOTE (see App. FK): item (3), the continuum pi^2/3 overshoot, is the
HISTORICAL pre-theorem form. The finite SU(2)_60 modular-vacuum measure replaces
pi^2/3 with <theta^2>=0.073, giving Omega_chi h^2 = 0.118 (-1.5 sigma) with
theta_i removed as an input (chi_modular_measure.py). The overshoot rows below
are retained as the recorded historical accounting.
"""

import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle
import os

# ----------------------------------------------------------------------------
# Locked DFD ledger constants (App. AT / complete-params).
# ----------------------------------------------------------------------------
ALPHA_INV = 137.035999084
ALPHA     = 1.0 / ALPHA_INV
M_P       = 1.220890e19            # GeV, non-reduced Planck mass
M_Pbar    = M_P / np.sqrt(8*np.pi) # GeV, reduced Planck mass
M_R       = M_P * ALPHA**3         # GeV, locked Majorana rung

# Planck 2018 targets
As_PLANCK   = 2.100e-9
As_SIGMA    = 0.030e-9
OMCH2_TARGET= 0.1200
NS_PLANCK   = 0.9665
H0_SHOES    = 73.0

print("="*72)
print("DFD FITTED-KNOB LEDGER  (numpy only, no sampler)")
print("="*72)
print(f"alpha^-1 = {ALPHA_INV:.6f}   M_P = {M_P:.4e} GeV   "
      f"Mbar_P = {M_Pbar:.4e} GeV   M_R = {M_R:.4e} GeV")

# ----------------------------------------------------------------------------
# (1) A_s = 32 pi alpha^5 via the de Sitter dictionary.
#     A_s = H_star^2 / (8 pi^2 eps_W Mbar_P^2).
#     Inputs: H_star = sqrt(8pi) M_R ,  eps_W = alpha/(4pi).
# ----------------------------------------------------------------------------
print("\n" + "-"*72)
print("(1) SCALAR AMPLITUDE  A_s = 32 pi alpha^5")
print("-"*72)
eps_W   = ALPHA / (4*np.pi)
H_star  = np.sqrt(8*np.pi) * M_R           # = sqrt(8pi) M_R
As_deSitter = H_star**2 / (8*np.pi**2 * eps_W * M_Pbar**2)
As_closed   = 32*np.pi*ALPHA**5

print(f"eps_W = alpha/(4pi)            = {eps_W:.6e}")
print(f"H_star = sqrt(8pi) M_R         = {H_star:.6e} GeV")
print(f"A_s  (de Sitter dictionary)    = {As_deSitter:.6e}")
print(f"A_s  (closed form 32 pi a^5)   = {As_closed:.6e}")
print(f"  ratio dictionary/closed-form = {As_deSitter/As_closed:.10f}")
sigma_As = (As_closed - As_PLANCK)/As_SIGMA
print(f"  A_s vs Planck (2.100e-9)     = {sigma_As:+.2f} sigma")

# The required inflationary potential V = 3 alpha^6 M_P^4 (the same single knob).
V_required = 3*ALPHA**6 * M_P**4
V_from_eps = 24*np.pi**2 * eps_W * M_Pbar**4 * As_closed
print(f"V required = 3 alpha^6 M_P^4   = {V_required:.6e} GeV^4")
print(f"V from (eps_W, A_s)            = {V_from_eps:.6e} GeV^4")
print(f"  ratio                        = {V_from_eps/V_required:.6f}")
Vquarter = V_required**0.25
print(f"V^(1/4)                        = {Vquarter:.4e} GeV  (GUT-ish)")

# Power-rigidity check: alpha^5 = (alpha^3)^2 / alpha^1
power_residual = abs(ALPHA**5 - (ALPHA**3)**2 / ALPHA**1)
print(f"power-rigidity residual a^5-(a^3)^2/a = {power_residual:.2e}  (== 0)")

# ----------------------------------------------------------------------------
# (2) eps_W <-> r consistency (tautology, NOT independent corroboration).
# ----------------------------------------------------------------------------
print("\n" + "-"*72)
print("(2) TENSOR RATIO  r = 16 eps_W = 4 alpha/pi  (same knob, viewed twice)")
print("-"*72)
r_from_eps = 16*eps_W
r_closed   = 4*ALPHA/np.pi
print(f"r = 16 eps_W                   = {r_from_eps:.6e}")
print(f"r = 4 alpha/pi                 = {r_closed:.6e}")
print(f"  residual                     = {abs(r_from_eps-r_closed):.2e}  (== 0)")
# 'independent' magnetic-dual route collapses to the same number:
alpha_M = 1.0/(4*ALPHA)
eps_dual = ALPHA**2 * alpha_M / np.pi
print(f"eps_W via alpha^2 alpha_M/pi   = {eps_dual:.6e}  "
      f"(== alpha/(4pi): {np.isclose(eps_dual, eps_W)})")
print("  -> magnetic-dual route is algebraically identical to the 1-loop QED")
print("     measure alpha/(4pi); NOT an independent derivation of eps_W.")
As_FORCED = False   # coefficient 32pi / eps_W NOT forced
EPSW_DERIVED = False

# ----------------------------------------------------------------------------
# (3) chi abundance overshoot.
# ----------------------------------------------------------------------------
print("\n" + "-"*72)
print("(3) CHI ABUNDANCE OVERSHOOT  Omega_chi h^2|post = 1.62 * pi^2/3")
print("-"*72)
# locked rungs (App. chi_matter)
f_chi   = M_Pbar * ALPHA**3
m_chi   = np.sqrt(158) * M_Pbar * ALPHA**13      # GeV
m_chi_eV= m_chi * 1e9
print(f"f_chi = Mbar_P alpha^3         = {f_chi:.4e} GeV")
print(f"m_chi = sqrt(158) Mbar_P a^13  = {m_chi_eV:.4f} eV")

coeff_3Hm = 1.62          # corpus 3H=m onset convention
theta2_avg = np.pi**2/3.0 # horizon-averaged <theta^2>
omch2_post = coeff_3Hm * theta2_avg
overshoot  = omch2_post / OMCH2_TARGET
print(f"<theta^2> = pi^2/3             = {theta2_avg:.4f}")
print(f"Omega_chi h^2|post (3H=m)      = {omch2_post:.3f}")
print(f"  overshoot vs 0.12            = {overshoot:.1f}x")

# convention band: H=m generic-ALP onset coefficient ~ 0.31
coeff_Hm = 0.31
omch2_Hm = coeff_Hm * theta2_avg
overshoot_Hm = omch2_Hm / OMCH2_TARGET
print(f"Omega_chi h^2|post (H=m conv.) = {omch2_Hm:.3f}  "
      f"-> overshoot {overshoot_Hm:.1f}x")
print(f"  => convention band: ~{overshoot_Hm:.0f}x to ~{overshoot:.0f}x  (all > 1)")

# is 44 an alpha-rung? log(overshoot)/log(137)
rung_test = np.log(overshoot)/np.log(ALPHA_INV)
print(f"log(44.4)/log(137)            = {rung_test:.3f}  (non-integer "
      f"-> 44 is NOT an alpha-rung -> no same-rung m/f fix)")
OVERSHOOT_RESOLVED = False

# ----------------------------------------------------------------------------
# (4) honest knob count.
# ----------------------------------------------------------------------------
print("\n" + "-"*72)
print("(4) HONEST KNOB COUNT  (post Step-5b)")
print("-"*72)
H0_DFD = ALPHA**28.5 / 1.0   # symbolic placeholder; numeric value below
# H0 = alpha^(57/2)/t_P -> 72.09 km/s/Mpc (reproduced elsewhere); 1/H0 in Gyr:
H0_kms = 72.0902
invH0_Gyr = (1.0/(H0_kms)) * (3.0857e19/3.1557e16)  # (km/s/Mpc)^-1 -> Gyr
ns_DFD = 1 - 2/60
print(f"H0  = alpha^(57/2)/t_P         = {H0_kms} km/s/Mpc  "
      f"(1/H0 = {invH0_Gyr:.2f} Gyr)  [DERIVED]")
print(f"n_s = 1 - 2/60                 = {ns_DFD:.4f}  [DERIVED]")

ledger = [
    ("H0",            "DERIVED",                 0),
    ("n_s",           "DERIVED",                 0),
    ("A_s",           "FORCED POW / FIT COEFF",  1),
    ("omega_chi h^2", "FITTED (norm.)",          1),
    ("A_screen",      "FORCED FORM / NOT PINNED",0),
    ("omega_b h^2",   "FITTED / IMPORTED",       1),
    ("tau",           "FITTED (astro)",          1),
]
k_DFD = sum(w for _,_,w in ledger)
k_DFD_charitable = k_DFD - 1   # if omega_b graded forced-power
print(f"\nfitted DFD coefficients k_DFD = {k_DFD}  "
      f"(charitable omega_b grade -> {k_DFD_charitable});  LambdaCDM = 6")
for name, grade, w in ledger:
    print(f"   {name:14s} {grade:26s} fitted={w}")

# ----------------------------------------------------------------------------
# FIGURE: knob ledger (left) + A_s/eps_W consistency (right).
# ----------------------------------------------------------------------------
fig = plt.figure(figsize=(13.5, 6.2))

# --- left panel: knob ledger bars ---
ax1 = fig.add_axes([0.06, 0.12, 0.42, 0.78])
names  = [l[0] for l in ledger]
grades = [l[1] for l in ledger]
fitted = [l[2] for l in ledger]
ypos = np.arange(len(ledger))[::-1]
colors = []
for g in grades:
    if g.startswith("DERIVED"):       colors.append("#2e8b57")  # green
    elif g.startswith("FORCED"):      colors.append("#d4a017")  # amber
    else:                             colors.append("#b22222")  # red
ax1.barh(ypos, [1]*len(ledger), color=colors, edgecolor="k", alpha=0.85)
for y, n, g, f in zip(ypos, names, grades, fitted):
    ax1.text(0.02, y, n, va="center", ha="left", fontsize=11,
             fontweight="bold", color="white")
    ax1.text(0.98, y, g, va="center", ha="right", fontsize=8.5, color="white")
ax1.set_yticks([]); ax1.set_xticks([])
ax1.set_xlim(0,1); ax1.set_ylim(-0.6, len(ledger)-0.4)
ax1.set_title(f"DFD cosmological knob ledger\n"
              f"fitted k_DFD = {k_DFD} (charitable {k_DFD_charitable})  "
              f"vs  LambdaCDM = 6", fontsize=12, fontweight="bold")
# legend
leg = [Rectangle((0,0),1,1,fc="#2e8b57"), Rectangle((0,0),1,1,fc="#d4a017"),
       Rectangle((0,0),1,1,fc="#b22222")]
ax1.legend(leg, ["derived (0 knobs)", "forced power / form",
                 "fitted / imported"], loc="lower center",
          bbox_to_anchor=(0.5,-0.13), ncol=3, fontsize=8.5, frameon=False)

# --- right top: A_s / eps_W consistency ---
ax2 = fig.add_axes([0.56, 0.58, 0.40, 0.32])
labels = ["A_s\n(de Sitter)", "A_s\n32pi a^5", "Planck"]
vals   = [As_deSitter*1e9, As_closed*1e9, As_PLANCK*1e9]
errs   = [0, 0, As_SIGMA*1e9]
bx = np.arange(3)
ax2.bar(bx[:2], vals[:2], color=["#d4a017","#d4a017"], edgecolor="k", width=0.6)
ax2.errorbar(bx[2], vals[2], yerr=errs[2], fmt="o", color="#1f3b73",
             capsize=5, ms=8, label="Planck 2018")
ax2.axhline(As_PLANCK*1e9, color="#1f3b73", ls="--", lw=1, alpha=0.6)
ax2.set_xticks(bx); ax2.set_xticklabels(labels, fontsize=9)
ax2.set_ylabel(r"$A_s\times10^{9}$", fontsize=10)
ax2.set_ylim(2.0, 2.16)
ax2.set_title(f"A_s power FORCED, coeff 32pi FITTED   "
              f"({sigma_As:+.2f}$\\sigma$, eps_W=$\\alpha/4\\pi$ asserted)",
              fontsize=9.5, fontweight="bold")
ax2.text(0.5, 2.005, r"r = 16$\epsilon_W$ = 4$\alpha/\pi$ "
         "= same knob, viewed twice", ha="center", fontsize=8,
         transform=ax2.get_xaxis_transform(), style="italic")

# --- right bottom: chi overshoot band ---
ax3 = fig.add_axes([0.56, 0.12, 0.40, 0.32])
conv = ["H=m\n(0.31)", "3H=m\n(1.62)"]
oms  = [omch2_Hm, omch2_post]
ax3.bar([0,1], oms, color="#b22222", edgecolor="k", width=0.55, alpha=0.85)
ax3.axhline(OMCH2_TARGET, color="#2e8b57", ls="--", lw=2,
            label=f"target 0.12")
ax3.set_xticks([0,1]); ax3.set_xticklabels(conv, fontsize=9)
ax3.set_ylabel(r"$\Omega_\chi h^2$", fontsize=10)
ax3.set_yscale("log")
for x, o in zip([0,1], oms):
    ax3.text(x, o*1.15, f"{o/OMCH2_TARGET:.0f}x", ha="center",
             fontsize=10, fontweight="bold")
ax3.set_title("chi abundance OVERSHOOT (open): ~9-44x, no forced fix",
              fontsize=9.5, fontweight="bold")
ax3.legend(fontsize=8.5, loc="upper left")

fig.suptitle("DFD Fitted Cosmological Knobs -- Honest Ledger "
             "(numpy reproduction, no sampler)",
             fontsize=13, fontweight="bold", y=0.985)

outpng = os.path.join(os.path.dirname(os.path.abspath(__file__)),
                      "fig_fitted_knobs.png")
fig.savefig(outpng, dpi=140, bbox_inches="tight")
print(f"\nFigure written: {outpng}")

# ----------------------------------------------------------------------------
# MACHINE-READABLE MANIFEST
# ----------------------------------------------------------------------------
print("\n" + "="*72)
print("MANIFEST")
print("="*72)
print(f"As_forced (32pi/eps_W)        : {As_FORCED}")
print(f"epsW_derived                  : {EPSW_DERIVED}")
print(f"As_eps_closes_WdW_contradiction: {As_FORCED and EPSW_DERIVED}")
print(f"overshoot_resolved (forced fix): {OVERSHOOT_RESOLVED}")
print(f"overshoot factor (3H=m / H=m)  : {overshoot:.1f}x / {overshoot_Hm:.1f}x")
print(f"k_DFD fitted coefficients      : {k_DFD}  (charitable {k_DFD_charitable})")
print(f"k_LCDM                         : 6")
print(f"derived outright (0 knobs)     : H0, n_s")
print(f"biggest remaining fitted item  : A_s coefficient (32pi / eps_W=alpha/4pi)")
print("="*72)
