# 2,000 simulated studies with the binding.csv design (5 experiments, control and treated membranes, 8 concentrations,
# triplicate wells). Python 3.13, NumPy, SciPy. Seed fixed before running.
import numpy as np
from scipy import stats
from scipy.optimize import curve_fit
from model import experiment, L, WELLS
rng = np.random.default_rng(20261112); R = 2_000
one_site = lambda l, bmax, kd: bmax * l / (kd + l)

def specific(rows):
    tot = np.array([r[2] for r in rows]).reshape(len(L), WELLS); nsb = np.array([r[3] for r in rows]).reshape(len(L), WELLS)
    return np.repeat(L, WELLS), (tot - nsb.mean(1, keepdims=True)).ravel()
def one(bmax_ratio, kd_ratio):
    per, pooled = [], []
    for e in range(5):
        rows = experiment(rng, bmax_ratio=bmax_ratio, kd_ratio=kd_ratio)
        fits = {}
        for c in ("control", "treated"):
            l, s = specific([r[1:] for r in rows if r[0] == c]); p, _ = curve_fit(one_site, l, s, p0=[s.max(), 1.0], maxfev=10_000)
            fits[c] = p; pooled.append((l, s, len(pooled), int(c == "treated")))
        per.append([np.log(fits["treated"][0] / fits["control"][0]), np.log(fits["treated"][1] / fits["control"][1])])
    per = np.array(per)
    l = np.concatenate([x[0] for x in pooled]); s = np.concatenate([x[1] for x in pooled])
    g = np.concatenate([np.full(len(x[0]), x[2]) for x in pooled]); t = np.concatenate([np.full(len(x[0]), x[3]) for x in pooled]); ng = len(pooled)
    full = lambda X, *q: np.array(q)[g] * X / (np.array(q)[ng + t] + X); red = lambda X, *q: np.array(q)[g] * X / (q[ng] + X)
    pf, _ = curve_fit(full, l, s, p0=[500] * ng + [1, 1], maxfev=20_000); pr, _ = curve_fit(red, l, s, p0=[500] * ng + [1], maxfev=20_000)
    ssf, ssr = np.sum((s - full(l, *pf)) ** 2), np.sum((s - red(l, *pr)) ** 2); dfe = len(s) - ng - 2
    return (stats.ttest_1samp(per[:, 0], 0).pvalue < 0.05, stats.ttest_1samp(per[:, 1], 0).pvalue < 0.05,
            stats.f.sf((ssr - ssf) / (ssf / dfe), 1, dfe) < 0.05)
for lab, br, kr in (("Bmax 30% lower, Kd unchanged", 0.7, 1.0), ("Kd 1.5 times higher, Bmax unchanged", 1.0, 1.5)):
    r = np.array([one(br, kr) for _ in range(R)])
    print(f"{lab}: p<0.05 for Bmax, one ratio per experiment {r[:, 0].mean():.1%}; for Kd, one ratio per experiment {r[:, 1].mean():.1%}; "
          f"for Kd, global fit pooling wells {r[:, 2].mean():.1%}")