#!/usr/bin/env python3
"""E5b: moderate tier (still conservative, but per-hit windows not
resolution-capped) + Poisson expected-hits cross-check. Brackets the E5
verdict so it does not rest on a single maximally-stacked corner."""
import json, math
import numpy as np
from statistics import NormalDist
Z = lambda p: NormalDist().inv_cdf(1-p) if 0 < p < 1 else 0.0
rng = np.random.default_rng(20260702)
out = {}

# MODERATE charged p-values (each still conservative, derivation in comments):
# s2w: rational menu q<=30 (natural 'simple fraction' universe): 1/56
# theta12: ansatz menu 1/14
# H0/Lam/a0 block: AU grammar 3/135 x SH0ES refinement (1.9/51 span) x2 dataset choice
# Omega_chi: kept at E4 hostile 0.136 (NOT softened)
p_mod = {'s2w': 1/56, 'theta12': 1/14, 'H0LamA0': (3/135)*(2*0.95/51)*2, 'Omega_chi': 0.136}
T_mod = math.prod(p_mod.values())
out['moderate_survivors'] = p_mod
out['T_mod'] = T_mod
out['Z_mod_no_selection'] = Z(T_mod)

def mc(M, k, T, nsim=6_000_000, chunk=500_000):
    hits = 0
    done = 0
    while done < nsim:
        n = min(chunk, nsim - done)
        u = rng.random((n, M))
        u.sort(axis=1)
        pr = np.prod(u[:, :k], axis=1)
        hits += int((pr <= T).sum())
        done += n
    return max(hits, 1) / nsim

res = {}
for M in (12, 20, 30):
    res[f'M={M}'] = {'p_sel_mod_k4': mc(M, 4, T_mod)}
    res[f'M={M}']['p_sel_mod_k4_kscan_x3'] = min(1.0, 3*res[f'M={M}']['p_sel_mod_k4'])
out['selection_null_moderate'] = res

# Poisson cross-check: expected number of >=0.136-grade hits among M
# comparisons if every comparison had hit-prob equal to the LARGEST
# surviving window (hostile to DFD), lambda = M x 0.136; and a moderate
# lambda = sum of assigned charged p's over an enumerated 20-comparison
# ledger (counted hits at their charges, all others granted 0.05 each).
from math import exp, factorial
def pois_ge(lam, k):
    return 1 - sum(exp(-lam)*lam**i/factorial(i) for i in range(k))
lam_host = 20*0.136
lam_mod = sum(p_mod.values()) + 16*0.05
out['poisson_check'] = {
 'lam_hostile(all 20 at 0.136)': lam_host, 'P(>=4)_hostile': pois_ge(lam_host, 4),
 'lam_moderate(sum charges +16x0.05)': lam_mod, 'P(>=4)_moderate': pois_ge(lam_mod, 4),
 'note': 'counts-only check; ignores HOW small the survivor windows are, so it understates evidence; order-stat MC is the operative number'}

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