# 10,000 simulated studies with the seahorse.csv design (5 plates, 10 WT and 10 KO wells per plate; KO wells hold 25% fewer
# 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 = 10_000

def one(ko_rate):
    rows = np.array([(r[0], r[1] == "KO", r[3], r[6] - r[7]) for r in study(rng, ko_rate=ko_rate)], float)
    plate, ko, cells, mx = rows.T; ko = ko.astype(bool); raw = np.log(mx); norm = np.log(mx / cells * 1000)
    out = [stats.ttest_ind(raw[ko], raw[~ko], equal_var=False).pvalue < 0.05, stats.ttest_ind(norm[ko], norm[~ko], equal_var=False).pvalue < 0.05]
    for v in (raw, norm):
        r = [v[(plate == p) & ko].mean() - v[(plate == p) & ~ko].mean() for p in np.unique(plate)]
        out.append(stats.ttest_1samp(r, 0).pvalue < 0.05)
    return out
for lab, rate in (("Same respiration per cell", 1.0), ("KO maximal respiration per cell 20% lower", 0.8)):
    r = np.array([one(rate) for _ in range(R)])
    print(f"{lab}: p<0.05 raw OCR wells as n {r[:, 0].mean():.1%}; per cell, wells as n {r[:, 1].mean():.1%}; "
          f"raw OCR plates as n {r[:, 2].mean():.1%}; per cell, plates as n {r[:, 3].mean():.1%}")