# 2,000 simulated studies with the log2_expression.csv design (4 vs 4 samples, 2,000 genes, the same true fold changes).
# Python 3.13, NumPy, SciPy. Seed fixed before running.
import numpy as np
from scipy import stats
from model import study, true_lfc
rng = np.random.default_rng(20261112); R = 2_000; truth = true_lfc(); tau = 1.0

def one(n):
    ctl, trt = study(rng) if n == 4 else study_n(n)
    df = 2 * n - 2; lfc = trt.mean(0) - ctl.mean(0)
    se = np.sqrt((ctl.var(0, ddof=1) + trt.var(0, ddof=1)) / 2 * (2 / n))
    q = stats.false_discovery_control(2 * stats.t.sf(np.abs(lfc) / se, df))
    qt = stats.false_discovery_control(stats.t.sf((np.abs(lfc) - tau) / se, df) + stats.t.sf((np.abs(lfc) + tau) / se, df))
    out = []
    for called in (q < 0.05, (q < 0.05) & (np.abs(lfc) > tau), qt < 0.05):
        k = called.sum(); out += [k, np.sum(called & (np.abs(truth) < tau)) / max(k, 1), np.sum(called & (np.abs(truth) > tau)) / np.sum(np.abs(truth) > tau)]
    return out
def study_n(n):
    import model; old = model.N; model.N = n
    try: return model.study(rng)
    finally: model.N = old
for n in (4, 8):
    r = np.array([one(n) for _ in range(R)])
    labs = ("padj < 0.05", "padj < 0.05 and |LFC| > 1 after testing", "TREAT |LFC| > 1, adjusted p < 0.05")
    print(f"{n} vs {n} samples: " + "; ".join(f"{l}: mean {r[:, 3 * i].mean():.0f} genes, share with true change below 2-fold {r[:, 3 * i + 1].mean():.1%}, "
          f"genes truly above 2-fold found {r[:, 3 * i + 2].mean():.1%}" for i, l in enumerate(labs)))