# (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).
import json, numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

plt.rcParams.update({"font.size": 11, "figure.dpi": 150, "axes.grid": True,
                     "grid.alpha": 0.3})
R = json.load(open("results.json"))
C = json.load(open("cutoff_clean.json"))
ns = [10,20,30,50,100,200]

# ---------- Figure 1: Shapiro-Wilk pass rate vs n, per CV ----------
fig, ax = plt.subplots(figsize=(7,4.5))
sw = R["shapiro_pass_rate"]
cmap = plt.cm.viridis(np.linspace(0,0.9,len(sw)))
for col,(cv,d) in zip(cmap, sw.items()):
    y=[d[str(n)] for n in ns]
    ax.plot(ns, y, "-o", color=col, label=f"CV={float(cv):.2g}", ms=4)
ax.axhline(0.5, ls="--", c="grey", lw=1)
ax.set_xlabel("sample size per group, n")
ax.set_ylabel("P(Shapiro–Wilk passes, p>0.05)")
ax.set_title("Non-normality hides at small n\n(lognormal data 'passing' as normal)")
ax.set_xscale("log"); ax.set_xticks(ns); ax.set_xticklabels(ns)
ax.legend(title="CV = σ/μ", fontsize=8, ncol=2)
fig.tight_layout(); fig.savefig("fig1_shapiro.png"); plt.close()

# ---------- Figure 2: power raw vs log at true log-shift ln(1.5) ----------
fig, ax = plt.subplots(figsize=(7,4.5))
P = R["power_logshift_ln1.5"]
cvs = sorted(P.keys(), key=float)
x = np.arange(len(cvs)); w=0.35
n_show = "50"
raw=[P[cv][n_show]["raw_welch"] for cv in cvs]
log=[P[cv][n_show]["log"] for cv in cvs]
ax.bar(x-w/2, raw, w, label="raw-scale t-test", color="#c44e52")
ax.bar(x+w/2, log, w, label="log-scale t-test", color="#4c72b0")
ax.set_xticks(x); ax.set_xticklabels([f"{float(c):.2g}" for c in cvs])
ax.set_xlabel("CV = σ/μ"); ax.set_ylabel(f"power (n={n_show}/group)")
ax.set_title("Analyzing in the natural domain forfeits power\n(true median ratio 1.5)")
ax.legend()
fig.tight_layout(); fig.savefig("fig2_power.png"); plt.close()

# ---------- Figure 3: the cutoff (2 panels) ----------
fig, (a1,a2) = plt.subplots(1,2, figsize=(11,4.5))
cvg = sorted([float(c) for c in next(iter(C.values())).keys()])
cmap = plt.cm.plasma(np.linspace(0,0.85,len(C)))
for col,(n,d) in zip(cmap, C.items()):
    a1.plot(cvg, [d[str(c) if str(c) in d else f"{c}"]["sw_flag"] if False else d[str(c)]["sw_flag"] for c in cvg],
            "-o", color=col, ms=4, label=f"n={n}")
    a2.plot(cvg, [d[str(c)]["log_gain"] for c in cvg], "-o", color=col, ms=4, label=f"n={n}")
for a in (a1,a2):
    a.axvspan(0.5,1.0, color="gold", alpha=0.15)
    a.set_xlabel("CV = σ/μ (each group)")
a1.axhline(0.8, ls="--", c="grey", lw=1)
a1.set_ylabel("P(Shapiro–Wilk flags ≥1 group)"); a1.set_title("Detectability"); a1.legend(fontsize=9)
a2.axhline(0.05, ls="--", c="grey", lw=1)
a2.set_ylabel("power gained by log-transform"); a2.set_title("Cost of staying in natural domain"); a2.legend(fontsize=9)
fig.suptitle("The cutoff region: CV ≈ 0.5–1.0 (shaded)", y=1.02)
fig.tight_layout(); fig.savefig("fig3_cutoff.png", bbox_inches="tight"); plt.close()

# ---------- Figure 4: type I error vs CV ----------
fig, ax = plt.subplots(figsize=(7,4.5))
tu=R["typeI_unequal_n_30_10"]; te=R["typeI_equal_n_20_20"]
cvs=sorted(tu.keys(), key=float); x=[float(c) for c in cvs]
ax.plot(x,[tu[c]["raw_welch"] for c in cvs],"-o",c="#c44e52",label="raw Welch, n=30 vs 10")
ax.plot(x,[te[c]["raw_student"] for c in cvs],"-s",c="#dd8452",label="raw Student, n=20 vs 20")
ax.plot(x,[te[c]["log"] for c in cvs],"-^",c="#4c72b0",label="log-scale, n=20 vs 20")
ax.axhline(0.05, ls="--", c="grey", lw=1)
ax.set_xlabel("CV = σ/μ"); ax.set_ylabel("type I error (nominal 0.05)")
ax.set_title("Under the null, raw-scale t-tests miscalibrate as CV grows")
ax.legend(fontsize=9)
fig.tight_layout(); fig.savefig("fig4_typeI.png"); plt.close()
print("figures done")
