# 10,000 simulated studies with the wound.csv design (5 experiments, 3 wells per condition, 3 fields per well, images at
# 0, 6, 12 and 24 h). Python 3.13, NumPy, SciPy. Seed fixed before running.
import numpy as np, warnings
warnings.filterwarnings("ignore", category=RuntimeWarning)
from scipy import stats
from model import experiment, TIMES, WELLS, FIELDS
rng = np.random.default_rng(20261112); R = 10_000

def one(ratio, sd_cond=0.08):
    sp_w, cl_w, cond_w, sp_e, cl_e = [], [], [], [], []
    for e in range(5):
        rows = experiment(rng, drug_ratio=ratio, sd_cond=sd_cond)
        y = np.array([r[4] for r in rows]).reshape(2, WELLS, FIELDS, len(TIMES))
        ok = y > 20                                                     # linear phase only: points before closure
        n = ok.sum(-1); tm = (ok * TIMES).sum(-1) / n; ym = (ok * y).sum(-1) / n
        speed = -((ok * (TIMES - tm[..., None]) * (y - ym[..., None])).sum(-1) / (ok * (TIMES - tm[..., None]) ** 2).sum(-1)) / 2; closure = 100 * (y[..., 0] - y[..., -1]) / y[..., 0]
        sw, cw = np.nanmean(speed, -1), closure.mean(-1)                # per well (fields with fewer than 2 open time points skipped)
        sp_w += list(sw.ravel()); cl_w += list(cw.ravel()); cond_w += [0] * WELLS + [1] * WELLS
        sp_e.append(np.log(sw[1].mean() / sw[0].mean())); cl_e.append(cw[1].mean() - cw[0].mean())
    sp_w, cond_w = np.log(sp_w), np.array(cond_w, bool)
    return (stats.ttest_ind(sp_w[cond_w], sp_w[~cond_w], equal_var=False).pvalue < 0.05,
            stats.ttest_1samp(sp_e, 0).pvalue < 0.05, stats.ttest_1samp(cl_e, 0).pvalue < 0.05)
for lab, ratio, sc in (("Drug has no effect", 1.0, 0.08), ("Drug has no effect, conditions vary more between experiments (SD 0.15)", 1.0, 0.15), ("Drug slows migration by 20%", 0.8, 0.08)):
    r = np.array([one(ratio, sc) for _ in range(R)])
    print(f"{lab}: p<0.05 for migration speed with wells as n {r[:, 0].mean():.1%}; migration speed with experiments as n {r[:, 1].mean():.1%}; "
          f"24-h percent closure with experiments as n {r[:, 2].mean():.1%}")