# 10,000 simulated experiments per scenario, 8 mice per group, same model as the worked example (model.py).
# Python 3.13, NumPy, SciPy. Seed fixed before running.
import numpy as np
from scipy import stats
from model import simulate, DAYS
rng = np.random.default_rng(20261111); R = 10_000; n = 8

def slopes(v):                                         # per-mouse least-squares slope of log volume on day
    lv = np.log(v); ok = ~np.isnan(lv); x = np.where(ok, DAYS, np.nan)
    xm = np.nanmean(x, 1, keepdims=True); ym = np.nanmean(lv, 1, keepdims=True)
    return np.nansum((x - xm) * (lv - ym), 1) / np.nansum((x - xm) ** 2, 1)

def run(rate_v, rate_d):
    sig_slope = sig_last = sig_any = no_test = 0; tc = []
    for _ in range(R):
        v, d = simulate(rng, n, rate_v), simulate(rng, n, rate_d)
        sig_slope += stats.ttest_ind(slopes(d), slopes(v), equal_var=False).pvalue < 0.05
        ps = []
        for j in range(len(DAYS)):
            a, b = np.log(v[:, j]), np.log(d[:, j]); a, b = a[~np.isnan(a)], b[~np.isnan(b)]
            ps.append(stats.ttest_ind(b, a, equal_var=False).pvalue if min(len(a), len(b)) >= 2 else np.nan)
        if np.isnan(ps[-1]):
            no_test += 1
        else:
            sig_last += ps[-1] < 0.05
            a, b = np.log(v[:, -1]), np.log(d[:, -1]); tc.append(np.exp(np.nanmean(b) - np.nanmean(a)))
        sig_any += np.nanmin(ps[1:]) < 0.05                  # any day after day 0
    return sig_slope / R, sig_last / R, sig_any / R, no_test / R, np.median(tc)

for label, rv, rd in [("No drug effect (both 0.15 per day)", 0.15, 0.15), ("Drug slows growth (0.15 vs 0.10 per day)", 0.15, 0.10)]:
    s, l, a, nt, tc = run(rv, rd)
    print(f"{label}: growth-rate comparison p<0.05 {s:.1%}; day-21 survivors p<0.05 {l:.1%} "
          f"(no test possible in {nt:.1%}); any single day p<0.05 {a:.1%}; median day-21 T/C among survivors {tc:.2f}")
print(f"True day-21 T/C without culling for the drug scenario: {np.exp((0.10 - 0.15) * 21):.2f}")