# 10,000 simulated studies with the ipsc.csv design (3 patient and 3 control donors, 2 clones each, 3 differentiations
# per clone), plus an isogenic design (3 patient donors, each with a gene-corrected line). Python 3.13, NumPy, SciPy.
# Seed fixed before running.
import numpy as np
from scipy import stats
from model import case_control, isogenic
rng = np.random.default_rng(20261112); R = 10_000

def cc(effect, donors):
    rows = case_control(rng, effect=effect, donors=donors)
    st = np.array([r[0] == "patient" for r in rows]); v = np.array([r[4] for r in rows])
    donor = np.array([r[1] for r in rows]); clone = np.array([r[2] for r in rows])
    out = [stats.ttest_ind(v[st], v[~st], equal_var=False).pvalue < 0.05]
    for key in (clone, donor):
        u = list(dict.fromkeys(key)); m = np.array([v[key == k].mean() for k in u]); s = np.array([st[key == k][0] for k in u])
        out.append(stats.ttest_ind(m[s], m[~s], equal_var=False).pvalue < 0.05)
    return out
def iso(effect, donors):
    rows = isogenic(rng, effect=effect, donors=donors)
    line = np.array([r[0] for r in rows]); donor = np.array([r[1] for r in rows]); v = np.array([r[3] for r in rows])
    u = list(dict.fromkeys(donor))
    diff = [v[(donor == k) & (line == "patient")].mean() - v[(donor == k) & (line == "corrected")].mean() for k in u]
    return stats.ttest_1samp(diff, 0).pvalue < 0.05
for lab, eff in (("No disease effect", 0.0), ("Patient lines 30% lower (log effect -0.36)", np.log(0.7))):
    for donors in (3, 6):
        a = np.array([cc(eff, donors) for _ in range(R)]); b = np.mean([iso(eff, donors) for _ in range(R)])
        print(f"{lab}, {donors} donors per group: p<0.05 with differentiations as n {a[:, 0].mean():.1%}; clones as n {a[:, 1].mean():.1%}; "
              f"donors as n {a[:, 2].mean():.1%}; isogenic pairs ({donors} patient donors, paired) {b:.1%}")