# How often each analysis gives p < 0.05: 5,000 simulated studies of 6 + 6 mice with 1-3 slices each, with no knockout
# effect and with LTP potentiation reduced by 40% (about 150% to 130% of baseline). Python 3.13, NumPy, SciPy.
# Seed 20261112. Column D (a test at every minute with mice as n) and the 12 + 12 mice scenarios were added after
# review; they come last or use no random numbers, so the earlier rows are unchanged.
import numpy as np
from scipy import stats
from model import study, MINUTES
rng = np.random.default_rng(20261112); R = 5_000
post_idx = MINUTES > 0; late_idx = (MINUTES >= 51) & (MINUTES <= 60)

def one(ratio, mice=6):
    rows = study(rng, ratio, mice=mice); n = len(MINUTES)
    slope = np.array([r[4] for r in rows]).reshape(-1, n); ko = np.array([r[0] == "knockout" for r in rows[::n]]); mouse = np.array([r[1] for r in rows[::n]])
    pct = 100 * slope / slope[:, ~post_idx].mean(axis=1, keepdims=True)
    ltp = pct[:, late_idx].mean(axis=1)
    pa = stats.ttest_ind(ltp[ko], ltp[~ko], equal_var=False).pvalue
    pb = stats.ttest_ind(pct[ko][:, post_idx], pct[~ko][:, post_idx], axis=0, equal_var=False).pvalue.min()
    mice, inv = np.unique(mouse, return_inverse=True); mltp = np.bincount(inv, ltp) / np.bincount(inv); mko = np.bincount(inv, ko) > 0
    pc = stats.ttest_ind(mltp[mko], mltp[~mko], equal_var=False).pvalue
    mpct = np.array([pct[inv == i].mean(axis=0) for i in range(len(mice))])
    pd_ = stats.ttest_ind(mpct[mko][:, post_idx], mpct[~mko][:, post_idx], axis=0, equal_var=False).pvalue.min()
    return np.array([pa, pb, pc, pd_]) < 0.05
for lab, ratio, mice in (("No effect", 1.0, 6), ("Potentiation reduced by 40%", 0.6, 6),
                         ("No effect, 12 + 12 mice", 1.0, 12), ("Potentiation reduced by 40%, 12 + 12 mice", 0.6, 12)):
    r = np.array([one(ratio, mice) for _ in range(R)]).mean(axis=0)
    print(f"{lab}: A slices as n {r[0]:.1%}; B any of 60 minutes, slices as n {r[1]:.1%}; C mice as n {r[2]:.1%}; "
          f"D any of 60 minutes, mice as n {r[3]:.1%}")