# 10,000 simulated studies with the mepsc.csv design (5 WT and 5 KO mice, 4 neurons per mouse, 200 events per neuron).
# Python 3.13, NumPy, SciPy. Seed fixed before running.
import numpy as np
from scipy import stats
from model import study
rng = np.random.default_rng(20261112); R = 10_000

def one(ratio, mice):
    rows = study(rng, ko_ratio=ratio, mice=mice)
    ko = np.array([r[0] == "KO" for r in rows]); amp = [r[3] for r in rows]
    ev_ko = np.concatenate([a for a, k in zip(amp, ko) if k]); ev_wt = np.concatenate([a for a, k in zip(amp, ko) if not k])
    p_ks = stats.ks_2samp(ev_ko, ev_wt).pvalue
    med = np.log([np.median(a) for a in amp]); p_cell = stats.ttest_ind(med[ko], med[~ko], equal_var=False).pvalue
    mouse = np.array([r[1] for r in rows]); um = list(dict.fromkeys(mouse))
    mm = np.array([med[mouse == m].mean() for m in um]); kmm = np.array([m.startswith("KO") for m in um])
    p_mouse = stats.ttest_ind(mm[kmm], mm[~kmm], equal_var=False).pvalue
    return p_ks < 0.05, p_cell < 0.05, p_mouse < 0.05
for lab, ratio in (("No genotype difference", 1.0), ("KO amplitude 15% smaller", 0.85)):
    for mice in (5, 8):
        r = np.array([one(ratio, mice) for _ in range(R)])
        print(f"{lab}, {mice} mice per group: p<0.05 with all events pooled (KS) {r[:, 0].mean():.1%}; neurons as n {r[:, 1].mean():.1%}; mice as n {r[:, 2].mean():.1%}")