# How often each analysis gives p < 0.05: 5,000 simulated studies of 4 transfections, with no effect, with the treatment
# doubling firefly, and with the treatment raising only the Renilla control by 50%; then the doubling with
# 3 transfections. Python 3.13, NumPy, SciPy.
# Seed 20261112.
import numpy as np
from scipy import stats
from model import day, WELLS
rng = np.random.default_rng(20261112); R = 5_000; T = 4

def one(ff, rl, T=T):
    x = np.array([[r[2:] for r in day(rng, ff, rl)] for _ in range(T)], float).reshape(T, 2, WELLS, 2)   # [transfection, cond, well, reporter]
    ratio = x[..., 0] / x[..., 1]; m_ff, m_rl, m_rt = x[..., 0].mean(2), x[..., 1].mean(2), ratio.mean(2)
    p1 = lambda v: stats.ttest_1samp(np.log2(v[:, 1] / v[:, 0]), 0).pvalue
    pb = stats.ttest_ind(ratio[:, 1].ravel(), ratio[:, 0].ravel(), equal_var=False).pvalue
    lc = np.log2(m_rt[:, 1] / m_rt[:, 0])
    return np.array([p1(m_ff) < 0.05, pb < 0.05, p1(m_rt) < 0.05, p1(m_rl) < 0.05, (p1(m_rt) < 0.05) and lc.mean() < 0])
for lab, ff, rl, nt in (("No effect", 1.0, 1.0, 4), ("Treatment doubles firefly", 2.0, 1.0, 4), ("Treatment raises Renilla by 50% only", 1.0, 1.5, 4),
                        ("Treatment doubles firefly, 3 transfections", 2.0, 1.0, 3)):
    r = np.array([one(ff, rl, nt) for _ in range(R)]).mean(axis=0)
    print(f"{lab}: A firefly alone {r[0]:.1%}; B ratio, wells as n {r[1]:.1%}; C ratio, paired log2 {r[2]:.1%} "
          f"(of which 'repression' {r[4]:.1%}); D Renilla check {r[3]:.1%}")