# (c) 2026 William J. Dwyer. Free and open-source software under the MIT License (see LICENSE).
# If you use or build upon this software, please retain this notice and cite the accompanying
# methodology paper together with the archived software (Zenodo DOI: 10.5281/zenodo.20709963).
"""Clean cutoff experiment: both groups lognormal with the SAME CV = c
(no CV-ratio confound). Control mean 1.0, disease mean 1.5 (a real 50%
increase, satisfying the paper's eq-14 restriction). Vary c.

Isolates the effect of the absolute coefficient of variation (SD/mean) on:
  - Shapiro-Wilk detectability (prob >=1 group flags) by n
  - raw Welch t-test power vs log-scale t-test power to detect the difference
This answers: is there a SD/M (CV) limit beyond which the raw t-test fails?
"""
import numpy as np, json
from scipy import stats
rng = np.random.default_rng(7)

def ln_params(mu_n, cv):
    s2 = np.log(cv**2 + 1.0); return np.log(mu_n) - s2/2.0, np.sqrt(s2)
def draw(mu_n, cv, n, reps):
    me, se = ln_params(mu_n, cv); return rng.lognormal(me, se, size=(reps, n))

cv_grid = [0.2, 0.3, 0.5, 0.75, 1.0, 1.5, 2.0, 2.5, 3.0]
out = {}
for n in [15, 20, 30, 50]:
    out[n] = {}
    for c in cv_grid:
        go  = draw(1.0, c, n, 6000)
        gdz = draw(1.5, c, n, 6000)   # true 50% mean increase
        swo  = np.array([stats.shapiro(r).pvalue for r in go])  < .05
        swdz = np.array([stats.shapiro(r).pvalue for r in gdz]) < .05
        flag = float(np.mean(swo | swdz))
        pw_raw = float(np.mean(stats.ttest_ind(go, gdz, axis=1, equal_var=False).pvalue < .05))
        pw_log = float(np.mean(stats.ttest_ind(np.log(go), np.log(gdz), axis=1, equal_var=False).pvalue < .05))
        out[n][c] = dict(sw_flag=flag, power_raw=pw_raw, power_log=pw_log,
                         log_gain=round(pw_log - pw_raw, 3))
json.dump(out, open("cutoff_clean.json", "w"), indent=2)
print(' n   CV   SWflag  pow_raw  pow_log  log_gain')
for n,d in out.items():
    for c,e in d.items():
        print(f'{n:>3}  {c:>4}   {e["sw_flag"]:.2f}    {e["power_raw"]:.2f}    {e["power_log"]:.2f}    {e["log_gain"]:+.2f}')
print('done')
