#!/usr/bin/env python3
"""
DFD-CLASS -- chi relic abundance from the closed-S^3 modular-vacuum measure.

This is the REPRODUCIBLE backing for the A+ derivation in appendix_chi_matter.tex
(steps 1-2 after thm:chi-stateprep). It supersedes the continuum pi^2/3
"44x overshoot" calculation of chi_postinflation_check.py, which used the WRONG
(continuum classical-angle) measure that the surgery/no-boundary argument excludes.

Physics:
  The chi state space is NOT S^1 with a uniform classical-angle measure; it is the
  finite SU(2)_k Chern-Simons/WZW Hilbert space of (k+1) integrable sectors
  j = 0, 1/2, ..., k/2 at level k=60. The no-boundary closed-S^3 vacuum, prepared
  by Heegaard surgery (two solid tori glued by the modular S generator), has Born
  weights equal to the Plancherel/modular row
        p_j = |S_0j|^2,   S_0j = sqrt(2/K) sin((2j+1) pi / K),   K = k+2.
  Self-normalized: sum_j |S_0j|^2 = (S S^dagger)_00 = 1 (no ad hoc division).
  The flux's quadratic invariant is Q_j = C_2(j)/[k(k+2)], C_2(j)=j(j+1).
  Abundance:  Omega_chi h^2 = 1.62 * Q_CS,   Q_CS = sum_j p_j Q_j.

Checks that reproduce the appendix numbers exactly:
  (A) modular/Plancherel measure  -> Q_CS = 0.0729566 -> Omega_chi h^2 = 0.1182  (FORCED, -1.5 sigma)
  (B) linear-d_j (Cardy) rival     -> Q     = 0.0767    -> Omega_chi h^2 = 0.1243  (EXCLUDED: needs a boundary)
  (C) collapse identity d_j^2/D^2 == |S_0j|^2 (the squared quantum-dim measure IS Plancherel)
  (D) edge-lock C_2(k/2)/[k(k+2)] = 1/4 exactly.

Uniqueness gates (thm:chi-measure-unique):
  (E) T^a S T^b PRESENTATION INVARIANCE: every genus-1 gluing of S^3 is (up to phase)
      T^a S T^b; T is diagonal-unitary, so the Born weights of T^a S T^b |0> are EXACTLY
      |S_0j|^2 for ALL (a,b) — the measure does not depend on the Heegaard presentation.
      Gate: max deviation over 20 random (a,b) pairs < 1e-12.
  (F) TREFOIL-CUT rival (knotted preparation): p_j prop S_0j^2 |J_j(3_1; q)|^2 with the
      colored Jones of the trefoil at q = e^{2 pi i/K} (Habiro cyclotomic form, self-
      validated at N=2 against the classical Jones polynomial V(3_1)=t+t^3-t^4).
      Gives Q_trefoil = 0.0737 != 0.0730 (Omega_chi h^2 = 0.1194): a DIFFERENT number,
      excluded by the STATED two-solid-tori (unknotted) preparation, chi_matter:677.
  (G) MERIDIAN-dual rival: p_j = delta_j0 => Q = C_2(0)/[k(k+2)] = 0 => Omega_chi h^2 = 0
      EXACTLY (zero dark matter; data-excluded => selection-grade kill of the rival basis).

Honesty-guard: the SELECTION of (A) over (B) is conditional on DFD's no-boundary
Wheeler-DeWitt postulate (a load-bearing DFD assumption, not an independent fact);
GRANTING it, |S_0j|^2 is forced and the linear d_j is excluded. We report numbers,
not adjectives.
"""
import math
import cmath
import random

k = 60
K = k + 2                       # = 62
PLANCK_OMEGA = 0.1200           # Planck 2018 omega_cdm target
CONV = 1.62                     # Omega_chi h^2 = CONV * Q_CS  (appendix_chi_matter)

# integrable reps j = 0, 1/2, ..., k/2  ->  twoj = 0,1,...,k  (k+1 = 61 sectors)
twoj = list(range(0, k + 1))
jvals = [t / 2.0 for t in twoj]

def S0(j):
    return math.sqrt(2.0 / K) * math.sin((2.0 * j + 1.0) * math.pi / K)

S0j = [S0(j) for j in jvals]
p_plinng = [s * s for s in S0j]                 # |S_0j|^2  (Plancherel/Born)
C2 = [j * (j + 1.0) for j in jvals]
Q_j = [c / (k * (k + 2.0)) for c in C2]          # C_2(j)/[k(k+2)]

# --- normalization sanity ---
norm = sum(p_plinng)

# --- (A) modular/Plancherel measure ---
Q_CS = sum(p * q for p, q in zip(p_plinng, Q_j))
Omega_A = CONV * Q_CS

# --- (B) linear-d_j (Cardy) rival ---
S00 = S0j[0]
d_j = [s / S00 for s in S0j]                      # quantum dimensions
p_cardy = [d / sum(d_j) for d in d_j]
Q_cardy = sum(p * q for p, q in zip(p_cardy, Q_j))
Omega_B = CONV * Q_cardy

# --- (C) collapse identity: d_j^2 / D^2 == |S_0j|^2 ---
D2 = sum(d * d for d in d_j)                      # = 1/S00^2
collapse = [ (d * d) / D2 for d in d_j ]
max_collapse_err = max(abs(c - p) for c, p in zip(collapse, p_plinng))
D2_vs_inv_S00sq = abs(D2 - 1.0 / (S00 * S00))

# --- (D) edge-lock ---
edge = C2[-1] / (k * (k + 2.0))                   # j = k/2

# =====================================================================
# Uniqueness gates (E), (F), (G)
# =====================================================================

# --- (E) T^a S T^b presentation invariance ---
# Full SU(2)_k modular data on the (k+1)-dim space (index by twoj = 0..k):
#   S_{ab} = sqrt(2/K) sin((a+1)(b+1) pi / K)          (a,b = twoj)
#   T_a    = exp(2 pi i (h_a - c/24)),  h_a = (a/2)(a/2+1)/K,  c = 3k/K
# T is DIAGONAL-unitary, so |<j| T^a S T^b |0>|^2 = |S_{j0}|^2 identically;
# this is the numerical content of the T^aST^b phase lemma (every genus-1
# gluing of S^3 gives the SAME Born weights => presentation-independence).
Smat = [[math.sqrt(2.0 / K) * math.sin((a + 1.0) * (b + 1.0) * math.pi / K)
         for b in twoj] for a in twoj]
c_central = 3.0 * k / K
Tdiag = [cmath.exp(2j * math.pi * ((a / 2.0) * (a / 2.0 + 1.0) / K - c_central / 24.0))
         for a in twoj]

random.seed(1234)                                  # reproducible
max_TST_dev = 0.0
for _ in range(20):
    a_pow = random.randint(1, 1000)
    b_pow = random.randint(1, 1000)
    # |0> = e_0 ;  T^b |0> = T_0^b e_0 ;  (S e_0)_j = S_{j0} ;  T^a phases rows
    v = [ (Tdiag[j] ** a_pow) * Smat[j][0] * (Tdiag[0] ** b_pow) for j in range(len(twoj)) ]
    dev = max(abs(abs(vj) ** 2 - p) for vj, p in zip(v, p_plinng))
    max_TST_dev = max(max_TST_dev, dev)

# --- (F) trefoil-cut rival: p_j prop S_0j^2 |J_j(trefoil; q)|^2 ---
# Colored Jones of the trefoil, N-dim color, Habiro cyclotomic expansion
# (Habiro 2002 / Masbaum AGT 2003), normalized J_1 = 1:
#   J_N(3_1; q) = sum_{n>=0} q^n (q^{1-N}; q)_n (q^{1+N}; q)_n   (terms vanish n>=N)
# Chirality only sends q -> 1/q = conj(q) at a root of unity, so |J|^2 is
# chirality-independent. Self-validation: N=2 must equal the classical Jones
# polynomial of the (right-handed) trefoil V(t) = t + t^3 - t^4.
def J_trefoil(N, q):
    tot = 0.0 + 0.0j
    prod = 1.0 + 0.0j
    for n in range(N):                             # (q^{1±N};q)_n kills n >= N
        tot += (q ** n) * prod
        prod *= (1.0 - q ** (1 - N + n)) * (1.0 - q ** (1 + N + n))
    return tot

t_test = complex(0.3, 0.4)                         # generic point, exact identity
jones_selfcheck = abs(J_trefoil(2, t_test) - (t_test + t_test ** 3 - t_test ** 4))

q_root = cmath.exp(2j * math.pi / K)               # SU(2)_k root of unity, K = k+2
w_tre = [p * abs(J_trefoil(t + 1, q_root)) ** 2 for p, t in zip(p_plinng, twoj)]
Q_tre = sum(w * q for w, q in zip(w_tre, Q_j)) / sum(w_tre)
Omega_tre = CONV * Q_tre

# --- (G) meridian-dual rival: p_j = delta_j0 ---
# Born weight of the SAME state read in the meridian basis: all weight on j=0,
# and C_2(0) = 0, so Q = 0 and Omega_chi h^2 = 0 EXACTLY (zero dark matter).
p_meridian = [1.0] + [0.0] * k
Q_meridian = sum(p * q for p, q in zip(p_meridian, Q_j))
Omega_meridian = CONV * Q_meridian

print("="*72)
print(" chi relic abundance from the closed-S^3 modular-vacuum measure")
print("="*72)
print(f"level k = {k},  K = k+2 = {K},  integrable sectors = {len(jvals)} (j=0..{k//2})")
print(f"normalization  sum_j |S_0j|^2            = {norm:.10f}   (target 1)")
print()
print("(A) MODULAR / PLANCKEREL measure p_j = |S_0j|^2  [no-boundary closed-S^3 vacuum]")
print(f"    Q_CS = sum_j |S_0j|^2 C_2(j)/[k(k+2)] = {Q_CS:.7f}   (appendix 0.0729566)")
print(f"    Omega_chi h^2 = 1.62 * Q_CS           = {Omega_A:.4f}     (target {PLANCK_OMEGA}, FORCED)")
print()
print("(B) LINEAR-d_j (Cardy boundary) rival -- EXCLUDED (no internal boundary on closed S^3)")
print(f"    Q_Cardy                               = {Q_cardy:.7f}   (appendix 0.0767)")
print(f"    Omega_chi h^2 (Cardy)                 = {Omega_B:.4f}     (the +3.6 sigma overshoot)")
print()
print("(C) collapse identity  d_j^2/D^2 == |S_0j|^2  (squared quantum-dim = Plancherel)")
print(f"    max |d_j^2/D^2 - |S_0j|^2|            = {max_collapse_err:.2e}   (should be ~0)")
print(f"    |D^2 - 1/S_00^2|                      = {D2_vs_inv_S00sq:.2e}   (should be ~0)")
print()
print("(D) edge-lock  C_2(k/2)/[k(k+2)]")
print(f"    value                                 = {edge:.10f}   (exact 1/4 = 0.25)")
print()
print("(E) T^a S T^b PRESENTATION INVARIANCE (phase lemma, 20 random (a,b))")
print(f"    max | |<j|T^a S T^b|0>|^2 - |S_0j|^2 |= {max_TST_dev:.2e}   (gate < 1e-12)")
print()
print("(F) TREFOIL-CUT rival  p_j prop S_0j^2 |J_j(3_1; q)|^2,  q = e^(2 pi i/62)")
print(f"    Habiro N=2 self-check vs V(3_1)=t+t^3-t^4 = {jones_selfcheck:.2e}   (should be ~0)")
print(f"    Q_trefoil                             = {Q_tre:.7f}   (reference 0.0737, != Q_CS {Q_CS:.4f})")
print(f"    Omega_chi h^2 (trefoil)               = {Omega_tre:.4f}     (reference 0.1194; EXCLUDED by the")
print( "                                            STATED unknotted two-solid-tori preparation)")
print()
print("(G) MERIDIAN-dual rival  p_j = delta_j0")
print(f"    Q_meridian                            = {Q_meridian:.1f}         (exact 0)")
print(f"    Omega_chi h^2 (meridian)              = {Omega_meridian:.1f}         (ZERO dark matter, data-excluded)")
print()

# --- pass/fail gates (for reproduce.py / golden_master_gates.py) ---
ok_A    = abs(Q_CS - 0.0729566) < 5e-7 and abs(Omega_A - 0.1182) < 5e-4
ok_B    = abs(Q_cardy - 0.0767) < 5e-4
ok_C    = max_collapse_err < 1e-12 and D2_vs_inv_S00sq < 1e-9   # 61-term float roundoff
ok_D    = abs(edge - 0.25) < 1e-12
ok_norm = abs(norm - 1.0) < 1e-12
ok_E    = max_TST_dev < 1e-12
ok_F    = (jones_selfcheck < 1e-12 and abs(Q_tre - 0.0737) < 1e-3
           and abs(Q_tre - Q_CS) > 5e-4 and abs(Omega_tre - 0.1194) < 1.5e-3)
ok_G    = (Q_meridian == 0.0) and (Omega_meridian == 0.0)
allok = ok_A and ok_B and ok_C and ok_D and ok_norm and ok_E and ok_F and ok_G
print(f"CHECKS  A:{ok_A}  B:{ok_B}  C:{ok_C}  D:{ok_D}  norm:{ok_norm}  E:{ok_E}  F:{ok_F}  G:{ok_G}")
print(f"RESULT: {'ALL PASS' if allok else 'FAIL'} -- modular measure forces Omega_chi h^2=0.1182;"
      f" linear-d_j (0.1243) excluded by closed-S^3 (no boundary);")
print( "        presentation-invariance exact (E); trefoil rival = different number, excluded by stated")
print( "        preparation (F); meridian rival = zero dark matter, data-excluded (G).")
print("Honesty-guard: selection conditional on DFD no-boundary WdW postulate (P1-P4 of thm:chi-measure-unique);")
print("           numbers reported, not adjectives.")
import sys
sys.exit(0 if allok else 1)
