# Simulated experiments with the kinetics.csv design (8 substrate concentrations, duplicate wells, 5 experiments).
# Error model (model.py): constant well-to-well SD of 3 rate units; log Km varies between days (SD 0.10) and
# preparations (SD 0.10); Vmax varies between preparations (log SD 0.15).
# 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 S, WELLS, SD_RATE, experiment
rng = np.random.default_rng(20261112); R = 10_000; NEXP = 5
mm = lambda s, vmax, km: vmax * s / (km + s)
s_rep = np.repeat(S, WELLS)

def nlr(s, v):
    p, _ = curve_fit(mm, s, v, p0=[v.max(), np.median(s)], maxfev=10_000); return p[1], np.sum((v - mm(s, *p)) ** 2)
def lb(s, v):
    slope, icpt = np.polyfit(1 / s, 1 / v, 1); return slope / icpt

# 1. Accuracy of Km from one experiment (true Km 10 uM, Vmax 100, no preparation-to-preparation variation)
est = np.array([(nlr(s_rep, v)[0], lb(s_rep, v)) for v in
                (100 * s_rep / (10 + s_rep) + rng.normal(0, SD_RATE, s_rep.size) for _ in range(R))])
for j, lab in ((0, "Nonlinear regression"), (1, "Lineweaver-Burk")):
    k = est[:, j]; k = k[np.isfinite(k)]
    print(f"{lab}: median Km {np.median(k):.2f} uM; middle 95% of estimates {np.percentile(k, 2.5):.2f} to {np.percentile(k, 97.5):.2f}; "
          f"within 20% of the true 10 uM in {np.mean(np.abs(k - 10) <= 2):.1%}; negative or above 100 uM in {np.mean((est[:, j] < 0) | (est[:, j] > 100)):.1%}")

# 2. Comparing Km between enzymes: one ratio per experiment versus a global fit sharing Km across experiments
def global_ss(s, v, grp, kgrp, nk):                  # a Vmax for each preparation, nk Km values
    ng = 2 * NEXP; f = lambda X, *q: np.array(q)[grp] * X / (np.array(q)[ng + kgrp] + X)
    q, _ = curve_fit(f, s, v, p0=[100] * ng + [10] * nk, maxfev=20_000); return np.sum((v - f(s, *q)) ** 2)
def study(km_mut):
    rows = [(e, *r) for e in range(NEXP) for r in experiment(rng, 10.0, km_mut)]
    e = np.array([r[0] for r in rows]); mut = np.array([r[1] == "mutant" for r in rows]).astype(int)
    s = np.array([r[2] for r in rows]); v = np.array([r[4] for r in rows]); grp = e * 2 + mut
    lr = [np.log(nlr(s[(e == i) & (mut == 1)], v[(e == i) & (mut == 1)])[0] / nlr(s[(e == i) & (mut == 0)], v[(e == i) & (mut == 0)])[0])
          for i in range(NEXP)]
    p_ratio = stats.ttest_1samp(lr, 0).pvalue
    ss_full, ss_red = global_ss(s, v, grp, mut, 2), global_ss(s, v, grp, 0 * mut, 1); dfe = len(v) - 2 * NEXP - 2
    p_global = stats.f.sf((ss_red - ss_full) / (ss_full / dfe), 1, dfe)
    return p_ratio < 0.05, p_global < 0.05
for lab, km in (("Same Km in both enzymes (10 uM)", 10.0), ("Mutant Km 1.5 times higher (15 uM)", 15.0)):
    r = np.array([study(km) for _ in range(R)])
    print(f"{lab}: p<0.05 with one Km ratio per experiment (paired t-test, n = 5) {r[:, 0].mean():.1%}; "
          f"with a global fit sharing Km across experiments (extra-sum-of-squares F-test) {r[:, 1].mean():.1%}")