# 10,000 simulated studies with the events.csv design (3 staining days, 2 WT and 2 KO mice per day, 1,000 events per mouse).
# 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):
    rows = study(rng, ratio=ratio)
    ko = np.array([g == "KO" for _, _, g, _ in rows]); day = np.array([d for d, _, _, _ in rows])
    ev_ko = np.concatenate([f for _, _, g, f in rows if g == "KO"]); ev_wt = np.concatenate([f for _, _, g, f in rows if g == "WT"])
    p_events = stats.ttest_ind(ev_ko, ev_wt, equal_var=False).pvalue
    lm = np.log([np.median(f) for _, _, _, f in rows])
    p_mouse = stats.ttest_ind(lm[ko], lm[~ko], equal_var=False).pvalue
    X = np.column_stack([np.ones(len(lm)), ko, day == 2, day == 3]).astype(float)
    b, *_ = np.linalg.lstsq(X, lm, rcond=None); r = lm - X @ b; dfe = len(lm) - 4
    se = np.sqrt(r @ r / dfe * np.linalg.inv(X.T @ X)[1, 1]); p_block = 2 * stats.t.sf(abs(b[1] / se), dfe)
    return p_events < 0.05, p_mouse < 0.05, p_block < 0.05
for lab, ratio in (("No difference (ratio 1)", 1.0), ("KO 20% lower (ratio 0.8)", 0.8)):
    r = np.array([one(ratio) for _ in range(R)])
    print(f"{lab}: p<0.05 with events as n {r[:, 0].mean():.1%}; per-mouse median, ignoring day {r[:, 1].mean():.1%}; "
          f"per-mouse median, day as a block {r[:, 2].mean():.1%}")