# How often each analysis gives p < 0.05: 5,000 simulated studies of 4 experiments, with no drug effect and with the
# drug halving binding, and the halving with 6 and with 3 experiments. Python 3.13, NumPy, SciPy. Seed 20261112.
import numpy as np
from scipy import stats
from model import experiment, REPS
rng = np.random.default_rng(20261112); R = 5_000; DF = np.log2(100)

def one(ratio, E):
    ct = np.array([[r[3] for r in experiment(rng, drug_ratio=ratio)] for _ in range(E)]).reshape(E, 2, 3, REPS)   # [exp, cond, sample, well]
    m = ct.mean(3); cin, cip, cig = m[:, :, 0], m[:, :, 1], m[:, :, 2]
    pct = 100 * 2 ** ((cin - DF) - cip); fold = 2 ** (cig - cip)
    wells = 100 * 2 ** ((cin[:, :, None] - DF) - ct[:, :, 1, :])
    pa = stats.ttest_ind(wells[:, 1].ravel(), wells[:, 0].ravel(), equal_var=False).pvalue
    pb = stats.ttest_ind(np.log2(pct[:, 1]), np.log2(pct[:, 0]), equal_var=False).pvalue
    pc = stats.ttest_1samp(np.log2(pct[:, 1] / pct[:, 0]), 0).pvalue
    pd_ = stats.ttest_1samp(np.log2(fold[:, 1] / fold[:, 0]), 0).pvalue
    return np.array([pa, pb, pc, pd_]) < 0.05
for lab, ratio, E in (("No drug effect, 4 experiments", 1.0, 4), ("Drug halves binding, 4 experiments", 0.5, 4), ("Drug halves binding, 6 experiments", 0.5, 6),
                   ("Drug halves binding, 3 experiments", 0.5, 3)):
    r = np.array([one(ratio, E) for _ in range(R)]).mean(axis=0)
    print(f"{lab}: A wells as n {r[0]:.1%}; B experiments unpaired (log2) {r[1]:.1%}; C log2 ratio of % input, paired {r[2]:.1%}; D log2 ratio of fold over IgG, paired {r[3]:.1%}")