"""False-positive rates for common analyses of normalised (fold-change) data.
Page: /guides/statistics-on-fold-change/. Python 3.13, NumPy, SciPy 1.18.1. Seed 20261110.

Design: k independent experiments (e.g. blots), each with one control and one treated lane.
Values are lognormal: log(value) = blot effect + lane noise (+ treatment effect, here 0).
With no true treatment effect, every rejection at p < 0.05 is a false positive.
Run from this directory: python3 simulation.py
"""
import json

import numpy as np
from scipy import stats

RUNS = 20_000
rng = np.random.default_rng(20261110)


def null_pvalues(k, sd_blot, sd_lane):
    blot = rng.normal(0, sd_blot, k)
    ctrl = np.exp(blot + rng.normal(0, sd_lane, k))
    trt = np.exp(blot + rng.normal(0, sd_lane, k))
    ratio = trt / ctrl  # treated as a fraction of its own control
    return {
        # A: controls set to exactly 1, Student t-test of treated ratios against the 1s
        "student_vs_ones": stats.ttest_ind(ratio, np.ones(k), equal_var=True).pvalue,
        # B: one-sample t-test of raw ratios against 1
        "one_sample_raw_ratio": stats.ttest_1samp(ratio, 1.0).pvalue,
        # C: one-sample t-test of log ratios against 0 (= paired t-test on log values)
        "paired_log": stats.ttest_1samp(np.log(ratio), 0.0).pvalue,
        # D: Welch t-test on log values, ignoring which lanes share a blot
        "welch_unpaired_log": stats.ttest_ind(np.log(trt), np.log(ctrl), equal_var=False).pvalue,
    }


SCENARIOS = [
    {"k": 4, "sd_blot": 0.5, "sd_lane": 0.2},
    {"k": 4, "sd_blot": 0.5, "sd_lane": 0.5},
    {"k": 3, "sd_blot": 0.5, "sd_lane": 0.5},
    {"k": 6, "sd_blot": 0.5, "sd_lane": 0.5},
]

out = {"runs_per_scenario": RUNS, "seed": 20261110, "alpha": 0.05, "scenarios": []}
for sc in SCENARIOS:
    hits = {}
    for _ in range(RUNS):
        for name, pv in null_pvalues(**sc).items():
            hits[name] = hits.get(name, 0) + int(pv < 0.05)
    rates = {name: round(h / RUNS, 4) for name, h in hits.items()}
    out["scenarios"].append({**sc, "false_positive_rate": rates})
    print(sc, rates)


def power(k, sd_blot, sd_lane, fold):
    hits = {"paired_log": 0, "welch_unpaired_log": 0}
    for _ in range(RUNS):
        blot = rng.normal(0, sd_blot, k)
        ctrl = np.exp(blot + rng.normal(0, sd_lane, k))
        trt = np.exp(blot + np.log(fold) + rng.normal(0, sd_lane, k))
        hits["paired_log"] += int(stats.ttest_1samp(np.log(trt / ctrl), 0.0).pvalue < 0.05)
        hits["welch_unpaired_log"] += int(stats.ttest_ind(np.log(trt), np.log(ctrl), equal_var=False).pvalue < 0.05)
    return {name: round(h / RUNS, 4) for name, h in hits.items()}


out["power_1_5_fold"] = {"k": 4, "sd_blot": 0.5, "sd_lane": 0.2, "power": power(4, 0.5, 0.2, 1.5)}
print("power at 1.5-fold:", out["power_1_5_fold"])

# Skew of ratios with no effect: the arithmetic mean of ratios sits above 1, the median does not
r = np.exp(rng.normal(0, 0.5 * np.sqrt(2), 200_000))
out["ratio_skew_no_effect_sd_lane_0_5"] = {
    "mean_ratio": round(float(r.mean()), 3),
    "median_ratio": round(float(np.median(r)), 3),
}
print(out["ratio_skew_no_effect_sd_lane_0_5"])

# Larger numbers of experiments (added after review; appended so earlier results keep their random stream)
out["more_experiments"] = []
for sc in [{"k": 20, "sd_blot": 0.5, "sd_lane": 0.2}, {"k": 20, "sd_blot": 0.5, "sd_lane": 0.5},
           {"k": 40, "sd_blot": 0.5, "sd_lane": 0.2}, {"k": 40, "sd_blot": 0.5, "sd_lane": 0.5}]:
    hits = {}
    for _ in range(RUNS):
        for name, pv in null_pvalues(**sc).items():
            hits[name] = hits.get(name, 0) + int(pv < 0.05)
    rates = {name: round(h / RUNS, 4) for name, h in hits.items()}
    out["more_experiments"].append({**sc, "false_positive_rate": rates})
    print(sc, rates)

with open("simulation_output.json", "w") as fh:
    json.dump(out, fh, indent=2)