# How often each analysis gives p < 0.05, and the average estimated ratio: 5,000 simulated studies of 4 experiments,
# with no knockdown effect and with a knockdown that lowers the cell signal by 30%. Python 3.13, NumPy, SciPy.
# Seed 20261112.
import numpy as np
from scipy import stats
from model import experiment
rng = np.random.default_rng(20261112); R = 5_000; E = 4

def one(ratio):
    cells = {"control": [], "knockdown": []}; imgs = {"control": [], "knockdown": []}; rc, rr = [], []
    for _ in range(E):
        rows = experiment(rng, ratio); per = {}
        for cond, i, c, area, iden, bgm in rows:
            lc, lr = np.log(iden - area * bgm), np.log(iden); cells[cond].append(lc); per.setdefault((cond, i), []).append((lc, lr))
        im = {k: np.mean(v, axis=0) for k, v in per.items()}
        for (cond, i), v in im.items(): imgs[cond].append(v[0])
        ex = {cond: np.mean([v for (cc, i), v in im.items() if cc == cond], axis=0) for cond in ("control", "knockdown")}
        rc.append(ex["knockdown"][0] - ex["control"][0]); rr.append(ex["knockdown"][1] - ex["control"][1])
    pa = stats.ttest_ind(cells["knockdown"], cells["control"], equal_var=False).pvalue
    pb = stats.ttest_ind(imgs["knockdown"], imgs["control"], equal_var=False).pvalue
    pc = stats.ttest_1samp(rc, 0).pvalue; pd_ = stats.ttest_1samp(rr, 0).pvalue
    return np.array([pa < 0.05, pb < 0.05, pc < 0.05, pd_ < 0.05, np.exp(np.mean(rc)), np.exp(np.mean(rr))])
for lab, ratio in (("No knockdown effect", 1.0), ("Knockdown lowers signal by 30%", 0.7)):
    r = np.array([one(ratio) for _ in range(R)]).mean(axis=0)
    print(f"{lab}: p<0.05 with A cells as n {r[0]:.1%}; B images as n {r[1]:.1%}; C experiments, CTCF {r[2]:.1%}; D experiments, no background subtraction {r[3]:.1%}; "
          f"average estimated ratio C {r[4]:.2f}, D {r[5]:.2f}")