"""Pilot power analysis (simulation) for the scenario-weighted paired contrast of acceptance probabilities.
Model: scenario random effect on the logit scale (sd = tau), arm effect delta (logit), R repeats, J scenarios; test = scenario-cluster
bootstrap percentile CI (approximated here by the paired t-test on scenario-level differences, which is conservative-equivalent at J = 10).
Usage: python power_simulation.py --out ../../analysis/tables/power_pilot.csv"""
from __future__ import annotations
import argparse, itertools, numpy as np, pandas as pd
from scipy import stats

def simulate(J, R, p_base, delta_pp, tau, n_sim=4000, seed=7, alpha=0.05):
    rng = np.random.default_rng(seed); rej = 0
    logit = lambda p: np.log(p / (1 - p)); expit = lambda x: 1 / (1 + np.exp(-x))
    b0 = logit(p_base); d = logit(min(0.999, p_base + delta_pp / 100)) - b0
    for _ in range(n_sim):
        u = rng.normal(0, tau, J)
        p0 = expit(b0 + u); p1 = expit(b0 + d + u)
        y0 = rng.binomial(R, p0) / R; y1 = rng.binomial(R, p1) / R
        diff = y1 - y0
        if diff.std(ddof=1) == 0: rej += int(diff.mean() > 0); continue
        t = diff.mean() / (diff.std(ddof=1) / np.sqrt(J)); pval = 2 * (1 - stats.t.cdf(abs(t), J - 1))
        rej += int(pval < alpha and diff.mean() > 0)
    return rej / n_sim

if __name__ == "__main__":
    ap = argparse.ArgumentParser(); ap.add_argument("--out", required=True); a = ap.parse_args()
    rows = []
    for J, R, p, dpp, tau in itertools.product([10, 20, 40], [2, 3, 5], [0.3, 0.5], [15, 25, 35], [0.8]):
        rows.append({"J_scenarios": J, "R_repeats": R, "p_base": p, "delta_pp": dpp, "tau_logit": tau, "power": simulate(J, R, p, dpp, tau)})
    df = pd.DataFrame(rows); df.to_csv(a.out, index=False); print(df.to_string())