# 1,000 simulated studies with the intensities.csv design (1,000 proteins, 5 vs 5 samples, 100 changed 2-fold, missing
# values more likely at low intensity). Python 3.13, NumPy, SciPy. Seed fixed before running.
import numpy as np
from scipy import stats
from model import study
rng = np.random.default_rng(20261112); R = 1_000

def evaluate(p, keep, truth):
    q = np.full(len(truth), np.nan); q[keep] = stats.false_discovery_control(p); called = keep & (q < 0.05)
    return np.sum(called & (truth == 0)) / max(called.sum(), 1), np.sum(called & (truth != 0)) / np.sum(truth != 0), called.sum()
def one(lfc):
    x, truth = study(rng, lfc=lfc); x = x.T; ctl, trt = x[:, :5], x[:, 5:]
    xm = np.where(np.isnan(x), np.nanmin(x, axis=0), x)
    xs = np.where(np.isnan(x), np.nanmean(x, axis=0) - 1.8 * np.nanstd(x, axis=0, ddof=1), x)
    all_ = np.ones(len(truth), bool)
    out = []
    for xi in (xm, xs):
        with np.errstate(all="ignore"): p = stats.ttest_ind(xi[:, 5:], xi[:, :5], axis=1, equal_var=False).pvalue
        p = np.where(np.isnan(p), 1.0, p); out += evaluate(p, all_, truth)
    keep = (np.sum(~np.isnan(ctl), 1) >= 3) & (np.sum(~np.isnan(trt), 1) >= 3)
    p = np.array([stats.ttest_ind(trt[i][~np.isnan(trt[i])], ctl[i][~np.isnan(ctl[i])], equal_var=False).pvalue for i in np.where(keep)[0]])
    out += evaluate(p, keep, truth)
    return out
for lab, lfc in (("100 proteins changed 2-fold", 1.0), ("No protein changed", 0.0)):
    r = np.array([one(lfc) for _ in range(R)])
    names = ("impute sample minimum", "impute down-shifted value", "filter >= 3 valid per group, no imputation")
    print(f"{lab}: " + "; ".join(f"{n}: mean false discovery proportion {r[:, 3 * i].mean():.1%}, changed proteins found {r[:, 3 * i + 1].mean():.1%}, "
          f"studies with any call {np.mean(r[:, 3 * i + 2] > 0):.1%}" for i, n in enumerate(names)))