# How often each analysis gives p < 0.05: 5,000 simulated studies of 4 batches x 3 wells per condition, with no drug
# effect and with organoids 20% smaller; for the 20% effect, also the average estimated ratio when organoids below
# 15,000 um2 are excluded (a size threshold); then the 20% effect with 8 batches. Python 3.13, NumPy, SciPy. Seed 20261112.
import numpy as np
from scipy import stats
from model import batch
rng = np.random.default_rng(20261112); R = 5_000; B = 4; THRESH = 15000

def summaries(rows, thresh=0):
    keep = [(c, w, np.log(a)) for c, w, a in rows if a >= thresh]
    wm = {}
    for c, w, la in keep: wm.setdefault((c, w), []).append(la)
    return keep, {k: np.mean(v) for k, v in wm.items()}

def one(ratio, B=B):
    excl = {"control": 0, "drug": 0}
    org = {"control": [], "drug": []}; wells = {"control": [], "drug": []}; lr = []; lr_t = []
    for _ in range(B):
        rows = batch(rng, ratio); keep, wm = summaries(rows); _, wmt = summaries(rows, THRESH)
        for c, w, a in rows: excl[c] += a < THRESH
        for c, w, la in keep: org[c].append(la)
        for (c, w), v in wm.items(): wells[c].append(v)
        bm = {c: np.mean([v for (cc, w), v in wm.items() if cc == c]) for c in ("control", "drug")}
        bt = {c: np.mean([v for (cc, w), v in wmt.items() if cc == c]) for c in ("control", "drug")}
        lr.append(bm["drug"] - bm["control"]); lr_t.append(bt["drug"] - bt["control"])
    pa = stats.ttest_ind(org["drug"], org["control"], equal_var=False).pvalue
    pb = stats.ttest_ind(wells["drug"], wells["control"], equal_var=False).pvalue
    pc = stats.ttest_1samp(lr, 0).pvalue
    return np.array([pa < 0.05, pb < 0.05, pc < 0.05, np.exp(np.mean(lr)), np.exp(np.mean(lr_t)), excl["control"], excl["drug"]])
for lab, ratio, nb in (("No drug effect", 1.0, 4), ("Drug makes organoids 20% smaller", 0.8, 4), ("Drug makes organoids 20% smaller, 8 batches", 0.8, 8)):
    r = np.array([one(ratio, nb) for _ in range(R)]).mean(axis=0)
    print(f"{lab}: p<0.05 with A organoids as n {r[0]:.1%}; B wells as n {r[1]:.1%}; C batches as n, paired {r[2]:.1%}; "
          f"average estimated ratio {r[3]:.3f} with all organoids, {r[4]:.3f} after excluding organoids below 15,000 um2 "
          f"(on average {r[5]:.1f} control and {r[6]:.1f} drug organoids excluded per study)")