# How often each analysis gives p < 0.05: 5,000 simulated studies, with no knockout effect and with 20% fewer spines,
# for 6 mice x 6 neurons per genotype and for 12 mice x 3 neurons (the same number of neurons). Python 3.13, NumPy,
# SciPy. Seed 20261112.
import numpy as np
from scipy import stats
from model import study, SEGMENTS
rng = np.random.default_rng(20261112); R = 5_000

def one(ratio, mice, neurons):
    rows = study(rng, ratio, mice=mice, neurons=neurons)
    g = np.array([r[0] == "knockout" for r in rows]); mouse = np.array([r[1] for r in rows]); cell = np.array([f"{r[1]}-{r[2]}" for r in rows])
    L = np.array([r[4] for r in rows]); y = np.array([r[5] for r in rows], float)
    dens = 10 * y / L
    pa = stats.ttest_ind(dens[g], dens[~g], equal_var=False).pvalue
    def per(key):
        keys, inv = np.unique(key, return_inverse=True); sy = np.bincount(inv, y); sl = np.bincount(inv, L); gg = np.bincount(inv, g) > 0
        return 10 * sy / sl, gg
    dn, gn = per(cell); pb = stats.ttest_ind(dn[gn], dn[~gn], equal_var=False).pvalue
    dm, gm = per(mouse); pc = stats.ttest_ind(dm[gm], dm[~gm], equal_var=False).pvalue
    # Poisson model with one binary covariate and an offset has a closed form: rate ratio and likelihood-ratio statistic
    r1, r0, rall = y[g].sum() / L[g].sum(), y[~g].sum() / L[~g].sum(), y.sum() / L.sum()
    ll = lambda lam, yy, LL: np.sum(yy * np.log(lam * LL) - lam * LL)
    lr = 2 * (ll(r1, y[g], L[g]) + ll(r0, y[~g], L[~g]) - ll(rall, y, L))
    return np.array([pa, pb, pc, stats.chi2.sf(lr, 1)]) < 0.05
for lab, ratio, mice, neurons in (("No effect, 6 mice x 6 neurons", 1.0, 6, 6), ("20% fewer spines, 6 mice x 6 neurons", 0.8, 6, 6),
                                  ("No effect, 12 mice x 3 neurons", 1.0, 12, 3), ("20% fewer spines, 12 mice x 3 neurons", 0.8, 12, 3)):
    r = np.array([one(ratio, mice, neurons) for _ in range(R)]).mean(axis=0)
    print(f"{lab}: A segments as n {r[0]:.1%}; B neurons as n {r[1]:.1%}; C mice as n {r[2]:.1%}; D Poisson, all segments {r[3]:.1%}")