# 100,000 simulated two-group experiments per scenario (normal outcome, SD 1; Student's t-test, two-sided).
# Python 3.13, NumPy, SciPy. Seed fixed before running.
import numpy as np
from scipy import stats
rng = np.random.default_rng(20261111); R = 100_000; NMAX = 20

def p_at(c, t, n): return stats.ttest_ind(t[:, :n], c[:, :n], axis=1).pvalue

def p_block(c, t):                                    # y ~ group + cohort; cohorts of 8 + 8 and 4 + 4 (orthogonal)
    y = np.hstack([c[:, :12], t[:, :12]]); g = np.r_[np.zeros(12), np.ones(12)]; k = np.r_[np.zeros(8), np.ones(4)]
    k = np.r_[k, k]
    gm = y.mean(1, keepdims=True)
    ss_tot = ((y - gm) ** 2).sum(1)
    diff = y[:, g == 1].mean(1) - y[:, g == 0].mean(1)
    ss_g = 12 * ((y[:, g == 1].mean(1) - gm[:, 0]) ** 2 + (y[:, g == 0].mean(1) - gm[:, 0]) ** 2)
    ss_k = 16 * (y[:, k == 0].mean(1) - gm[:, 0]) ** 2 + 8 * (y[:, k == 1].mean(1) - gm[:, 0]) ** 2
    mse = (ss_tot - ss_g - ss_k) / 21
    return 2 * stats.t.sf(np.abs(diff) / np.sqrt(mse * (2 / 12)), 21)

def strategies(delta):
    c, t = rng.normal(0, 1, (R, NMAX)), rng.normal(delta, 1, (R, NMAX))
    p8, p12, p16, p20 = (p_at(c, t, n) for n in (8, 12, 16, 20))
    out = {}
    out["Fixed n = 8"] = (p8 < 0.05, np.full(R, 8))
    out["Fixed n = 16"] = (p16 < 0.05, np.full(R, 16))
    promising = (p8 >= 0.05) & (p8 < 0.10)
    out["n = 8; if 0.05 <= p < 0.10, add 4 and test again"] = ((p8 < 0.05) | (promising & (p12 < 0.05)), np.where(promising, 12, 8))
    out["The same, second test with cohort as a block"] = ((p8 < 0.05) | (promising & (p_block(c, t) < 0.05)), np.where(promising, 12, 8))
    sig = p8 < 0.05; n_used = np.full(R, 8)
    for n, p in ((12, p12), (16, p16), (20, p20)):
        go = ~sig; n_used = np.where(go, n, n_used); sig = sig | (go & (p < 0.05))
    out["n = 8; while p >= 0.05 add 4, up to 20"] = (sig, n_used)
    stop = p8 < 0.0294
    out["Planned two-stage design (8, then 16), Pocock p < 0.0294 at each look"] = (stop | (~stop & (p16 < 0.0294)), np.where(stop, 8, 16))
    return out

null, eff = strategies(0.0), strategies(1.0)
for k in null:
    print(f"{k}: false positives {null[k][0].mean():.1%}; power at d = 1 {eff[k][0].mean():.1%}; "
          f"mean mice per group at d = 1 {eff[k][1].mean():.1f}")