# 10,000 simulated studies with the responses.csv design (5 WT and 5 KO mice, 40 neurons each, 20 trials per neuron).
# 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(shift):
    rows = study(rng, ko_shift=shift)
    y = np.array([r[3] for r in rows]); ko = np.array([r[1] == "KO" for r in rows]); mouse = np.array([r[0] for r in rows])
    nmean = y.mean(1)
    p_neuron = stats.ttest_ind(nmean[ko], nmean[~ko], equal_var=False).pvalue
    um = list(dict.fromkeys(mouse)); mm = np.array([nmean[mouse == m].mean() for m in um]); km = np.array([m.startswith("KO") for m in um])
    p_mouse = stats.ttest_ind(mm[km], mm[~km], equal_var=False).pvalue
    odd, even = y[:, 0::2], y[:, 1::2]
    s = (stats.ttest_1samp(odd, 0, axis=1).pvalue < 0.05) & (odd.mean(1) > 0)
    return p_neuron < 0.05, p_mouse < 0.05, odd[s].mean() if s.any() else np.nan, even[s].mean() if s.any() else np.nan, nmean.mean()
for lab, shift in (("No genotype difference", 0.0), ("KO response lower by 0.04 (0.05 vs 0.01)", -0.04)):
    r = np.array([one(shift) for _ in range(R)])
    print(f"{lab}: p<0.05 with neurons as n {r[:, 0].mean():.1%}; with mice as n {r[:, 1].mean():.1%}; "
          f"neurons selected as responsive on odd trials: mean ΔF/F on the same trials {np.nanmean(r[:, 2]):.3f}, on held-out trials {np.nanmean(r[:, 3]):.3f}; "
          f"all neurons {r[:, 4].mean():.3f}")