# 10,000 simulated studies like make_data.py but with NO true difference (both compounds IC50 100 nM).
# Share of studies declared "different" (p < 0.05) by the pooled-wells F test and by the paired t-test on pIC50.
import numpy as np, warnings
from scipy import stats
from scipy.optimize import curve_fit
warnings.filterwarnings("ignore")
rng = np.random.default_rng(5)
lx = np.repeat(np.arange(-9, -4.9, 0.5), 3)
def f4(x, b, t, i, h): return b + (t - b) / (1 + 10 ** ((x - i) * h))
def full(X, bA, tA, iA, hA, bB, tB, iB, hB):
    x, g = X; return np.where(g, f4(x, bB, tB, iB, hB), f4(x, bA, tA, iA, hA))
def shared(X, bA, tA, hA, bB, tB, hB, i):
    x, g = X; return np.where(g, f4(x, bB, tB, i, hB), f4(x, bA, tA, i, hA))
runs, hit_f, hit_t = 10_000, 0, 0
for _ in range(runs):
    X, Y, G, pa, pb = [], [], [], [], []
    for e in range(4):
        shift = rng.normal(0, 0.15)
        for g in (0, 1):
            li = -7 + shift + rng.normal(0, 0.08)
            y = f4(lx, 2, 100, li, 1) + rng.normal(0, 5, lx.size)
            X.append(lx); Y.append(y); G.append(np.full(lx.size, g))
            ym = y.reshape(-1, 3).mean(1)
            p, _ = curve_fit(f4, lx[::3], ym, p0=[0, 100, -7, 1], maxfev=10000)
            (pa if g == 0 else pb).append(-p[2])
    x, y, g = np.concatenate(X), np.concatenate(Y), np.concatenate(G).astype(bool)
    pf, _ = curve_fit(full, (x, g), y, p0=[0, 100, -7, 1, 0, 100, -7, 1], maxfev=20000)
    ps, _ = curve_fit(shared, (x, g), y, p0=[0, 100, 1, 0, 100, 1, -7], maxfev=20000)
    rf, rs = np.sum((y - full((x, g), *pf)) ** 2), np.sum((y - shared((x, g), *ps)) ** 2)
    df2 = y.size - 8
    hit_f += stats.f.sf((rs - rf) / (rf / df2), 1, df2) < 0.05
    hit_t += stats.ttest_1samp(np.array(pa) - np.array(pb), 0).pvalue < 0.05
print(f"runs = {runs}; declared different with no true difference: pooled-wells F test {hit_f / runs:.1%}, paired t-test on pIC50 {hit_t / runs:.1%}")