# 1,000 simulated studies with the cells.csv design (250 cells per mouse, 100 genes, 5 truly changed when there is an effect).
# Fewer runs than other guides because each study has 2,000 or more cells. 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 one(effect, mice):
    rows = study(rng, effect=effect, mice=mice)
    x = np.vstack([r[2] for r in rows]).astype(float); tr = np.concatenate([[r[1] == "treated"] * r[2].shape[0] for r in rows])
    ln = np.log1p(x / x.sum(1, keepdims=True) * 1e4)
    q_cell = stats.false_discovery_control(stats.mannwhitneyu(ln[tr], ln[~tr], alternative="two-sided", axis=0).pvalue)
    pb = np.array([r[2].sum(0) for r in rows], float); tm = np.array([r[1] == "treated" for r in rows])
    cpm = np.log2(pb / pb.sum(1, keepdims=True) * 1e6 + 1)
    q_pb = stats.false_discovery_control(stats.ttest_ind(cpm[tm], cpm[~tm], equal_var=False).pvalue)
    return [(q_cell[5:] < 0.05).any(), (q_cell[5:] < 0.05).sum(), (q_cell[:5] < 0.05).mean(),
            (q_pb[5:] < 0.05).any(), (q_pb[:5] < 0.05).mean()]
for lab, effect, mice in (("No true differences, 4 mice per group", 0.0, 4), ("5 genes doubled, 4 mice per group", 1.0, 4),
                          ("5 genes doubled, 8 mice per group", 1.0, 8)):
    r = np.array([one(effect, mice) for _ in range(R)], float)
    print(f"{lab}: cells as replicates: at least one unchanged gene called {r[:, 0].mean():.1%} (mean {r[:, 1].mean():.1f} unchanged genes called), "
          f"changed genes detected {r[:, 2].mean():.1%}; pseudobulk: at least one unchanged gene called {r[:, 3].mean():.1%}, changed genes detected {r[:, 4].mean():.1%}")