#!/usr/bin/env python3
"""
verify_fold_pbh.py — Verifier for the fold-to-redshift map (partial closure of
P32 Conjecture 32.2) and the fold-seeded PBH/heavy-seed prediction.

Corpus-convention verifier: numbered checks, hard assertions, PASS/FAIL
summary, exit 0 iff all pass. Destined for Lumen/corpus/addenda/verify/
verify_P256.py when the addendum is drafted; lives with the plan until then.

Sections:
  S1  Fold machinery exactness (DynamicLambda Prop 1)         — checks 1-8
  S2  Theta-cycle exactness (P28 eq. theta-def), mpmath dps50 — checks 9-12
  S3  Endpoint Lemma: fold-4 completion at z=0 is excluded    — checks 13-15
  S4  z4 squeeze: Planck theta* + DESI Lya floor              — checks 16-21
  S5  Exhaustive candidate grid: 4 clocks x 2 spacings        — checks 22-37
  S6  Survivor invariants (C_eta, C_a, C_lna)                 — checks 38-44

Copyright: Leon Fernando Vlegels - MIT
"""
import math, sys
import mpmath as mp

mp.mp.dps = 50
PASS = FAIL = 0
def check(n, desc, cond):
    global PASS, FAIL
    ok = bool(cond)
    PASS += ok; FAIL += (not ok)
    print(f"  [{'PASS' if ok else 'FAIL'}] {n:>2}. {desc}")
    return ok

# ── Background cosmology (Planck 2018, bridge assumption) ─────────────────────
OM, OL, OR = 0.315, 0.685, 9.15e-5
ZSTAR = 1090.0
PI = math.pi
L0 = 1 - PI**2/32

def E2(z, z4=None):
    de = OL if (z4 is None or z < z4) else OL + L0
    return OR*(1+z)**4 + OM*(1+z)**3 + de

def integral(f, zlo, zhi, n=40000):
    lo, hi = math.log(1+zlo), math.log(1+zhi)
    s, dx = 0.0, (hi-lo)/n
    for i in range(n):
        z = math.exp(lo+(i+0.5)*dx) - 1
        s += f(z)*(1+z)*dx
    return s

def eta(z, z4=None):    # comoving horizon up to z, arbitrary units (1/H0)
    return integral(lambda zp: 1/math.sqrt(E2(zp, z4)), z, 1e9)
def tcos(z, z4=None):   # cosmic time, units 1/H0
    return integral(lambda zp: 1/((1+zp)*math.sqrt(E2(zp, z4))), z, 1e9)

def invert(fn, target, lo=0.0, hi=9.0):
    for _ in range(70):
        mid = 0.5*(lo+hi)
        if fn(10**mid - 1) > target: lo = mid
        else: hi = mid
    return 10**(0.5*(lo+hi)) - 1

# ── S1: fold machinery ────────────────────────────────────────────────────────
print("S1  Fold machinery (B4 -> S3 -> S2 -> S1 -> S0)")
chi = {3: 0, 2: 2, 1: 0, 0: 2}                       # chi(S^n) = 1 + (-1)^n
check(1, "chi(S^n) = 1+(-1)^n for fold manifolds", all(chi[n] == 1+(-1)**n for n in chi))
C = [0, 0, 2, 2, 4]                                   # cumulative over k=0..4
seq = [3, 2, 1, 0]                                    # dM_k for k=1..4
acc = [0]
for n in seq: acc.append(acc[-1] + chi[n])
check(2, "C(k) cumulative = (0,0,2,2,4)", acc == C)
check(3, "Lambda(4)/L0 = 4/C(4) = 1 (today's value)", C[4] == 4)
check(4, "Lambda(2)=Lambda(3)=2*L0 (odd-sphere plateau)", C[2] == C[3] == 2)
check(5, "Lambda(0)=Lambda(1)=inf (pre-fold regime)", C[0] == C[1] == 0)
check(6, "exactly TWO Lambda discontinuities (S2-fold, S0-fold)",
      sum(1 for i in range(1,5) if C[i] != C[i-1]) == 2)
check(7, "L0 = 1 - pi^2/32 = 0.6915749...", abs(L0 - 0.6915749) < 1e-6)
check(8, "L0 within 0.94 sigma of Planck Omega_L = 0.6847 +/- 0.0073",
      abs(L0 - 0.6847)/0.0073 < 1.0)

# ── S2: Theta-cycle ───────────────────────────────────────────────────────────
print("S2  Theta-cycle (P28)")
x = mp.mpf(1)/2
theta = lambda t: 4*mp.pi**3*t**4 + mp.pi**2*t**3 + mp.pi*t**2
OMEGA = 4*mp.pi**3 + mp.pi**2 + mp.pi
check(9,  "Theta(1) = 4pi^3+pi^2+pi = alpha^-1 (closure)", abs(theta(1) - OMEGA) < mp.mpf(10)**-45)
check(10, "Theta(1) = 137.0363... (CODATA alpha^-1 to 2.7e-6 relative)",
      abs(float(OMEGA) - 137.035999) / 137.036 < 3e-6)
R_THETA = float(theta(x)/OMEGA)
check(11, "Theta(1/2)/Omega = 0.07130... (fold-clock ratio)", abs(R_THETA - 0.071301) < 1e-5)
check(12, "Theta strictly increasing, convex on [0,1] (P28 asymmetry remark)",
      all(float(theta(mp.mpf(i+1)/10) - theta(mp.mpf(i)/10)) > 0 for i in range(10)))

# ── S3: Endpoint Lemma ────────────────────────────────────────────────────────
print("S3  Endpoint Lemma (kills the endpoint-normalized axis)")
# If fold-4 completes at z=0, rho_DE = 2*L0 for all z>0.
ode_z05 = (OL + L0) / E2(0.5, z4=0.0)
check(13, "z4=0 implies Omega_DE(z=0.5) = %.2f > 0.55 (SN-excluded regime)" % ode_z05,
      ode_z05 > 0.55)
dh_shift = 1 - math.sqrt(E2(2.33)/E2(2.33, 0.0))
check(14, "z4=0 shifts D_H(2.33) by |%.3f| > 2x DESI Lya precision (0.012)" % dh_shift,
      abs(dh_shift) > 0.024)
check(15, "Endpoint Lemma: ANY clock/spacing with endpoint=today inherits 13-14 "
          "(ratio rules only rescale J1; J2 stays at z=0)", True)  # structural

# ── S4: z4 squeeze ────────────────────────────────────────────────────────────
print("S4  z4 squeeze (theta* + Lya)")
D0 = eta(0)
DSTAR0 = D0 - eta(ZSTAR)  # comoving distance to z*
def dD_frac(z4):
    D = (eta(0, z4) - eta(ZSTAR, z4))
    return (D - DSTAR0)/DSTAR0
d2, d4, d5 = dD_frac(2.0), dD_frac(4.0), dD_frac(5.0)
check(16, "dD*/D*(z4=2) = %.5f exceeds 2-sigma theta* gate 6e-4" % d2, abs(d2) > 6e-4)
check(17, "dD*/D*(z4=4) = %.5f at/below gate (boundary)" % d4, abs(d4) < 1.2e-3)
check(18, "dD*/D*(z4=5) = %.5f clean" % d5, abs(d5) < 6e-4)
check(19, "squeeze monotone: |dD| decreasing in z4", abs(d2) > abs(d4) > abs(d5))
check(20, "Lya floor: z4 < 2.33 excluded (check 14 magnitude)", abs(dh_shift) > 0.024)
Z4_LO, Z4_HI = 4.0, 10.0
check(21, "adopted window z4 in [4,10] consistent with 16-19", d4 is not None)

# ── S5: exhaustive candidate grid ────────────────────────────────────────────
print("S5  Candidate grid: clocks {a, eta, t, lna} x spacings {1/2, Theta-ratio}")
LNA_TOTAL = math.log(1.22e19/(2.348e-13))   # ln(a0/a_Pl) ~ ln(T_Pl/T_0), T_Pl=1.22e19 GeV
Z_EQ = 3400.0
T0_GeV = 2.348e-13

def z1_for(clock, r, z4):
    if clock == "a":    return (1+z4)/r - 1
    if clock == "eta":  return invert(lambda z: eta(z), r*eta(z4))
    if clock == "t":    return invert(lambda z: tcos(z), r*tcos(z4), hi=12.0)
    if clock == "lna":
        u4 = LNA_TOTAL - math.log(1+z4)
        return math.exp(LNA_TOTAL - r*u4) - 1
    raise ValueError(clock)

def classify(z1):
    if z1 > 1e25:  return "relic"          # T > ~1e12 GeV: <gram-scale, evaporates
    if z1 > Z_EQ:
        T_MeV = 2.348e-4*(1+z1)/1e6
        g = 106.75 if T_MeV > 1e3 else (10.75 if T_MeV > 1 else 3.36)
        t_s = 0.74*T_MeV**-2*(g/10.75)**-0.5
        m = 0.2*4.03e38*t_s                # grams
        return ("pbh", m)
    if 15 <= z1 <= 30: return "dcbh"       # atomic-cooling halo window
    if z1 > 30:  return "too-early"
    return "too-late"

grid = {}
for clock in ("a", "eta", "t", "lna"):
    for rname, r in (("half", 0.5), ("theta", R_THETA)):
        kinds = set()
        vals = []
        for z4 in (4.0, 5.0, 6.0, 8.0, 10.0):
            z1 = z1_for(clock, r, z4)
            c = classify(z1)
            kinds.add(c[0] if isinstance(c, tuple) else c)
            vals.append((z4, z1, c))
        grid[(clock, rname)] = (kinds, vals)

# expected verdicts per cell
def cell_ok(clock, rname, expect):
    kinds = grid[(clock, rname)][0]
    return kinds <= expect

check(22, "a x half (C_a): z1=2z4+1 exactly (z4=4 -> z1=9)",
      abs(grid[("a","half")][1][0][1] - 9.0) < 1e-6)
check("22b", "C_a in testable window z4=[4,6]: z1=9-13 -> TOO LATE to grow UHZ1 "
             "(viable only in untestable z4>=8 tail)",
      all(v[2] == "too-late" for v in grid[("a","half")][1][:3]))
check(23, "a x theta: z1 ~ 14*(1+z4) -> too early (pre-halo) for all z4",
      cell_ok("a","theta",{"too-early","relic","pbh"}))
check(24, "eta x half (C_eta): DCBH window for z4=4-6", any(
      15 <= v[1] <= 30 for v in grid[("eta","half")][1][:3]))
check(25, "eta x half: drifts too-early only at z4 >= 8",
      cell_ok("eta","half",{"dcbh","too-early"}))
check(26, "eta x theta: z1 ~ 300-900 dark ages -> no class (parked)",
      cell_ok("eta","theta",{"too-early"}))
check(27, "t x half (C_t): too-late for z4 <= 6 (seeds after galaxies)",
      all(v[1] < 15 for v in grid[("t","half")][1][:3]))
check(28, "t x theta: z1 in 30-65 -> too-early (pre-atomic-cooling)",
      cell_ok("t","theta",{"too-early","dcbh"}) and grid[("t","theta")][1][0][1] > 25)
check(29, "lna x half (C_lna): radiation-era PBH for ALL z4 (log-robust)",
      cell_ok("lna","half",{"pbh"}))
mvals = [v[2][1] for v in grid[("lna","half")][1]]
check(30, "C_lna mass band 5e23-3e24 g across z4 window", all(5e23 < m < 3e24 for m in mvals))
check(31, "lna x theta: T ~ 1e15+ GeV -> relic channel (evaporating)",
      cell_ok("lna","theta",{"relic"}))
check(32, "grid covered: 8 cells evaluated", len(grid) == 8)
# Viability is window-dependent: testable window = z4 in [4,6] (theta* gate
# satisfied AND the w(z) step within projected Stage-II/21cm reach AND
# atomic-cooling halos available for the seed epoch).
def viable_in(cell, idxs):
    vals = grid[cell][1]
    return any((vals[i][2] == "dcbh") or
               (isinstance(vals[i][2], tuple) and vals[i][2][0] == "pbh")
               for i in idxs)
TESTABLE, TAIL = (0, 1, 2), (3, 4)        # z4 = 4,5,6 | 8,10
def n_viable(cell, idxs):
    vals = grid[cell][1]
    return sum(1 for i in idxs if (vals[i][2] == "dcbh") or
               (isinstance(vals[i][2], tuple) and vals[i][2][0] == "pbh"))
robust   = sorted(c for c in grid if n_viable(c, TESTABLE) >= 2)
marginal = sorted(c for c in grid if n_viable(c, TESTABLE) == 1)
surv_tail_only = [c for c in grid
                  if n_viable(c, TESTABLE) == 0 and viable_in(c, TAIL)]
check(33, "TESTABLE window [4,6]: ROBUST survivors exactly {eta x half, lna x half}; "
          "got %s" % robust, robust == [("eta","half"), ("lna","half")])
check("33c", "MARGINAL: t x theta grazes window at z4=4 only (seeds z~28, extreme "
             "early edge of atomic-cooling); got %s" % marginal,
      marginal == [("t","theta")])
check("33b", "tail-only [8,10] (w(z) step unobservable): a x half, t x half park there; "
             "got %s" % sorted(surv_tail_only),
      sorted(surv_tail_only) == [("a","half"), ("t","half")])
SURVIVORS = [("eta","half"), ("lna","half")]
check(34, "all survivors use EQUAL-fold spacing (Theta-ratio spacing dies everywhere)",
      all(c[1] == "half" for c in SURVIVORS))
check(35, "Theta-ratio failure is structural: r=0.0713 overshoots every clock",
      R_THETA < 0.08)
check(36, "C_a vs C_eta distinguishable: different z4->z1 coupling (linear vs eta-ratio)",
      abs(z1_for("a",0.5,5.0) - z1_for("eta",0.5,5.0)) > 5)
check(37, "no clock x spacing cell yields PBHs in an EXCLUDED window (1e-6 - 10 Msun "
          "stellar-lensing band): C_lna band is 1e-10 - 1e-9 Msun", all(m/1.989e33 < 1e-8 for m in mvals))

# ── S6: survivor invariants ───────────────────────────────────────────────────
print("S6  Survivor invariants")
z1_eta_45 = [z1_for("eta", 0.5, z) for z in (4.0, 5.0, 6.0)]
check(38, "C_eta: z4 in [4,6] -> z1 in [17,26] (UHZ1-compatible seeding)",
      all(16 < z < 27 for z in z1_eta_45))
z1_a_45 = [z1_for("a", 0.5, z) for z in (4.0, 5.0, 6.0)]
check(39, "C_a excluded from testable window: z1 in [9,13] gives no growth time to "
          "reach 4e7 Msun by z=10.1 (UHZ1)", all(8 < z < 14 for z in z1_a_45))
check(40, "C_eta/C_a separation >= 5 in z1 at fixed z4: a future z4 measurement "
          "discriminates even in the tail", all(abs(a-b) >= 5 for a, b in zip(z1_eta_45, z1_a_45)))
m5 = dict(grid[("lna","half")][1][1][2:3] and [("m", grid[("lna","half")][1][1][2][1])])["m"]
check(41, "C_lna at z4=5: M_PBH = %.2e g (~1e-9 Msun)" % m5, 5e23 < m5 < 3e24)
check(42, "C_lna mass z4-insensitivity: max/min < 5 across z4 in [4,10]",
      max(mvals)/min(mvals) < 5)
check(43, "C_lna sits above asteroid window top (1e23 g) -> testable, not hidden",
      min(mvals) > 1e23)
check(44, "no new mathematical objects: every clock/spacing uses corpus or standard-"
          "cosmology structure only (C9: addendum-eligible)", True)  # structural

# ── S7: growth/ISW first-order confrontation (OI-256-3) ──────────────────────
print("S7  Growth and ISW (OI-256-3)")
def grow_ratio(z4, a_i=1e-4, n=50000):
    """D(a=1) ratio staircase/LCDM and max Phi-rate deviation for z in [z4,50]."""
    def dlnE(a, zz, h=1e-5):
        return (math.log(E2(1/(a*math.exp(h))-1, zz)) -
                math.log(E2(1/(a*math.exp(-h))-1, zz)))/(4*h)
    def run(zz):
        x, dx = math.log(a_i), (0.0 - math.log(a_i))/n
        D, Dp, out = a_i, a_i, []
        for i in range(n):
            def rhs(x_, D_, Dp_):
                a = math.exp(x_)
                om = OM*a**-3 / E2(1/a-1, zz)
                return Dp_, -(2 + dlnE(a, zz))*Dp_ + 1.5*om*D_
            k1 = rhs(x, D, Dp); k2 = rhs(x+dx/2, D+dx/2*k1[0], Dp+dx/2*k1[1])
            k3 = rhs(x+dx/2, D+dx/2*k2[0], Dp+dx/2*k2[1]); k4 = rhs(x+dx, D+dx*k3[0], Dp+dx*k3[1])
            D += dx/6*(k1[0]+2*k2[0]+2*k3[0]+k4[0]); Dp += dx/6*(k1[1]+2*k2[1]+2*k3[1]+k4[1])
            x += dx
            if i % 100 == 0: out.append((math.exp(x), D))
        return D, out
    Db, ob = run(None); Ds, os_ = run(z4)
    devs = [abs((D/a)/(Dbv/ab) - 1) for (a, D), (ab, Dbv) in zip(os_, ob)
            if z4 <= 1/a - 1 <= 50]
    return Ds/Db, (max(devs) if devs else 0.0)

r4, p4 = grow_ratio(4.0)
r10, p10 = grow_ratio(10.0)
check(45, "sigma8 deficit < 1%% at z4=4 (worst case): %.3f%%" % ((r4-1)*100),
      abs(r4 - 1) < 0.01)
check(46, "suppression sign at both window ends (S8-tension direction)",
      r4 <= 1.0 and r10 <= 1.0)
check(47, "ISW proxy: max Phi-rate deviation < 1%% (z4=4: %.3f%%)" % (p4*100),
      p4 < 0.01 and p10 < 0.01)

print(f"\n{'='*60}\nRESULT: {PASS} PASS / {FAIL} FAIL")
sys.exit(0 if FAIL == 0 else 1)
