# (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).
"""Figure 1 - two real-data illustrations from the Cushing's syndrome urinary
steroid data (Aitchison & Dunsmore 1975).
Layout per row: Natural-domain box (raw t-test) | Log-domain box (log t-test AND
its permutation backstop) | Mann-Whitney shown once at the bottom (scale-invariant
-- identical on either scale).
  Row 1 (pregnanetriol, a vs c): logs normal -> log-t and permutation-t agree and
        recover the effect; raw t and Mann-Whitney miss it.
  Row 2 (tetrahydrocortisone, b vs c): logs stay non-normal -> the parametric
        log-t is unreliable (and misses it), while the permutation-t and the
        scale-invariant rank test detect it.
"""
import pandas as pd, numpy as np
from scipy import stats
import matplotlib; matplotlib.use("Agg")
import matplotlib.pyplot as plt

d = pd.read_csv("Cushings.csv")
C1, C2, FLAG = "#2c6e9c", "#b45f06", "#b45f06"

def perm_log(x, y, B=50000, seed=1):
    rng = np.random.default_rng(seed); x = np.asarray(x, float); y = np.asarray(y, float)
    obs = abs(np.log(x).mean() - np.log(y).mean()); pool = np.log(np.concatenate([x, y])); n1 = len(x); c = 0
    for _ in range(B):
        p = rng.permutation(pool); c += abs(p[:n1].mean() - p[n1:].mean()) >= obs - 1e-12
    return (c + 1) / (B + 1)

def stat(x, y):
    x = np.asarray(x, float); y = np.asarray(y, float)
    return dict(p_raw=stats.ttest_ind(x, y, equal_var=False).pvalue,
                p_log=stats.ttest_ind(np.log(x), np.log(y), equal_var=False).pvalue,
                p_perm=perm_log(x, y), p_mw=stats.mannwhitneyu(x, y, alternative="two-sided").pvalue,
                cv1=x.std(ddof=1)/x.mean(), cv2=y.std(ddof=1)/y.mean(),
                swl1=stats.shapiro(np.log(x)).pvalue, swl2=stats.shapiro(np.log(y)).pvalue,
                gmr=np.exp(np.mean(np.log(y)))/np.exp(np.mean(np.log(x))))

A = d[d.Type == "a"].Pregnanetriol.values; Cc = d[d.Type == "c"].Pregnanetriol.values
B = d[d.Type == "b"].Tetrahydrocortisone.values; Cc2 = d[d.Type == "c"].Tetrahydrocortisone.values
R1 = stat(A, Cc); R2 = stat(B, Cc2)

def strip(ax, ya, yc, labels):
    rng = np.random.default_rng(60)
    for i, (yy, col) in enumerate(zip([ya, yc], [C1, C2]), start=1):
        ax.scatter(i + rng.uniform(-.10, .10, len(yy)), yy, s=42, color=col, edgecolor="white", zorder=3)
        ax.plot([i-.17, i+.17], [yy.mean(), yy.mean()], color=col, lw=2.3)
    ax.set_xlim(.5, 2.5); ax.set_xticks([1, 2]); ax.set_xticklabels(labels, fontsize=9); ax.grid(axis="y", alpha=.3)

fig = plt.figure(figsize=(10.6, 10.2))
# explicit panel rectangles [left, bottom, width, height]
axn1 = fig.add_axes([0.08, 0.605, 0.40, 0.255]); axl1 = fig.add_axes([0.57, 0.605, 0.40, 0.255])
axn2 = fig.add_axes([0.08, 0.155, 0.40, 0.255]); axl2 = fig.add_axes([0.57, 0.155, 0.40, 0.255])

# ---- Row 1: pregnanetriol ----
strip(axn1, A, Cc, ["adenoma\n(n=6)", "carcinoma\n(n=5)"]); axn1.set_ylabel("Pregnanetriol (mg/24h)"); axn1.set_title("Natural domain", fontweight="bold")
axn1.text(.03, .90, f"CV  {R1['cv1']:.2f} ⚑   |   {R1['cv2']:.2f}", transform=axn1.transAxes, va="top", fontsize=9, color=FLAG, fontweight="bold")
axn1.text(.03, .80, f"raw-scale t  p = {R1['p_raw']:.2f}\n(not significant)", transform=axn1.transAxes, va="top", fontsize=9, bbox=dict(fc="#dbe7f3", ec="#2c6e9c"))
axn1.annotate("outlier", xy=(0.965, 11.7), xytext=(1.30, 10.4), fontsize=8, color="#666", arrowprops=dict(arrowstyle="->", color="#999", mutation_scale=9, shrinkA=1, shrinkB=3))
strip(axl1, np.log(A), np.log(Cc), ["adenoma", "carcinoma"]); axl1.set_ylabel("log Pregnanetriol"); axl1.set_title("Log domain", fontweight="bold")
axl1.text(.03, .935, f"Shapiro after log  p = {R1['swl1']:.2f} | {R1['swl2']:.2f}  → no departure flagged", transform=axl1.transAxes, va="top", fontsize=8.3, color="#38761d")
axl1.text(.03, .86, f"log-scale t  p = {R1['p_log']:.3f}\npermutation t  p = {R1['p_perm']:.3f}\n(agree → trust the parametric test)", transform=axl1.transAxes, va="top", fontsize=8.7, bbox=dict(fc="#d9ead3", ec="#38761d"))

# ---- Row 2: tetrahydrocortisone ----
strip(axn2, B, Cc2, ["hyperplasia\n(n=10)", "carcinoma\n(n=5)"]); axn2.set_ylabel("Tetrahydrocortisone (mg/24h)"); axn2.set_title("Natural domain", fontweight="bold")
axn2.text(.03, .97, f"CV  {R2['cv1']:.2f}   |   {R2['cv2']:.2f} ⚑", transform=axn2.transAxes, va="top", fontsize=9, color=FLAG, fontweight="bold")
axn2.text(.03, .87, f"raw-scale t  p = {R2['p_raw']:.2f}\n(not significant)", transform=axn2.transAxes, va="top", fontsize=9, bbox=dict(fc="#dbe7f3", ec="#2c6e9c"))
axn2.annotate("outlier", xy=(2.056, 53.8), xytext=(1.40, 48), fontsize=8, color="#666", arrowprops=dict(arrowstyle="->", color="#999", mutation_scale=9, shrinkA=1, shrinkB=3))
strip(axl2, np.log(B), np.log(Cc2), ["hyperplasia", "carcinoma"]); axl2.set_ylabel("log Tetrahydrocortisone"); axl2.set_title("Log domain", fontweight="bold")
axl2.text(.03, .935, f"Shapiro after log  p = {R2['swl1']:.2f} | {R2['swl2']:.3f}  → departure flagged ✗", transform=axl2.transAxes, va="top", fontsize=8.3, color="#b45f06")
axl2.text(.03, .84, f"log-scale t  p = {R2['p_log']:.2f}  (n.s.)\npermutation t  p = {R2['p_perm']:.3f}  (sig.)\n(disagree → parametric unreliable)", transform=axl2.transAxes, va="top", fontsize=8.7, bbox=dict(fc="#fce5cd", ec="#b45f06"))

# ---- Mann-Whitney at the bottom of each row (scale-invariant) ----
mw_box = dict(boxstyle="round", fc="#fce5cd", ec="#b45f06")
fig.text(0.5, 0.528, f"Mann–Whitney U  (rank-based; identical on either scale):  p = {R1['p_mw']:.3f}   →  not significant",
         ha="center", fontsize=10, bbox=mw_box)
fig.text(0.5, 0.075, f"Mann–Whitney U  (rank-based; identical on either scale):  p = {R2['p_mw']:.3f}   →  significant",
         ha="center", fontsize=10, bbox=mw_box)

# ---- row verdicts ----
fig.text(0.5, 0.905, f"Logs are normal → the parametric log-t and its permutation backstop agree and recover the effect "
                      f"({R1['gmr']:.0f}× geometric-mean ratio); the rank test misses it",
         ha="center", fontsize=10.3, color="#274e13", fontweight="bold")
fig.text(0.5, 0.455, "Logs stay non-normal → the parametric log-t is unreliable and misses it; its permutation backstop "
                     "and the scale-invariant rank test detect it",
         ha="center", fontsize=10.3, color="#7a3f04", fontweight="bold")
fig.text(0.5, 0.965, "Cushing's syndrome urinary steroids (Aitchison & Dunsmore 1975): the CV screen flags both;\n"
                     "the post-transformation normality check decides whether the parametric log-t can be trusted",
         ha="center", fontsize=12, fontweight="bold")
fig.savefig("fig1_cushings.png", dpi=200, bbox_inches="tight", metadata={"Author": "William J. Dwyer (ORCID 0009-0004-0855-7222)", "Title": "fig1_cushings", "Subject": "Generated by fig_cushings.py for the Comparable-Strata Designer; MIT License; cite Zenodo DOI 10.5281/zenodo.20709963"})
fig.savefig("Dwyer_Figure1.tiff", dpi=600, bbox_inches="tight", pil_kwargs={"compression":"tiff_lzw"})
print("saved.  R1 raw=%.3f log=%.3f perm=%.3f mw=%.3f gmr=%.1f | R2 raw=%.3f log=%.3f perm=%.3f mw=%.3f"
      % (R1['p_raw'], R1['p_log'], R1['p_perm'], R1['p_mw'], R1['gmr'], R2['p_raw'], R2['p_log'], R2['p_perm'], R2['p_mw']))
