#!/usr/bin/env python3
"""LENS-2 independent verification of the G8d null ladder (2026-07-02).
Independent seed, independent MC implementation, plus EXACT analytics
for the product-of-uniforms null so MC is cross-checked against closed form.
"""
import math, json
import numpy as np
from statistics import NormalDist

nd = NormalDist()
Z = lambda p: nd.inv_cdf(1 - p) if 0 < p < 1 else 0.0
rng = np.random.default_rng(987654321)  # different seed from G8d (20260702)

out = {}

# ------------------------------------------------------------------
# 0. Re-derive the flagship per-hit charges from corpus anchors
# ------------------------------------------------------------------
# s2w rational menu: reduced p/q in [0.15,0.35], window |3/13-0.23122|
meas, win = 0.23122, abs(3/13 - 0.23122)
def rmenu(qmax):
    tot = hit = 0
    for q in range(2, qmax+1):
        for p in range(1, q):
            if math.gcd(p, q) == 1 and 0.15 <= p/q <= 0.35:
                tot += 1
                hit += (abs(p/q - meas) <= win)
    return hit, tot
out['s2w_menus'] = {q: rmenu(q) for q in (13, 20, 30, 60)}
# conservative charge = worst ratio (qmax=13 -> 1/12); moderate = qmax=30 -> 1/56

# AU grammar scan reproduction (135 candidates, band [1.8,2.8] meV)
alpha = 1/137.035999084
E_P_meV = 1.220890e28 * 1e3
prefs = [1.0, 3/(8*math.pi), 1/(8*math.pi), 3/(4*math.pi), 1/(4*math.pi), 0.5, 2.0, math.pi, 1/math.pi]
au_hits = sum(1 for pv in prefs for n in range(50, 65)
              if 1.8 <= (pv**0.25)*(alpha**(n/4))*E_P_meV <= 2.8)
out['AU_scan'] = {'hits': au_hits, 'candidates': 9*15}

# flagship products
p_cons = (1/12) * (3/135) * 0.136
p_mod  = (1/56) * ((3/135)*0.074) * 0.136
out['flagship'] = {'cons': {'p': p_cons, 'Z': Z(p_cons)},
                   'mod':  {'p': p_mod,  'Z': Z(p_mod)}}
# 4-hit fair set (theta12 readmitted at ansatz-menu 1/14)
p_cons4 = p_cons * (1/14)
p_mod4  = p_mod  * (1/14)
out['flagship_plus_theta12'] = {'cons4': {'p': p_cons4, 'Z': Z(p_cons4)},
                                'mod4':  {'p': p_mod4,  'Z': Z(p_mod4)}}

# ------------------------------------------------------------------
# 1. EXACT null: P(product of ALL n iid U(0,1) <= t) = t * sum_{j<n} (-ln t)^j/j!
#    (pre-registered full-scorecard product test, misses charged at p=1)
# ------------------------------------------------------------------
def p_product_all(t, n):
    L = -math.log(t)
    return t * sum(L**j / math.factorial(j) for j in range(n))
out['full_scorecard_product_M8'] = {
    'cons (obs=cons3*theta12, misses@1)': p_product_all(p_cons4, 8),
    'mod  (obs=mod3*theta12, misses@1)':  p_product_all(p_mod4, 8)}

# ------------------------------------------------------------------
# 2. MC: order-statistics selection null (product of k smallest of M uniforms)
#    independent implementation (sort-based, vectorized)
# ------------------------------------------------------------------
def sel_p(t_obs, M, k, nsim=4_000_000):
    u = rng.random((nsim, M))
    u.sort(axis=1)
    stat = np.prod(u[:, :k], axis=1)
    return float(np.mean(stat <= t_obs))

sel = {}
for tag, t, k in [('cons3', p_cons, 3), ('mod3', p_mod, 3),
                  ('cons4', p_cons4, 4), ('mod4', p_mod4, 4)]:
    for M in (8, 10, 12, 20):
        if M >= k:
            p = sel_p(t, M, k)
            sel[f'{tag}_M{M}_k{k}'] = {'p': p, 'Z': Z(p)}
out['selection_null'] = sel

# ------------------------------------------------------------------
# 3. Poisson count cross-check (hit-capable comparisons only)
#    lam = sum of per-comparison hit probabilities under the null
# ------------------------------------------------------------------
from math import exp
def pois_ge(k, lam):
    return 1 - sum(exp(-lam)*lam**j/math.factorial(j) for j in range(k))
# moderate: 4 real windows + (M-4) pool comparisons at generic menu-charge 0.05
for M in (10, 12, 20):
    lam = (1/56) + (3/135)*0.074 + 0.136 + (1/14) + (M-4)*0.05
    out[f'poisson_M{M}_mod_P(>=4)'] = pois_ge(4, lam)

# ------------------------------------------------------------------
# 4. Corpus App CE floor-4 arithmetic re-check + honest-row repair
# ------------------------------------------------------------------
ce = {'s2w': 0.004, 'a0': 0.005, 'theta12': 0.074, 'sin2beta': 0.031}
p_floor_ce = math.prod(ce.values())
out['appCE_floor4_as_written'] = {'p': p_floor_ce, 'Z': Z(p_floor_ce)}
# repair: sin2beta -> 1.0 (function of 4 postulate-selected integers, AO anti-thm);
# a0 -> 0.20 honest (empirical a0 known only to ~10%: 2*0.10 window on O(1))
ce_fix = {'s2w': 0.004, 'a0': 0.20, 'theta12': 0.074}
p_floor_fix = math.prod(ce_fix.values())
out['appCE_floor_repaired_own_null'] = {'p': p_floor_fix, 'Z': Z(p_floor_fix)}

# maximal-6 as written
ce6 = [0.004, 0.005, 0.004, 0.009, 0.074, 0.031]
p6 = math.prod(ce6)
out['appCE_max6_as_written'] = {'p': p6, 'Z': Z(p6)}

print(json.dumps(out, indent=1, default=str))
with open('verify3x/lens2_results.json', 'w') as f:
    json.dump(out, f, indent=1, default=str)
