# 10,000 simulated studies with the counts.csv design (16 WT and 16 KO mice, 4 cages per genotype, 30 taxa).
# 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 = 10_000

def one(bloom):
    rows = study(rng, bloom=bloom)
    x = np.array([r[3] for r in rows], float); ko = np.array([r[1] == "KO" for r in rows]); cage = np.array([r[2] for r in rows])
    rel = x / x.sum(1, keepdims=True); lx = np.log(x + 0.5); clr = lx - lx.mean(1, keepdims=True)
    p_rel = stats.ttest_ind(rel[ko], rel[~ko], equal_var=False).pvalue
    q_mouse = stats.false_discovery_control(stats.ttest_ind(clr[ko], clr[~ko], equal_var=False).pvalue)
    uc = list(dict.fromkeys(cage)); cm = np.array([clr[cage == c].mean(0) for c in uc]); kc = np.array([c.startswith("KO") for c in uc])
    q_cage = stats.false_discovery_control(stats.ttest_ind(cm[kc], cm[~kc], equal_var=False).pvalue)
    other = slice(1, None)                                   # taxa 2-30 are unchanged in absolute terms
    return [(p_rel[other] < 0.05).any(), (q_mouse[other] < 0.05).any(), (q_cage[other] < 0.05).any(), q_cage[0] < 0.05]
for lab, bloom in (("No genotype effect", 1.0), ("KO raises taxon 1 four-fold", 4.0)):
    r = np.array([one(bloom) for _ in range(R)])
    print(f"{lab}: at least one unchanged taxon called different: relative abundance, mice as n, uncorrected {r[:, 0].mean():.1%}; "
          f"CLR, mice as n, BH {r[:, 1].mean():.1%}; CLR, cages as n, BH {r[:, 2].mean():.1%}; taxon 1 detected (CLR, cages as n, BH) {r[:, 3].mean():.1%}")