# How often each analysis gives p < 0.05: 5,000 simulated studies of 5 + 5 mice with 80-120 axons each, with no
# difference, with knockout axons 20% smaller but normal myelin, and with myelin 15% thinner. Python 3.13, NumPy,
# SciPy. Seed 20261112. Analysis D first used a line of g-ratio on axon diameter (13.1% false positives when only
# axon size differed, because the g-ratio is not linear in axon diameter); it now uses a line of fibre on axon diameter.
# That first version is kept as column E, which uses no random numbers, so columns A-D are unchanged.
import numpy as np
from scipy import stats
from model import study
rng = np.random.default_rng(20261112); R = 5_000; D0 = 1.0

def one(myelin_ratio, diam_ratio):
    rows = study(rng, myelin_ratio, diam_ratio)
    ko = np.array([r[0] == "knockout" for r in rows]); mouse = np.array([r[1] for r in rows]); a = np.array([r[3] for r in rows]); f = np.array([r[4] for r in rows])
    g = a / f
    pa = stats.ttest_ind(g[ko], g[~ko], equal_var=False).pvalue
    X = np.column_stack([np.ones(len(g)), ko, a]); b, *_ = np.linalg.lstsq(X, g, rcond=None); res = g - X @ b; dfr = len(g) - 3
    pb = 2 * stats.t.sf(abs(b[1] / np.sqrt(res @ res / dfr * np.linalg.inv(X.T @ X)[1, 1])), dfr)
    gm, g1, ge, mk = [], [], [], []
    for mo in np.unique(mouse):
        s = mouse == mo; c = np.polyfit(a[s], f[s], 1); gm.append(g[s].mean()); g1.append(D0 / (c[0] * D0 + c[1])); mk.append(ko[s][0])
        e = np.polyfit(a[s], g[s], 1); ge.append(e[0] * D0 + e[1])
    gm, g1, ge, mk = np.array(gm), np.array(g1), np.array(ge), np.array(mk)
    pc = stats.ttest_ind(gm[mk], gm[~mk], equal_var=False).pvalue; pd_ = stats.ttest_ind(g1[mk], g1[~mk], equal_var=False).pvalue
    pe = stats.ttest_ind(ge[mk], ge[~mk], equal_var=False).pvalue
    return np.array([pa, pb, pc, pd_, pe]) < 0.05
for lab, mr, dr in (("No difference", 1.0, 1.0), ("Knockout axons 20% smaller, myelin normal", 1.0, 0.8), ("Myelin 15% thinner", 0.85, 1.0)):
    r = np.array([one(mr, dr) for _ in range(R)]).mean(axis=0)
    print(f"{lab}: p<0.05 with A axons as n {r[0]:.1%}; B regression, axons as n {r[1]:.1%}; C mean g-ratio per mouse {r[2]:.1%}; D g-ratio at 1 um per mouse {r[3]:.1%}; "
          f"E (first version of D) line of g-ratio on axon diameter, at 1 um, per mouse {r[4]:.1%}")