# 10,000 simulated studies with the viability.csv design (5 experiments, one plate each, 3 wells per condition).
# Python 3.13, NumPy, SciPy. Seed fixed before running.
import numpy as np
from scipy import stats
from model import experiment
rng = np.random.default_rng(20261112); R = 10_000

def study(sf_ab):
    lb, lh, plate1 = [], [], None
    for e in range(5):
        rows = experiment(rng, sf_ab=sf_ab)
        sig = {g: np.array([r[2] for r in rows if r[0] == g]) for g in ("control", "A", "B", "A+B")}
        c = sig["control"].mean(); sf = {g: sig[g] / c for g in ("A", "B", "A+B")}; m = {g: sf[g].mean() for g in sf}
        lb.append(np.log(m["A+B"] / (m["A"] * m["B"]))); lh.append(np.log(m["A+B"] / min(m["A"], m["B"])))
        if e == 0: plate1 = stats.ttest_1samp(sf["A+B"], m["A"] * m["B"]).pvalue < 0.05   # wells of one plate as n
    across = stats.ttest_1samp(lb, 0).pvalue < 0.05
    hsa_syn = stats.ttest_1samp(lh, 0).pvalue < 0.05 and np.mean(lh) < 0
    return plate1, across, hsa_syn
for lab, sf_ab in (("Independent drugs (A+B = 0.30, the Bliss expectation)", 0.30), ("True synergy (A+B = 0.225, 25% below the Bliss expectation)", 0.225)):
    r = np.array([study(sf_ab) for _ in range(R)])
    print(f"{lab}: p<0.05 for a Bliss difference with one plate's wells as n {r[:, 0].mean():.1%}; "
          f"with one ratio per experiment (n = 5) {r[:, 1].mean():.1%}; HSA 'synergy' across experiments {r[:, 2].mean():.1%}")