#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Class differences in modern bone/tissue collagen atomic C:N and their
amino-acid basis.

PRE-REGISTRATION (frozen before any test was run)
-------------------------------------------------
Source data: a published compilation of modern vertebrate collagen
amino-acid compositions and elemental data (Table S1, n=436 samples,
194 unique taxa), distributed under CC0 with the article
DOI 10.1111/2041-210X.13433. Amino-acid values are residues per 1000.
The atomic C:N ratio ("C/N" column) is the standard quality-control (QC)
gatekeeper used to accept or reject collagen for stable-isotope work.

Question (chosen from the real columns; not a restatement of the source
abstract): Within this modern-collagen compilation,
  (Q1) does the atomic C:N ratio differ between ray-finned fish
       (Actinopterygii) and mammals (Mammalia)? and
  (Q2) is variation in C:N associated with amino-acid composition, and in
       particular with glycine -- the most abundant residue in collagen and
       the one with the lowest carbon-to-nitrogen atom ratio (2 C : 1 N),
       so that glycine-richer collagen is expected to carry a *lower* C:N?

Frozen hypotheses, tests, expected outputs:
  T1  Fish vs mammal C:N. Mann-Whitney U, two-sided. Effect size: Cliff's
      delta and Hodges-Lehmann median shift with 95% bootstrap CI.
      Expected direction (from Q2 mechanism + higher fish glycine): fish C:N
      < mammal C:N. Output: fig-1, tbl-1.
  T2  Association of C:N with glycine. Spearman rho (C:N vs Gly) overall and
      within each of the two focal classes. Expected: negative rho.
      Output: fig-2, results.json.
  T3  Amino-acid-wide screen: Spearman rho of each of 19 residues vs C:N,
      Benjamini-Hochberg FDR across the 19 tests. Expected: glycine among the
      strongest negative, FDR-significant associations. Output: fig-3, tbl-2.
  Robustness (required): re-run T1 under (R1) batch/source stratification on
      the undocumented "Material" code; (R2) collapse to one row per unique
      taxon (pseudoreplication guard); (R3) broadened groups
      fish=Actinopterygii+Elasmobranchii vs amniotes=Mammalia+Aves;
      (R4) window-free effect size only. Output: fig-4, tbl-3.

Analysis discipline: nonparametric throughout (distributions are non-normal
by Shapiro-Wilk); effect sizes + n reported alongside p; multiplicity
controlled by BH-FDR in T3; robustness battery guards batch, pseudoreplication,
unequal groups and grouping choice. All randomness seeded (SEED=0); re-running
reproduces every number in results.json.
"""

import os
import json
import hashlib
import platform
import subprocess
import urllib.request
from datetime import datetime, timezone

import numpy as np
import pandas as pd
from scipy import stats

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib.lines import Line2D

# --------------------------------------------------------------------------- #
# 0. Determinism
# --------------------------------------------------------------------------- #
SEED = 0
RNG = np.random.default_rng(SEED)
N_BOOT = 10000

HERE = os.path.dirname(os.path.abspath(__file__))
DATA_DIR = os.path.join(HERE, "data")
FIG_DIR = os.path.join(HERE, "figures")
TBL_DIR = os.path.join(HERE, "tables")
for d in (DATA_DIR, FIG_DIR, TBL_DIR):
    os.makedirs(d, exist_ok=True)

# Colourblind-safe palette (Okabe-Ito), bound to entities and reused everywhere.
C_FISH   = "#0072B2"  # blue        -> Actinopterygii / fish
C_MAMMAL = "#D55E00"  # vermillion  -> Mammalia / amniotes
C_OTHER  = "#999999"  # grey        -> other / pooled
C_SIG    = "#009E73"  # green       -> FDR-significant
C_NS     = "#BBBBBB"  # light grey  -> not significant

# --------------------------------------------------------------------------- #
# 1. Data acquisition + integrity.
#    Source: the five Table_S CSVs of Dryad DOI 10.5061/dryad.ffbg79crm,
#    obtained from the mirror on Zenodo record 4327658. The MD5 digests below
#    are the expected values pinned in this script. acquire_and_verify()
#    downloads any file not already present locally, then recomputes each
#    file's MD5 and compares it to the pinned value; the run aborts (assert)
#    if any file fails to match. The per-run comparison is recorded in
#    results.json -> data_integrity (md5_expected, md5_observed, verified).
# --------------------------------------------------------------------------- #
FILES = {
    "Table_S1.csv": ("https://zenodo.org/api/records/4327658/files/Table_S1.csv/content",
                     "b5485f7b225022480c4961981098a33c"),
    "Table_S2.csv": ("https://zenodo.org/api/records/4327658/files/Table_S2.csv/content",
                     "71742d8057362a3ffc7456e0ae88df79"),
    "Table_S3.csv": ("https://zenodo.org/api/records/4327658/files/Table_S3.csv/content",
                     "8731d29e1a97d0e83007255a815d9685"),
    "Table_S4.csv": ("https://zenodo.org/api/records/4327658/files/Table_S4.csv/content",
                     "11f8c37ddede2999ec79c607fe7c83af"),
    "Table_S5.csv": ("https://zenodo.org/api/records/4327658/files/Table_S5.csv/content",
                     "9e727589a1501b91a7955b6a0581af21"),
}


def acquire_and_verify():
    """Download each pinned file if absent; verify MD5; return status records."""
    records = []
    for name, (url, md5_expected) in FILES.items():
        path = os.path.join(DATA_DIR, name)
        if not os.path.exists(path):
            with urllib.request.urlopen(url, timeout=120) as r:
                blob = r.read()
            with open(path, "wb") as fh:
                fh.write(blob)
        with open(path, "rb") as fh:
            blob = fh.read()
        md5_got = hashlib.md5(blob).hexdigest()
        records.append({
            "file": name, "url": url, "bytes": len(blob),
            "md5_expected": md5_expected, "md5_observed": md5_got,
            "verified": bool(md5_got == md5_expected),
        })
    return records


# --------------------------------------------------------------------------- #
# 2. Effect-size helpers
# --------------------------------------------------------------------------- #
def cliffs_delta(a, b):
    """Cliff's delta = P(a>b) - P(a<b), via ranks (O(n log n))."""
    a = np.asarray(a, float); b = np.asarray(b, float)
    na, nb = len(a), len(b)
    ranks = stats.rankdata(np.concatenate([a, b]))
    u_a = ranks[:na].sum() - na * (na + 1) / 2.0
    return float((2.0 * u_a) / (na * nb) - 1.0)


def hodges_lehmann(a, b):
    """HL estimator: median of all pairwise differences a_i - b_j."""
    a = np.asarray(a, float); b = np.asarray(b, float)
    return float(np.median((a[:, None] - b[None, :]).ravel()))


def hl_boot_ci(a, b, rng, n_boot=N_BOOT, alpha=0.05):
    """Percentile bootstrap CI for the HL shift (a - b)."""
    a = np.asarray(a, float); b = np.asarray(b, float)
    est = hodges_lehmann(a, b)
    na, nb = len(a), len(b)
    boots = np.empty(n_boot)
    for i in range(n_boot):
        sa = a[rng.integers(0, na, na)]
        sb = b[rng.integers(0, nb, nb)]
        boots[i] = np.median((sa[:, None] - sb[None, :]).ravel())
    lo, hi = np.percentile(boots, [100 * alpha / 2, 100 * (1 - alpha / 2)])
    return est, float(lo), float(hi)


def mwu(a, b):
    """Mann-Whitney U two-sided; returns U (for a), p."""
    u, p = stats.mannwhitneyu(a, b, alternative="two-sided")
    return float(u), float(p)


def bh_fdr(pvals):
    """Benjamini-Hochberg q-values."""
    p = np.asarray(pvals, float)
    n = len(p)
    order = np.argsort(p)
    ranked = p[order] * n / np.arange(1, n + 1)
    q_sorted = np.minimum.accumulate(ranked[::-1])[::-1]
    q = np.empty(n)
    q[order] = np.clip(q_sorted, 0, 1)
    return q


AA = ["Asp", "Hyp", "Thr", "Ser", "Glu", "Pro", "Gly", "Ala", "Cys", "Val",
      "Met", "Ile", "Leu", "Tyr", "Phe", "Hyl", "Lys", "His", "Arg"]
AA_FULL = {
    "Asp": "Aspartate", "Hyp": "Hydroxyproline", "Thr": "Threonine",
    "Ser": "Serine", "Glu": "Glutamate", "Pro": "Proline", "Gly": "Glycine",
    "Ala": "Alanine", "Cys": "Cysteine", "Val": "Valine", "Met": "Methionine",
    "Ile": "Isoleucine", "Leu": "Leucine", "Tyr": "Tyrosine",
    "Phe": "Phenylalanine", "Hyl": "Hydroxylysine", "Lys": "Lysine",
    "His": "Histidine", "Arg": "Arginine",
}
QC_LO, QC_HI = 2.9, 3.6  # conventional archaeological C:N acceptance window


def load_main():
    """Load Table_S1, keep populated rows, clean labels.

    Taxonomic name is stripped but genuine blanks are kept as NaN (not the
    string 'nan') so they are neither counted as a spurious taxon nor collapsed
    together as one pseudo-species in the R2 robustness check.
    """
    df = pd.read_csv(os.path.join(DATA_DIR, "Table_S1.csv"), encoding="latin-1")
    df = df[df["Gly"].notna()].copy()            # 436 populated sample rows
    df["Class"] = df["Class"].astype(str).str.strip()
    tn = df["Taxonomic name"].astype("string").str.strip()
    df["Taxonomic name"] = tn.mask(tn.isin(["", "nan", "NaN"]))  # blanks -> NA
    df["C_N"] = pd.to_numeric(df["C/N"], errors="coerce")
    return df


# --------------------------------------------------------------------------- #
# 3. Plot style
# --------------------------------------------------------------------------- #
def set_style():
    plt.rcParams.update({
        "figure.dpi": 200, "savefig.dpi": 300,
        "font.size": 9, "axes.titlesize": 9, "axes.labelsize": 9,
        "xtick.labelsize": 8, "ytick.labelsize": 8, "legend.fontsize": 7.5,
        "axes.spines.top": False, "axes.spines.right": False,
        "xtick.direction": "out", "ytick.direction": "out",
        "legend.frameon": False, "figure.autolayout": False,
        "font.family": "DejaVu Sans", "savefig.bbox": "tight",
    })


# =========================================================================== #
#  MAIN
# =========================================================================== #
def main():
    set_style()
    results = {
        "meta": {
            "title": "Modern collagen atomic C:N differs between fish and mammals "
                     "and tracks glycine content",
            "generated_utc": datetime.now(timezone.utc).strftime("%Y-%m-%dT%H:%M:%SZ"),
            "seed": SEED, "n_bootstrap": N_BOOT,
            "source_doi_article": "10.1111/2041-210X.13433",
            "source_doi_data": "10.5061/dryad.ffbg79crm",
            "source_mirror": "Zenodo record 4327658",
            "amino_acid_units": "residues per 1000",
            "cn_qc_window_archaeological": [QC_LO, QC_HI],
        }
    }

    # -- 1. integrity -------------------------------------------------------- #
    integ = acquire_and_verify()
    results["data_integrity"] = integ
    assert all(r["verified"] for r in integ), "checksum verification failed"

    df = load_main()
    fish = df[df["Class"] == "Actinopterygii"]
    mam = df[df["Class"] == "Mammalia"]
    cn_fish = fish["C_N"].dropna().values
    cn_mam = mam["C_N"].dropna().values

    results["dataset"] = {
        "n_samples_total": int(len(df)),
        "n_unique_taxa": int(df["Taxonomic name"].nunique()),  # excludes NA
        "n_rows_missing_taxon": int(df["Taxonomic name"].isna().sum()),
        "n_actinopterygii": int((df["Class"] == "Actinopterygii").sum()),
        "n_mammalia": int((df["Class"] == "Mammalia").sum()),
        "class_counts": {k: int(v) for k, v in df["Class"].value_counts().items()},
    }

    # normality justification (why nonparametric)
    results["normality_shapiro"] = {
        "cn_fish": {"n": int(len(cn_fish)), "W": float(stats.shapiro(cn_fish)[0]),
                    "p": float(stats.shapiro(cn_fish)[1])},
        "cn_mammal": {"n": int(len(cn_mam)), "W": float(stats.shapiro(cn_mam)[0]),
                      "p": float(stats.shapiro(cn_mam)[1])},
    }

    # ==================================================================== #
    #  T1  fish vs mammal C:N
    # ==================================================================== #
    U, p1 = mwu(cn_fish, cn_mam)
    delta1 = cliffs_delta(cn_fish, cn_mam)
    hl1, lo1, hi1 = hl_boot_ci(cn_fish, cn_mam, np.random.default_rng(SEED))

    def grp_stats(x):
        x = np.asarray(x, float)
        return {
            "n": int(len(x)), "median": float(np.median(x)),
            "iqr_lo": float(np.percentile(x, 25)), "iqr_hi": float(np.percentile(x, 75)),
            "mean": float(np.mean(x)), "sd": float(np.std(x, ddof=1)),
            "min": float(x.min()), "max": float(x.max()),
            "pct_in_archaeo_window": float(100 * np.mean((x >= QC_LO) & (x <= QC_HI))),
        }

    results["T1_fish_vs_mammal_CN"] = {
        "test": "Mann-Whitney U (two-sided)",
        "U": U, "p": p1,
        "cliffs_delta": delta1,
        "hodges_lehmann_shift_fish_minus_mammal": hl1,
        "hl_ci95": [lo1, hi1],
        "group_fish_actinopterygii": grp_stats(cn_fish),
        "group_mammal_mammalia": grp_stats(cn_mam),
        "expected_direction": "fish < mammal",
        "observed_direction": "fish < mammal" if hl1 < 0 else "fish >= mammal",
    }

    # tbl-1
    tbl1 = pd.DataFrame([
        {"group": "Actinopterygii (fish)", **grp_stats(cn_fish)},
        {"group": "Mammalia", **grp_stats(cn_mam)},
    ])
    tbl1.to_csv(os.path.join(TBL_DIR, "tbl-1-cn-by-class.csv"), index=False)

    # fig-1
    fig, ax = plt.subplots(figsize=(5.2, 4.0))
    groups = [("Actinopterygii\n(fish)", cn_fish, C_FISH),
              ("Mammalia", cn_mam, C_MAMMAL)]
    for i, (lab, vals, col) in enumerate(groups):
        jit = (np.random.default_rng(SEED + i).random(len(vals)) - 0.5) * 0.28
        ax.scatter(np.full(len(vals), i) + jit, vals, s=10, alpha=0.35,
                   color=col, edgecolors="none", zorder=2)
        med = np.median(vals)
        q1, q3 = np.percentile(vals, [25, 75])
        ax.plot([i - 0.22, i + 0.22], [med, med], color="black", lw=2.2, zorder=4)
        ax.plot([i - 0.16, i + 0.16], [q1, q1], color="black", lw=1.0, zorder=4)
        ax.plot([i - 0.16, i + 0.16], [q3, q3], color="black", lw=1.0, zorder=4)
        ax.plot([i, i], [q1, q3], color="black", lw=1.0, zorder=3)
    ax.axhspan(QC_LO, QC_HI, color=C_OTHER, alpha=0.12, zorder=0)
    ax.axhline(QC_LO, color=C_OTHER, lw=0.8, ls="--", zorder=1)
    ax.axhline(QC_HI, color=C_OTHER, lw=0.8, ls="--", zorder=1)
    ax.text(1.46, QC_HI, "archaeological\nQC window 2.9-3.6", fontsize=6.8,
            va="top", ha="right", color="#555555")
    ax.set_xticks([0, 1]); ax.set_xticklabels([g[0] for g in groups])
    ax.set_ylabel("Atomic C:N ratio (dimensionless)")
    ax.set_xlim(-0.6, 1.6)
    ax.set_title("Modern collagen C:N is lower in fish than in mammals")
    p_txt = "p < 0.001" if p1 < 1e-3 else f"p = {p1:.3f}"
    ax.annotate(f"Mann-Whitney U, {p_txt}\nHL shift = {hl1:+.03f} "
                f"[{lo1:+.03f}, {hi1:+.03f}]\nCliff's \u03b4 = {delta1:+.2f}",
                xy=(0.02, 0.02), xycoords="axes fraction", fontsize=6.8,
                va="bottom", ha="left")
    fig.savefig(os.path.join(FIG_DIR, "fig-1-cn-by-class.png"))
    plt.close(fig)

    # ==================================================================== #
    #  T2  C:N vs glycine
    # ==================================================================== #
    def spearman(sub):
        s = sub[["C_N", "Gly"]].dropna()
        rho, p = stats.spearmanr(s["C_N"], s["Gly"])
        return {"rho": float(rho), "p": float(p), "n": int(len(s))}

    results["T2_CN_vs_glycine"] = {
        "test": "Spearman rank correlation (C:N vs glycine, residues/1000)",
        "overall": spearman(df),
        "actinopterygii": spearman(fish),
        "mammalia": spearman(mam),
        "expected_direction": "negative",
    }

    # fig-2
    fig, ax = plt.subplots(figsize=(5.4, 4.0))
    for lab, sub, col in [("Actinopterygii (fish)", fish, C_FISH),
                          ("Mammalia", mam, C_MAMMAL),
                          ("Other classes", df[~df["Class"].isin(["Actinopterygii", "Mammalia"])], C_OTHER)]:
        s = sub[["C_N", "Gly"]].dropna()
        ax.scatter(s["Gly"], s["C_N"], s=12, alpha=0.4, color=col,
                   edgecolors="none", label=f"{lab} (n={len(s)})")
    # rank-consistent robust trend lines within the two focal classes
    for sub, col in [(fish, C_FISH), (mam, C_MAMMAL)]:
        s = sub[["C_N", "Gly"]].dropna().sort_values("Gly")
        if len(s) > 5:
            res = stats.theilslopes(s["C_N"].values, s["Gly"].values)
            xs = np.array([s["Gly"].min(), s["Gly"].max()])
            ax.plot(xs, res[1] + res[0] * xs, color=col, lw=1.8, zorder=5)
    ax.set_xlabel("Glycine content (residues per 1000)")
    ax.set_ylabel("Atomic C:N ratio (dimensionless)")
    ax.set_title("Higher glycine content tracks lower C:N within each class")
    ov = results["T2_CN_vs_glycine"]["overall"]
    ax.annotate(f"overall Spearman \u03c1 = {ov['rho']:+.2f} (n={ov['n']})",
                xy=(0.02, 0.02), xycoords="axes fraction", fontsize=7, va="bottom")
    ax.legend(loc="upper right", markerscale=1.4)
    fig.savefig(os.path.join(FIG_DIR, "fig-2-cn-vs-glycine.png"))
    plt.close(fig)

    # ==================================================================== #
    #  T3  amino-acid-wide screen + BH-FDR
    # ==================================================================== #
    rows = []
    for aa in AA:
        s = df[["C_N", aa]].dropna()
        if len(s) < 10:
            rows.append({"residue": aa, "residue_full": AA_FULL[aa],
                         "rho": np.nan, "p": np.nan, "n": int(len(s))})
            continue
        rho, p = stats.spearmanr(s["C_N"], s[aa])
        rows.append({"residue": aa, "residue_full": AA_FULL[aa],
                     "rho": float(rho), "p": float(p), "n": int(len(s))})
    scr = pd.DataFrame(rows)
    valid = scr["p"].notna()
    scr.loc[valid, "q_bh"] = bh_fdr(scr.loc[valid, "p"].values)
    scr["significant_q<0.05"] = scr["q_bh"] < 0.05
    scr = scr.sort_values("rho", na_position="last").reset_index(drop=True)
    scr.to_csv(os.path.join(TBL_DIR, "tbl-2-aa-cn-spearman-fdr.csv"), index=False)

    results["T3_aa_screen"] = {
        "test": "Spearman rho per residue vs C:N; Benjamini-Hochberg FDR over "
                f"{int(valid.sum())} residues",
        "n_residues_tested": int(valid.sum()),
        "results": [
            {k: (None if (isinstance(v, float) and np.isnan(v)) else v)
             for k, v in r.items()}
            for r in scr.to_dict("records")
        ],
        "glycine": next(r for r in scr.to_dict("records") if r["residue"] == "Gly"),
    }

    # fig-3
    plot_df = scr[scr["rho"].notna()].sort_values("rho")
    fig, ax = plt.subplots(figsize=(5.6, 5.6))
    colors = [C_SIG if s else C_NS for s in plot_df["significant_q<0.05"]]
    ypos = np.arange(len(plot_df))
    ax.barh(ypos, plot_df["rho"], color=colors, edgecolor="none")
    ax.axvline(0, color="black", lw=0.8)
    ax.set_yticks(ypos)
    ax.set_yticklabels(plot_df["residue"])
    ax.set_xlabel("Spearman \u03c1 (residue vs atomic C:N)")
    ax.set_title("Amino-acid associations with C:N (Benjamini-Hochberg FDR)")
    gy = ypos[plot_df["residue"].values == "Gly"][0]
    grho = plot_df.loc[plot_df["residue"] == "Gly", "rho"].values[0]
    ax.annotate("Glycine\n(most abundant residue)", xy=(grho, gy),
                xytext=(grho + 0.06, gy + 2.1),
                fontsize=7.5, va="center", ha="left", fontweight="bold",
                arrowprops=dict(arrowstyle="-", lw=0.7))
    handles = [Line2D([0], [0], marker="s", ls="", color=C_SIG, label="q < 0.05"),
               Line2D([0], [0], marker="s", ls="", color=C_NS, label="not significant")]
    ax.legend(handles=handles, loc="lower right")
    ax.margins(y=0.01)
    fig.savefig(os.path.join(FIG_DIR, "fig-3-aa-cn-fdr.png"))
    plt.close(fig)

    # ==================================================================== #
    #  Robustness battery (re-run T1)
    # ==================================================================== #
    rob = []

    def one(label, a, b, note=""):
        if len(a) < 3 or len(b) < 3:
            return {"variant": label, "n_fish": int(len(a)), "n_mammal": int(len(b)),
                    "hl_shift": None, "ci_lo": None, "ci_hi": None,
                    "cliffs_delta": None, "U": None, "p": None, "note": note + " (insufficient n)"}
        U_, p_ = mwu(a, b)
        d_ = cliffs_delta(a, b)
        hl_, lo_, hi_ = hl_boot_ci(a, b, np.random.default_rng(SEED))
        return {"variant": label, "n_fish": int(len(a)), "n_mammal": int(len(b)),
                "hl_shift": hl_, "ci_lo": lo_, "ci_hi": hi_,
                "cliffs_delta": d_, "U": U_, "p": p_, "note": note}

    rob.append(one("Primary (all samples)", cn_fish, cn_mam))

    # R1 batch: within each Material stratum containing both classes
    for m in sorted(df["Material"].dropna().unique()):
        a = df[(df["Material"] == m) & (df["Class"] == "Actinopterygii")]["C_N"].dropna().values
        b = df[(df["Material"] == m) & (df["Class"] == "Mammalia")]["C_N"].dropna().values
        rob.append(one(f"R1 batch: Material={int(m)}", a, b, note="within-source stratum"))

    # R2 pseudoreplication: one median row per unique *named* taxon; rows with a
    # missing taxonomic name are retained as individual observations rather than
    # being merged into one pseudo-species or dropped.
    sub = df.dropna(subset=["C_N"]).copy()
    named = sub[sub["Taxonomic name"].notna()]
    unnamed = sub[sub["Taxonomic name"].isna()]
    coll_named = (named.groupby(["Class", "Taxonomic name"], as_index=False)["C_N"]
                       .median())
    coll = pd.concat([coll_named[["Class", "C_N"]],
                      unnamed[["Class", "C_N"]]], ignore_index=True)
    a = coll[coll["Class"] == "Actinopterygii"]["C_N"].values
    b = coll[coll["Class"] == "Mammalia"]["C_N"].values
    rob.append(one("R2 taxon-collapsed (median/species)", a, b,
                   note="one median row per named taxon; 6 unnamed rows kept individual"))

    # R3 broadened groups
    a = df[df["Class"].isin(["Actinopterygii", "Elasmobranchii"])]["C_N"].dropna().values
    b = df[df["Class"].isin(["Mammalia", "Aves"])]["C_N"].dropna().values
    rob.append(one("R3 fish+elasmo vs mammal+bird", a, b, note="broadened grouping"))

    results["robustness"] = {
        "R4_window_free_cliffs_delta_primary": delta1,
        "variants": rob,
    }
    tbl3 = pd.DataFrame(rob)
    tbl3.to_csv(os.path.join(TBL_DIR, "tbl-3-robustness.csv"), index=False)

    # fig-4 forest plot
    fp = [r for r in rob if r["hl_shift"] is not None]
    fig, ax = plt.subplots(figsize=(6.2, 4.2))
    ypos = np.arange(len(fp))[::-1]
    for y, r in zip(ypos, fp):
        col = C_FISH if r["variant"].startswith("Primary") else "#444444"
        ax.plot([r["ci_lo"], r["ci_hi"]], [y, y], color=col, lw=1.6, zorder=2)
        ax.scatter([r["hl_shift"]], [y], s=42, color=col, zorder=3)
    ax.axvline(0, color=C_MAMMAL, lw=1.0, ls="--", zorder=1)
    ax.set_yticks(ypos)
    ax.set_yticklabels([r["variant"] for r in fp], fontsize=7)
    ax.set_xlabel("Hodges-Lehmann shift in C:N, fish \u2212 mammal (dimensionless)\n"
                  "negative = fish lower; dashed line = no difference")
    ax.set_title("Fish-lower-than-mammal C:N is stable across robustness variants")
    fig.savefig(os.path.join(FIG_DIR, "fig-4-robustness-forest.png"))
    plt.close(fig)

    # ==================================================================== #
    #  write results.json
    # ==================================================================== #
    with open(os.path.join(HERE, "results.json"), "w") as fh:
        json.dump(results, fh, indent=2, sort_keys=True)

    return results


if __name__ == "__main__":
    res = main()
    t1 = res["T1_fish_vs_mammal_CN"]
    print("T1 fish vs mammal C:N: U=%.1f p=%.3g delta=%.3f HL=%.4f CI[%.4f,%.4f]" % (
        t1["U"], t1["p"], t1["cliffs_delta"],
        t1["hodges_lehmann_shift_fish_minus_mammal"], *t1["hl_ci95"]))
    t2 = res["T2_CN_vs_glycine"]["overall"]
    print("T2 C:N~Gly overall: rho=%.3f p=%.3g n=%d" % (t2["rho"], t2["p"], t2["n"]))
    gly = res["T3_aa_screen"]["glycine"]
    print("T3 glycine: rho=%.3f q=%.3g n=%d" % (gly["rho"], gly["q_bh"], gly["n"]))