#!/usr/bin/env python3
"""
Monte Carlo v2 for 'The Graveyard Is the Product' — full suite.

Experiments:
  1. Pure null (zero alpha): false-acceptance rates of research behaviors.
  2. Era-locked alternative (real in-sample drift confined to one temporal
     group): generalization-failure rates.
  3. Power: temporally stable true alpha at annualized SR 0.5/0.75/1.0.
  4. Correlated variants (rho = 0.8): null size + power at SR 0.5.

Rules evaluated (candidate = best-of-K by full-sample t unless stated):
  naive       : accept if candidate t > 2
  reuse       : select best-of-K ON a 30% holdout, accept if its t > 2
  clean       : select best-of-K on 70% train, accept if untouched-holdout t > 2
  folds       : >= 80% of the 15 CPCV split test-means positive
  defl        : candidate t > t0 + 1.645, t0 = BLdP closed form for K trials
  protocol    : folds AND defl
  lbeo        : leave-best-era-out — drop the candidate's best temporal group,
                recompute t on the remaining 5/6 sample, require
                t_lbeo > (t0 + 1.645) * sqrt(5/6)
  lbeoproto   : defl AND lbeo
Analytic cross-checks printed where closed forms exist.
"""
import numpy as np
from itertools import combinations
from math import sqrt, log, e as EULER_E

# --- normal helpers (no scipy dependency) -----------------------------------
def norm_cdf(x):
    from math import erf
    return 0.5 * (1 + erf(x / sqrt(2)))

def norm_ppf(p, lo=-10, hi=10):
    for _ in range(80):
        mid = (lo + hi) / 2
        if norm_cdf(mid) < p: lo = mid
        else: hi = mid
    return (lo + hi) / 2

GAMMA = 0.5772156649
def bldp_t0(K):
    return (1 - GAMMA) * norm_ppf(1 - 1 / K) + GAMMA * norm_ppf(1 - 1 / (K * EULER_E))

# --- config -----------------------------------------------------------------
rng = np.random.default_rng(20260810)
T, K, SIG = 2016, 20, 0.01
G = 6
bounds = np.linspace(0, T, G + 1).astype(int)
gslice = [slice(bounds[i], bounds[i + 1]) for i in range(G)]
glen = np.array([bounds[i + 1] - bounds[i] for i in range(G)])
SPLITS = list(combinations(range(G), 2))
CUT = int(T * 0.7)

T0 = bldp_t0(K)
HURDLE = T0 + 1.645
LBEO_HURDLE = HURDLE * sqrt(5 / 6)

def run(n_fam, drift=None, sr_ann=0.0, rho=0.0, batch=500):
    """drift='era' for era-locked; sr_ann>0 for stable alpha; rho common factor."""
    mu = sr_ann * SIG / sqrt(252)
    counts = dict(naive=0, reuse=0, clean=0, folds=0, defl=0,
                  protocol=0, lbeo=0, lbeoproto=0)
    done = 0
    while done < n_fam:
        b = min(batch, n_fam - done); done += b
        if rho > 0:
            f = rng.normal(0, SIG, size=(b, 1, T))
            eps = rng.normal(0, SIG, size=(b, K, T))
            R = sqrt(rho) * f + sqrt(1 - rho) * eps
        else:
            R = rng.normal(0, SIG, size=(b, K, T))
        if mu: R += mu
        if drift == "era":
            MU = 2.5 * sqrt(T) / (T / G) * SIG
            g = rng.integers(0, G, size=(b, K))
            for i in range(b):
                for k in range(K):
                    R[i, k, gslice[g[i, k]]] += MU
        m = R.mean(2); s = R.std(2, ddof=1)
        t_full = m / (s / sqrt(T))
        sel = t_full.argmax(1)
        ar = np.arange(b)
        tc = t_full[ar, sel]                              # candidate t
        counts["naive"] += int((t_full.max(1) > 2).sum())
        # reuse: select on 30% tail, accept same t
        tail = R[:, :, CUT:]
        t_tail = tail.mean(2) / (tail.std(2, ddof=1) / sqrt(T - CUT))
        counts["reuse"] += int((t_tail.max(1) > 2).sum())
        # clean: select on train, evaluate holdout once
        tr = R[:, :, :CUT]
        t_tr = tr.mean(2) / (tr.std(2, ddof=1) / sqrt(CUT))
        sel_tr = t_tr.argmax(1)
        hold = R[ar, sel_tr, CUT:]
        t_hold = hold.mean(1) / (hold.std(1, ddof=1) / sqrt(T - CUT))
        counts["clean"] += int((t_hold > 2).sum())
        # folds on candidate: group means
        gm = np.stack([R[ar, sel, sl].mean(1) for sl in gslice], 1)  # b x G
        ok = np.zeros(b, dtype=int)
        for (i, j) in SPLITS:
            w = (gm[:, i] * glen[i] + gm[:, j] * glen[j]) / (glen[i] + glen[j])
            ok += (w > 0)
        folds_pass = ok >= 12                          # >=80% of 15
        counts["folds"] += int(folds_pass.sum())
        defl_pass = tc > HURDLE
        counts["defl"] += int(defl_pass.sum())
        counts["protocol"] += int((folds_pass & defl_pass).sum())
        # LBEO on candidate
        best_g = gm.argmax(1)
        lbeo_pass = np.zeros(b, dtype=bool)
        for i in range(b):
            keep = np.ones(T, dtype=bool)
            keep[gslice[best_g[i]]] = False
            r = R[i, sel[i], keep]
            t_l = r.mean() / (r.std(ddof=1) / sqrt(keep.sum()))
            lbeo_pass[i] = t_l > LBEO_HURDLE
        counts["lbeo"] += int(lbeo_pass.sum())
        counts["lbeoproto"] += int((lbeo_pass & defl_pass).sum())
    return {k: v / n_fam for k, v in counts.items()}

def pct(x): return f"{100*x:6.2f}%"

print(f"K={K}, T={T}, BLdP t0={T0:.4f}, hurdle={HURDLE:.4f}, LBEO hurdle={LBEO_HURDLE:.4f}")
print(f"analytic: naive null = {100*(1-norm_cdf(2)**K):.2f}%  "
      f"clean null = {100*(1-norm_cdf(2)):.2f}%  "
      f"defl null = {100*(1-norm_cdf(HURDLE)**K):.2f}%")
for Kx in (5, 20, 100):
    print(f"  K-sensitivity naive: K={Kx}: {100*(1-norm_cdf(2)**Kx):.2f}%")

N1 = 20000
r1 = run(N1)
print(f"\nEXP1 pure null (n={N1}):")
for k, v in r1.items(): print(f"  {k:10s} {pct(v)}")

N2 = 10000
r2 = run(N2, drift="era")
print(f"\nEXP2 era-locked alternative (n={N2}) — generalization-failure rates:")
for k, v in r2.items(): print(f"  {k:10s} {pct(v)}")

N3 = 10000
for sr in (0.5, 0.75, 1.0):
    r3 = run(N3, sr_ann=sr)
    print(f"\nEXP3 power, stable alpha SR={sr} (t≈{sr*sqrt(8):.2f}, n={N3}):")
    for k, v in r3.items(): print(f"  {k:10s} {pct(v)}")

N4 = 10000
r4n = run(N4, rho=0.8)
r4p = run(N4, rho=0.8, sr_ann=0.5)
print(f"\nEXP4 correlated variants rho=0.8 (n={N4}):")
print("  null size:")
for k in ("naive", "defl", "lbeoproto"): print(f"    {k:10s} {pct(r4n[k])}")
print("  power at SR=0.5:")
for k in ("naive", "clean", "defl", "protocol", "lbeoproto"): print(f"    {k:10s} {pct(r4p[k])}")