# How often each analysis gives p < 0.05: 5,000 simulated studies per row of 3, 4 or 6 independent experiments (one
# plate each: 3 blank, 6 vehicle and 6 treated wells), with no effect and with treated viability 85% of vehicle.
# Python 3.13, NumPy, SciPy. Seed 20261112.
import numpy as np
from scipy import stats
from model import study
rng = np.random.default_rng(20261112); R = 5_000

def one(n_exp, ratio):
    rows = study(rng, n_exp, ratio)
    e = np.array([r[0] for r in rows]); c = np.array([r[1] for r in rows]); a = np.array([r[3] for r in rows])
    E = np.unique(e); bl = {k: a[(e == k) & (c == "blank")].mean() for k in E}; s = a - np.array([bl[k] for k in e])
    veh = np.array([s[(e == k) & (c == "vehicle")].mean() for k in E]); trt = np.array([s[(e == k) & (c == "treated")].mean() for k in E])
    vm = dict(zip(E, veh)); pct = 100 * s / np.array([vm[k] for k in e])
    pa = stats.ttest_1samp(pct[c == "treated"], 100).pvalue
    pb = stats.ttest_ind(pct[c == "treated"], pct[c == "vehicle"], equal_var=False).pvalue
    pc = stats.ttest_1samp(100 * trt / veh, 100).pvalue
    pd_ = stats.ttest_1samp(np.log(trt / veh), 0).pvalue
    pe = stats.ttest_rel(trt, veh).pvalue
    return np.array([pa, pb, pc, pd_, pe]) < 0.05

for n_exp in (3, 4, 6):
    for ratio, lab in ((1.0, "No effect (false positives)"), (0.85, "Viability 85% of vehicle (power)")):
        r = np.array([one(n_exp, ratio) for _ in range(R)]).mean(axis=0)
        print(f"{lab}, {n_exp} experiments: A wells as n, vs 100% {r[0]:.1%}; B wells as n, vs vehicle wells {r[1]:.1%}; "
              f"C experiments, percent vs 100% {r[2]:.1%}; D experiments, log ratio {r[3]:.1%}; E experiments, paired raw signal {r[4]:.1%}")