# 10,000 simulated studies: 3 independent experiments x 2 lines x 3 wells, same model as the worked example.
# Python 3.13, NumPy, SciPy. Seed fixed before running.
import numpy as np
from scipy import stats
from model import experiment, TIMES
rng = np.random.default_rng(20261111); R = 10_000
def rates(x, upto):                                  # per-experiment rate from mean log counts, times up to `upto`
    m = np.log(x).mean(0); k = TIMES <= upto
    return np.polyfit(TIMES[k], m[k], 1)[0]
for lab, rk in (("No difference (both 0.035 per hour)", 0.035), ("KO slower (0.035 vs 0.030 per hour)", 0.030)):
    hits = np.zeros(2); dt = {48: [], 72: [], 96: []}
    for _ in range(R):
        ew = [experiment(rng, 0.035) for _ in range(3)]; ek = [experiment(rng, rk) for _ in range(3)]
        rw = np.array([rates(x, 72) for x in ew]); rkk = np.array([rates(x, 72) for x in ek])
        hits[0] += stats.ttest_rel(rkk, rw).pvalue < 0.05
        hits[1] += stats.ttest_ind(np.concatenate([x[:, -1] for x in ek]), np.concatenate([x[:, -1] for x in ew]), equal_var=False).pvalue < 0.05
        for u in dt: dt[u].append(np.log(2) / np.mean([rates(x, u) for x in ew]))
    print(f"{lab}: p<0.05 with the paired test on experiment growth rates {hits[0] / R:.1%}; with 96-h counts, wells as n {hits[1] / R:.1%}")
    if rk == 0.035:
        print("  WT doubling time (true 19.8 h), median estimate fitting up to 48 h: {:.1f} h; up to 72 h: {:.1f} h; up to 96 h: {:.1f} h".format(*[np.median(dt[u]) for u in (48, 72, 96)]))