# How often each analysis gives p < 0.05 (D: q < 0.05 at any radius): 5,000 simulated studies of 6 + 6 mice with
# 8 neurons each, with no knockout effect and with 15% fewer intersections at every radius; then the 15% effect with
# twice the neurons per mouse, or twice the mice. Python 3.13, NumPy, SciPy.
# Seed 20261112.
import numpy as np
from scipy import stats
from model import RADII, NEURONS, MICE, profile
rng = np.random.default_rng(20261112); R = 5_000

def counts(ratio, mice, neurons, nm=8.0, rc=60.0, sd_mouse=0.15, sd_neuron=0.25, sd_rc=0.10):   # same draws as model.study, as arrays
    nm_m = nm * ratio * np.exp(rng.normal(0, sd_mouse, (mice, 1, 1)))
    nm_n = nm_m * np.exp(rng.normal(0, sd_neuron, (mice, neurons, 1))); rc_n = rc * np.exp(rng.normal(0, sd_rc, (mice, neurons, 1)))
    x = RADII[None, None, :] / rc_n; return rng.poisson(nm_n * x ** 2.0 * np.exp(2.0 * (1 - x)))

def one(ko_ratio, mice=MICE, neurons=NEURONS):
    wt, ko = counts(1.0, mice, neurons), counts(ko_ratio, mice, neurons)                               # [mouse, neuron, radius]
    cw, ck = wt.reshape(-1, len(RADII)), ko.reshape(-1, len(RADII))
    pa = stats.ttest_ind(ck, cw, axis=0, equal_var=False).pvalue.min()
    pb = stats.ttest_ind(ck.sum(1), cw.sum(1), equal_var=False).pvalue
    mw, mk = wt.mean(1), ko.mean(1)
    pc = stats.ttest_ind(mk.sum(1), mw.sum(1), equal_var=False).pvalue
    qd = stats.false_discovery_control(stats.ttest_ind(mk, mw, axis=0, equal_var=False).pvalue).min()
    return np.array([pa, pb, pc, qd]) < 0.05
for lab, kr, mice, neurons in (("Knockout has no effect, 6 mice x 8 neurons", 1.0, 6, 8), ("15% fewer intersections, 6 mice x 8 neurons", 0.85, 6, 8),
                               ("15% fewer intersections, 6 mice x 16 neurons", 0.85, 6, 16), ("15% fewer intersections, 12 mice x 8 neurons", 0.85, 12, 8)):
    r = np.array([one(kr, mice, neurons) for _ in range(R)]).mean(axis=0)
    print(f"{lab}: A any of 20 radii, neurons as n {r[0]:.1%}; B total per neuron, neurons as n {r[1]:.1%}; "
          f"C total per mouse, mice as n {r[2]:.1%}; D radii on mouse means, Benjamini-Hochberg {r[3]:.1%}")