"""Publication figures. Fully opaque rendering: no alpha anywhere, solid white
canvas, solid tint fills, crisp hairline grid."""
import json, os
import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib.colors import to_rgb
from matplotlib.patches import FancyBboxPatch, FancyArrowPatch, Circle, Rectangle
from matplotlib.lines import Line2D
import ministats as stats
from data import MEANS, DAYS, TRT, PL_FRAC, R

M = json.load(open("models.json"))

WHITE = "#FFFFFF"
GRID = "#E4E7EA"
INK = "#1A1D21"
plt.rcParams.update({
    "font.family": "DejaVu Sans", "font.size": 8.0,
    "axes.linewidth": .8, "axes.edgecolor": "#4A5158",
    "axes.facecolor": WHITE, "figure.facecolor": WHITE, "savefig.facecolor": WHITE,
    "savefig.edgecolor": WHITE, "savefig.transparent": False,
    "xtick.major.width": .8, "ytick.major.width": .8,
    "xtick.major.size": 2.6, "ytick.major.size": 2.6,
    "xtick.color": "#3A4046", "ytick.color": "#3A4046",
    "axes.labelsize": 8.5, "axes.labelcolor": INK,
    "axes.titlesize": 9.0, "text.color": INK,
    "legend.fontsize": 8.0, "legend.frameon": False,
    "figure.dpi": 500, "savefig.dpi": 500, "savefig.bbox": "tight",
    "axes.spines.top": False, "axes.spines.right": False,
    "axes.grid": False, "grid.color": GRID, "grid.linewidth": .6,
    "grid.alpha": 1.0, "path.simplify": False,
    "pdf.fonttype": 42, "ps.fonttype": 42,   # TrueType, not Type 3
})


def tint(c, f):
    """Blend color c toward white by fraction f (0 = pure, 1 = white). Opaque."""
    r, g, b = to_rgb(c)
    return (r + (1 - r) * f, g + (1 - g) * f, b + (1 - b) * f)


BLUE, RED, GREEN, GRAY = "#1F5FA0", "#B02A2F", "#2E7D32", "#6E767D"
CT = ["#0E3450", "#1C6E86", "#3F9A8E", "#D98B3A", "#B02A2F"]
MK = ["o", "s", "^", "D", "v"]
LB = ["100:0", "75:25", "50:50", "25:75", "0:100"]
BAND_G = tint(GREEN, .90)          # solid pale green
BAND_K = tint(GRAY, .88)           # solid pale gray
MATRIX, TOXIC = ["Zn", "Cu", "Fe", "Mn", "Ni"], ["Cr", "Cd", "Pb"]
ALL = MATRIX + TOXIC
D = np.array(DAYS, float); DT = D - D[0]
U = "mg kg$^{-1}$"

MSE = {"Zn": [143.851, 164.889, 195.842, 124.354, 172.220, 188.949],
       "Cu": [124.7027, 76.738, 45.843, 66.958, 20.814, 17.637],
       "Fe": [8.284, 13.772, 5.580, 73.989, 5.500, 10.989],
       "Mn": [18.870, 50.173, 47.231, 7997.871, 40.185, 56.397],
       "Ni": [0.992, 0.405, 0.433, 1.297, 2.389, 1.534]}
SE = {k: np.sqrt(np.array(v) / R) for k, v in MSE.items()}
EF = lambda m, i: MEANS[m][i][5] / MEANS[m][i][0]
# solid lightness ramp for the five treatments within a group color
RAMP = [.46, .35, .24, .12, .0]


def panel(ax, tag, title):
    ax.set_title(f"({tag}) {title}", loc="left", fontweight="bold", pad=5, color=INK)


def grid_y(ax):
    """No background rules: journal figures are plain white."""
    ax.set_axisbelow(True)
    ax.grid(False)


def grid_xy(ax):
    ax.set_axisbelow(True)
    ax.grid(False)


def save(fig, name):
    """Write a 600-dpi raster for Word and a true vector PDF for the typeset PDF."""
    fig.savefig(name, dpi=600, facecolor=WHITE, edgecolor="none", transparent=False)
    fig.savefig(name.replace(".png", ".pdf"), facecolor=WHITE, edgecolor="none",
                transparent=False)
    plt.close(fig)


# ═════════════════════════ FIG 1 concentration trajectories ═════════════════
fig, axes = plt.subplots(4, 2, figsize=(6.85, 8.60))
for k, m in enumerate(ALL):
    ax = axes.flat[k]
    for i in range(5):
        y = np.array(MEANS[m][i])
        kw = dict(color=CT[i], marker=MK[i], ms=3.2, lw=1.25, mew=.6,
                  mec=WHITE, mfc=CT[i], zorder=3 + i)
        if m in SE:
            ax.errorbar(D, y, yerr=SE[m], capsize=2, elinewidth=.75,
                        ecolor=CT[i], label=LB[i], **kw)
        else:
            ax.plot(D, y, ls=(0, (4, 2)), label=LB[i], **kw)
    col = "#D98B3A" if m == "Fe" else BLUE if m in MATRIX else RED
    cls = ("mineral tracer" if m == "Fe"
           else "conservative" if m in MATRIX else "non-conservative")
    ax.set_title(f"({chr(97+k)}) {m}", loc="left", fontweight="bold", pad=5,
                 color=INK)
    ax.set_title(cls, loc="right", fontsize=8, style="italic", color=col, pad=5)
    ax.set_xticks(DAYS); ax.set_xlim(8, 97)
    grid_xy(ax)
    if k >= 6:                       # bottom row of the 4 x 2 grid
        ax.set_xlabel("Composting time (d)")
    if k % 2 == 0:                   # left column
        ax.set_ylabel(f"Concentration ({U})")
h, l = axes.flat[0].get_legend_handles_labels()
fig.legend(h, l, ncol=5, loc="upper center", bbox_to_anchor=(.5, 1.025),
           handlelength=1.7, columnspacing=1.4, handletextpad=.5,
           title="PL: VW (w/w)", title_fontsize=8)
fig.tight_layout(rect=[0, 0, 1, .955], h_pad=1.4, w_pad=1.8)
save(fig, "F1.png")

# ═════════════════════════ FIG 2 conservation test ══════════════════════════
fig = plt.figure(figsize=(6.85, 5.60))
gs = fig.add_gridspec(2, 2, height_ratios=[1.05, 1.0],
                      width_ratios=[.82, 1.0], wspace=.42, hspace=.40)

ax = fig.add_subplot(gs[0, :])
ax.set_ylim(.85, 46); ax.set_yscale("log")
for j, m in enumerate(ALL):
    base = BLUE if m in MATRIX else RED
    for i in range(5):
        ax.scatter(j + (i - 2) * .118, EF(m, i), s=29, color=tint(base, RAMP[i]),
                   ec=WHITE, lw=.6, zorder=3)
    ax.hlines(np.mean([EF(m, i) for i in range(5)]), j - .35, j + .35,
              color=base, lw=2.6, zorder=4)
ax.axhline(2.0, color=GREEN, lw=1.3, ls=(0, (5, 2)), zorder=2)
ax.axhline(1.0, color=GRAY, lw=.9, ls=(0, (1, 2)), zorder=2)
ax.text(7.45, 2.14, "ceiling at 50% dry-matter loss", fontsize=8, color=GREEN,
        ha="right", va="bottom", style="italic")
ax.set_yticks([1, 2, 5, 10, 20, 40]); ax.set_yticklabels("1 2 5 10 20 40".split())
ax.set_xticks(range(8)); ax.set_xticklabels(ALL); ax.set_xlim(-.6, 7.6)
ax.set_ylabel("Enrichment factor,  EF")
panel(ax, "a", "Apparent enrichment, day 15 to day 90")
grid_y(ax)
ax.legend(handles=[Line2D([], [], marker="o", ls="", color=BLUE, ms=5,
                          label="conservative set"),
                   Line2D([], [], marker="o", ls="", color=RED, ms=5,
                          label="Cr, Cd, Pb")], loc="upper left")

ax = fig.add_subplot(gs[1, 0])
mat = np.array([EF(m, i) for m in MATRIX if m != "Fe" for i in range(3)])
tox = np.array([EF(m, i) for m in TOXIC for i in range(3)])
Uu, pp = stats.mannwhitneyu(mat, tox)
ax.set_ylim(.85, 30); ax.set_yscale("log")
pr = ax.violinplot([mat, tox], showextrema=False, widths=.8)
for pc, c in zip(pr["bodies"], [BLUE, RED]):
    pc.set_facecolor(tint(c, .80)); pc.set_alpha(1.0)
    pc.set_edgecolor(c); pc.set_linewidth(1.1); pc.set_zorder(2)
rng = np.random.default_rng(1)
for q, (d, c) in enumerate(zip([mat, tox], [BLUE, RED])):
    ax.scatter(np.full(len(d), q + 1) + rng.uniform(-.07, .07, len(d)), d, s=23,
               color=c, ec=WHITE, lw=.6, zorder=4)
    ax.hlines(np.median(d), q + .70, q + 1.30, color=c, lw=2.6, zorder=5)
ax.axhline(2.0, color=GREEN, lw=1.3, ls=(0, (5, 2)), zorder=3)
ax.set_yticks([1, 2, 5, 10, 20]); ax.set_yticklabels("1 2 5 10 20".split())
ax.set_xticks([1, 2]); ax.set_xticklabels(["conservative", "Cr, Cd, Pb"])
ax.set_ylabel("Enrichment factor,  EF")
panel(ax, "b", "Group contrast")
ax.text(.60, 7.0, f"U = {Uu:.0f}\nP = 1.4 \u00d7 10$^{{-4}}$", ha="left", va="center", fontsize=8,
        bbox=dict(fc=WHITE, ec="#C7CDD2", lw=.6, pad=2.6), zorder=6)
grid_y(ax)

ax = fig.add_subplot(gs[1, 1])
req = sorted(((m, (1 - 1 / np.mean([EF(m, i) for i in range(3)])) * 100)
              for m in ALL), key=lambda z: z[1])
ypos = np.arange(len(req))
cols = [BLUE if m in MATRIX else RED for m, _ in req]
ax.set_xlim(0, 105)
ax.barh(ypos, [v for _, v in req], color=cols, height=.62, ec=WHITE, lw=.7, zorder=3)
ax.axvline(50, color=INK, lw=1.3, ls=(0, (5, 2)), zorder=4)
ax.text(50, len(req) - .30, "not physically attainable", fontsize=8,
        color="#3A4046", ha="center", va="bottom", zorder=6,
        bbox=dict(fc=WHITE, ec="none", pad=1.2))
for q, (m, v) in enumerate(req):
    ax.text(v + 1.8, q, f"{v:.0f}%", va="center", fontsize=8,
            color=cols[q], fontweight="bold", zorder=5)
ax.set_yticks(ypos); ax.set_yticklabels([m for m, _ in req])
ax.set_xticks([0, 20, 40, 60, 80, 100])
ax.set_xlabel("Dry-matter loss required to\nexplain the observed EF (%)")
panel(ax, "c", "Mass-loss ceiling test")
ax.set_axisbelow(True); ax.grid(False)
ax.set_ylim(-.6, len(req) + .35)
save(fig, "F2.png")

# ═════════════════════════ FIG 3 liberation factor ══════════════════════════
SC = list(M["lambda_scenarios"])
fig, axes = plt.subplots(2, 2, figsize=(6.85, 6.40),
                         gridspec_kw={"wspace": .30, "hspace": .34})
axes = axes.flat
for s, nm in enumerate(SC):
    ax = axes[s]
    lam = M["lambda_scenarios"][nm]
    dm = (1 - np.exp(-lam * DT[-1])) * 100
    ax.set_ylim(.35, 40); ax.set_yscale("log")
    for j, m in enumerate(ALL):
        base = BLUE if m in MATRIX else RED
        v = M["LIB"][nm][m]
        for i, x in enumerate(v):
            ax.scatter(j + (i - 2) * .12, x, s=21, color=tint(base, RAMP[i]),
                       ec=WHITE, lw=.5, zorder=3)
        ax.hlines(np.mean(v), j - .34, j + .34, color=base, lw=2.4, zorder=4)
    ax.axhline(1.0, color=INK, lw=1.2, ls=(0, (5, 2)), zorder=2)
    ax.set_yticks([.5, 1, 2, 5, 10, 20])
    ax.set_yticklabels("0.5 1 2 5 10 20".split())
    ax.set_xticks(range(8)); ax.set_xticklabels(ALL, fontsize=8)
    ax.set_xlim(-.6, 7.6)
    panel(ax, chr(97 + s), nm.split(None, 1)[1].replace(" %", "%"))
    ax.text(.5, .965, f"$\\lambda$ = {lam:.4f} d$^{{-1}}$   ({dm:.0f}% DM loss)",
            transform=ax.transAxes, ha="center", va="top", fontsize=8,
            color="#3A4046")
    if s == 0:
        ax.set_ylabel("Liberation factor,  $\\Lambda$")
    grid_y(ax)
axes[0].text(.035, .045, "conservative\nbehavior", transform=axes[0].transAxes,
             fontsize=8, color="#3A4046", va="bottom", style="italic")

ax = axes[3]
xs = np.arange(len(SC))
mv = [M["lambda_test"][nm]["matrix"] for nm in SC]
tv = [M["lambda_test"][nm]["toxic"] for nm in SC]
ax.set_ylim(0, 6.6)
ax.plot(xs, tv, "-s", color=RED, lw=1.9, ms=6.5, mec=WHITE, mew=.9,
        label="Cr, Cd, Pb", zorder=4)
ax.plot(xs, mv, "-o", color=BLUE, lw=1.9, ms=6.5, mec=WHITE, mew=.9,
        label="conservative set", zorder=4)
ax.axhline(1.0, color=INK, lw=1.2, ls=(0, (5, 2)), zorder=2)
for q in range(len(SC)):
    ax.annotate(f"{tv[q]:.2f}", (q, tv[q]), textcoords="offset points",
                xytext=(0, 9), ha="center", fontsize=8, color=RED,
                fontweight="bold", zorder=6, bbox=dict(fc=WHITE, ec="none", pad=1.2))
    ax.annotate(f"{mv[q]:.2f}", (q, mv[q]), textcoords="offset points",
                xytext=(0, -13), ha="center", fontsize=8, color=BLUE,
                fontweight="bold", zorder=6, bbox=dict(fc=WHITE, ec="none", pad=1.2))
ax.set_xticks(xs); ax.set_xlim(-.28, 2.28)
ax.set_xticklabels(["tracer\n(8% loss)", "30%\nloss", "50%\nloss"], fontsize=8)
ax.set_ylabel("Mean liberation factor,  $\\Lambda$")
panel(ax, "d", "Robustness to the mass-loss assumption")
ax.legend(loc="upper right", bbox_to_anchor=(1.0, .90))
grid_y(ax)
# short band tag only; the full statement is carried by the caption
ax.text(2.22, 1.10, "inadmissible below $\\Lambda$ = 1", ha="right", va="bottom",
        fontsize=8, color="#3A4046", style="italic", zorder=6)
save(fig, "F3.png")

# ═════════════════════════ FIG 4 kinetics ═══════════════════════════════════
fig, axes = plt.subplots(2, 2, figsize=(6.85, 5.90),
                         gridspec_kw={"wspace": .30, "hspace": .34})
axes = axes.flat
tt = np.linspace(0, DT[-1], 240)
KC = [CT[0], CT[1], CT[2]]
for s, m in enumerate(TOXIC):
    ax = axes[s]
    for i in range(3):
        d = M["kinetics"][m][TRT[i]]
        Li, k = d["Linf"], d["k"]
        ax.plot(tt, Li - (Li - 1) * np.exp(-k * tt), color=KC[i], lw=1.5, zorder=3)
        ax.scatter(DT, d["traj"], s=28, color=KC[i], ec=WHITE, lw=.7,
                   zorder=4, label=LB[i])
    ax.axhline(1.0, color=INK, lw=1.1, ls=(0, (5, 2)), zorder=2)
    panel(ax, chr(97 + s), m)
    ax.set_xlabel("Time since first sampling (d)")
    if s == 0:
        ax.set_ylabel("Liberation factor,  $\\Lambda$(t)")
    kk = np.mean([M["kinetics"][m][TRT[i]]["k"] for i in range(3)])
    ax.text(.045, .95, f"$k$ = {kk:.3f} d$^{{-1}}$\n$t_{{1/2}}$ = {np.log(2)/kk:.0f} d",
            transform=ax.transAxes, va="top", fontsize=8,
            bbox=dict(fc=WHITE, ec="#C7CDD2", lw=.6, pad=2.4), zorder=6)
    grid_xy(ax)

ax = axes[3]
for s, (m, c) in enumerate(zip(["Cu", "Ni"], [BLUE, "#3F9A8E"])):
    for i in range(3):
        d = M["kinetics"][m][TRT[i]]
        ax.scatter(DT, d["traj"], s=24, color=c, ec=WHITE, lw=.6, zorder=3,
                   label=m if i == 0 else None)
ax.axhline(1.0, color=INK, lw=1.1, ls=(0, (5, 2)), zorder=2)
ax.set_ylim(.55, 1.38)
panel(ax, "d", "Null case: Cu and Ni")
ax.set_xlabel("Time since first sampling (d)")
ax.set_ylabel("Liberation factor,  $\\Lambda$(t)")
ax.text(.5, .085, "no liberation detected\n($\\Lambda_\\infty$ = 1.00)",
        transform=ax.transAxes, ha="center", fontsize=8, color="#3A4046",
        style="italic")
ax.legend(loc="upper right")
grid_xy(ax)
h4, l4 = axes[0].get_legend_handles_labels()
fig.legend(h4, l4, ncol=3, loc="upper center", bbox_to_anchor=(.5, 1.005),
           handletextpad=.4, columnspacing=1.8, title="PL: VW", title_fontsize=8)
fig.tight_layout(rect=[0, 0, 1, .935], h_pad=1.6, w_pad=1.8)
save(fig, "F4.png")

# ═══════ FIG 7 conceptual model, rendered in isometric projection ══════════
# The figure carries the whole argument of the paper in one image: the metal
# inventory of the pile is fixed, the matrix barely shrinks, but the fraction
# the digestion can recover climbs steeply, so the reported concentration
# rises and the compliance verdict for an unchanged pile inverts.
from matplotlib.patches import Polygon

LAM = M["lambda_scenarios"]["S1 tracer-based"]
COS30, SIN30 = np.cos(np.pi / 6), np.sin(np.pi / 6)
A_ = B_ = 13.0            # footprint of the heap, isometric units
C_ = 9.4                  # height of the heap
DIN = 3.5                 # inset of the top face, giving a heap its batter
FRESH, HUMUS = "#CBAC72", "#5B4830"


def iso(p, ox, oy, s):
    """Project a 3-D point onto the page in a 30 degree isometric view."""
    x, y, z = p
    return (ox + (x - y) * COS30 * s, oy + (x + y) * SIN30 * s + z * s)


def dark(c, f):
    r, g, b = to_rgb(c)
    return (r * (1 - f), g * (1 - f), b * (1 - f))


def lerp(c1, c2, t):
    a_, b_ = np.array(to_rgb(c1)), np.array(to_rgb(c2))
    return tuple(a_ + (b_ - a_) * t)


def heap(ax, ox, oy, s, base, nfib, ngran, nocc, nlib, seed):
    """Draw one composting heap as a solid isometric frustum, textured to show
    its degree of humification and carrying a fixed number of metal units."""
    rng = np.random.default_rng(seed)
    top, rgt, lft = lerp(base, WHITE, .30), base, dark(base, .28)
    P = lambda p: iso(p, ox, oy, s)
    b00, b10, b01 = (0, 0, 0), (A_, 0, 0), (0, B_, 0)
    t00, t10, t11, t01 = ((DIN, DIN, C_), (A_ - DIN, DIN, C_),
                          (A_ - DIN, B_ - DIN, C_), (DIN, B_ - DIN, C_))
    for face, col, zo in (([b00, b01, t01, t00], lft, 2),
                          ([b00, b10, t10, t00], rgt, 2),
                          ([t00, t10, t11, t01], top, 3)):
        ax.add_patch(Polygon([P(p) for p in face], closed=True, fc=col,
                             ec=dark(base, .48), lw=.9, zorder=zo,
                             joinstyle="round"))

    def pts(n):
        """Scatter n points over the three visible faces, in proportion to the
        area each presents to the viewer."""
        out = []
        for f in rng.choice(["R", "L", "T"], size=n, p=[.37, .34, .29]):
            u = rng.uniform(.12, .86)
            if f == "R":
                out.append((rng.uniform(DIN * u + 1.0, A_ - DIN * u - 1.0),
                            0.0, u * C_))
            elif f == "L":
                out.append((0.0, rng.uniform(DIN * u + 1.0, B_ - DIN * u - 1.0),
                            u * C_))
            else:
                out.append((rng.uniform(DIN + 1.1, A_ - DIN - 1.1),
                            rng.uniform(DIN + 1.1, B_ - DIN - 1.1), C_))
        return out

    for p in pts(nfib):                                   # lignocellulose
        X, Y = P(p); a = rng.uniform(0, np.pi)
        ax.plot([X - .80 * s * np.cos(a), X + .80 * s * np.cos(a)],
                [Y - .80 * s * np.sin(a), Y + .80 * s * np.sin(a)],
                color=dark(base, .34), lw=1.4 * s, solid_capstyle="round",
                zorder=4)
    for p in pts(ngran):                                  # humic granules
        X, Y = P(p)
        ax.add_patch(Circle((X, Y), rng.uniform(.36, .62) * s,
                            fc=lerp(base, WHITE, .40), ec="none", zorder=4))
    for p in pts(nocc):                                   # occluded metal
        X, Y = P(p)
        ax.add_patch(Circle((X, Y), .60 * s, fc="#F3D9D9", ec=RED, lw=1.15,
                            zorder=6))
    for p in pts(nlib):                                   # acid-extractable
        X, Y = P(p)
        ax.add_patch(Circle((X, Y), .55 * s, fc=RED, ec=WHITE, lw=.55, zorder=7))
        ax.add_patch(Circle((X - .18 * s, Y + .18 * s), .16 * s, fc="#F6C9C9",
                            ec="none", zorder=8))


def badge(ax, x, y, n, ms=10.0):
    """A numbered marker. Drawn as a plot marker so that it stays circular
    whatever the aspect ratio of the host axes."""
    ax.plot([x], [y], marker="o", ms=ms, mfc=INK, mec="none", ls="",
            zorder=9, clip_on=False)
    ax.text(x, y, str(n), ha="center", va="center", color=WHITE,
            fontsize=8, fontweight="bold", zorder=10, clip_on=False)


fig = plt.figure(figsize=(6.85, 8.90))
gs = fig.add_gridspec(2, 1, height_ratios=[1.7, 1.0], hspace=.12,
                      left=.055, right=.985, top=.985, bottom=.075)
axT = fig.add_subplot(gs[0]); axB = fig.add_subplot(gs[1])
# the isometric projection only reads correctly at unit aspect, so the data
# limits of the upper panel are matched to its physical shape on the page
# the isometric projection only reads correctly at unit aspect, so the data
# limits of the upper panel are matched to its physical shape on the page
axT.set_xlim(0, 100); axT.set_ylim(0, 75.5); axT.axis("off")
axT.set_facecolor(WHITE)

# stage definitions: day, humification progress, matrix retained, metal split
STAGE = [(15, 0.00, 1.000, 14, 4, "1.0"),
         (45, 0.52, 0.967, 9, 9, "2.1"),
         (90, 1.00, 0.919, 1, 17, "3.8")]
CX, BASEY, SC = [18.0, 50.0, 82.0], 33.0, 1.05

for i, ((day, prog, keep, nocc, nlib, lab), cx) in enumerate(zip(STAGE, CX)):
    s = SC * keep ** (1 / 3)                   # true volumetric shrinkage
    heap(axT, cx, BASEY, s, lerp(FRESH, HUMUS, prog),
         nfib=int(round(26 - 22 * prog)), ngran=int(round(2 + 24 * prog)),
         nocc=nocc, nlib=nlib, seed=11 + i)
    badge(axT, cx - 7.5, 66.5, i + 1)
    axT.text(cx - 5.2, 66.5, f"Day {day}", ha="left", va="center",
             fontsize=9.5, fontweight="bold")
    axT.text(cx, 60.5, f"$\\Lambda$ = {lab}", ha="center", fontsize=9,
             color=RED, fontweight="bold")
    axT.text(cx, 29.0, f"{nlib} of 18 extractable", ha="center", fontsize=8,
             color="#3A4046")

for x0, x1 in ((31.5, 36.5), (63.5, 68.5)):
    axT.annotate("", xy=(x1, 44), xytext=(x0, 44),
                 arrowprops=dict(arrowstyle="-|>", lw=1.9, color=INK,
                                 mutation_scale=15, shrinkA=0, shrinkB=0))

# the two competing rate processes, stated once beneath the sequence
axT.add_patch(FancyBboxPatch((5.0, 19.5), 90.0, 5.4,
                             boxstyle="round,pad=0.4,rounding_size=1.2",
                             fc="#F4F6F8", ec="#C3CDD5", lw=.9, zorder=1))
axT.text(28.0, 22.2, "humification  $k$ = 0.0143 d$^{-1}$",
         ha="center", va="center", fontsize=8, color=RED, zorder=3)
axT.text(72.0, 22.2, "matrix loss  $\\lambda$ = 0.0011 d$^{-1}$",
         ha="center", va="center", fontsize=8, color=BLUE, zorder=3)

# the governing identity, set across the full width beneath the sequence
axT.add_patch(Rectangle((5.0, 2.0), 90.0, 14.5, fc=WHITE, ec="#C3CDD5",
                        lw=1.0, zorder=1))
axT.text(20.0, 9.2, "$C(t)=\\dfrac{M\\,f(t)}{m(t)}$", ha="center",
         va="center", fontsize=12, zorder=3)
axT.plot([37.0, 37.0], [3.5, 15.0], color="#C3CDD5", lw=1.0, zorder=3)
axT.text(45.0, 12.2, "$M$ fixed", ha="left", va="center", fontsize=8,
         color="#3A4046", zorder=3)
axT.text(45.0, 9.2, "$f\\uparrow$  ×3.8", ha="left", va="center", fontsize=8,
         color=RED, zorder=3)
axT.text(45.0, 6.2, "$m\\downarrow$  −8%", ha="left", va="center", fontsize=8,
         color=BLUE, zorder=3)
axT.plot([66.0, 66.0], [3.5, 15.0], color="#C3CDD5", lw=1.0, zorder=3)
axT.text(80.5, 9.2, "$C$ rises fourfold with\nno gain in metal content",
         ha="center", va="center", fontsize=8, fontweight="bold", zorder=3)

axT.legend(handles=[
    Line2D([], [], marker="o", ls="", mfc=RED, mec=WHITE, ms=7,
           label="acid-extractable metal"),
    Line2D([], [], marker="o", ls="", mfc="#F3D9D9", mec=RED, mew=1.3, ms=7.5,
           label="occluded metal"),
    Line2D([], [], color=dark(FRESH, .34), lw=2.2, label="lignocellulose"),
    Line2D([], [], marker="o", ls="", mfc=lerp(HUMUS, WHITE, .40), mec="none",
           ms=7, label="humic matrix")],
    loc="lower center", ncol=4, bbox_to_anchor=(.5, -.055),
    handletextpad=.4, columnspacing=1.2, fontsize=8)

# ── lower panel: the same story told by the measured data ──────────────────
cd = np.array(MEANS["Cd"][0], float)
dd = np.array(DAYS, float)
cons = cd[0] * np.exp(LAM * (dd - dd[0]))
axB.set_xlim(9, 96); axB.set_ylim(0, 4.2)
axB.add_patch(Rectangle((9, 0), 87, 1.5, fc=tint(GREEN, .90), ec="none",
                        zorder=0))
axB.add_patch(Rectangle((9, 1.5), 87, 2.7, fc=tint(RED, .91), ec="none",
                        zorder=0))
axB.axhline(1.5, color=INK, lw=1.3, ls=(0, (6, 3)), zorder=4)
axB.plot(dd, cons, ls=(0, (4, 2.5)), color=BLUE, lw=1.9, zorder=5,
         label="expected if Cd were conservative (matrix loss only)")
axB.plot(dd, cd, "-o", color=RED, lw=2.1, ms=5.4, mfc=RED, mec=WHITE, mew=.9,
         zorder=6, label="measured Cd, 100:0 pile")
xc = 30 + 15 * (1.5 - cd[1]) / (cd[2] - cd[1])
axB.plot([xc, xc], [0, 1.5], color=INK, lw=.9, ls=":", zorder=5)
axB.annotate(f"verdict inverts at {xc:.0f} d", xy=(xc, 1.5),
             xytext=(xc + 4.5, .46), fontsize=8, fontweight="bold",
             ha="left", va="center",
             arrowprops=dict(arrowstyle="-|>", lw=.9, color=INK,
                             mutation_scale=10, shrinkA=1, shrinkB=2))
axB.text(94.5, 1.36, "compliant", ha="right", va="top", fontsize=8,
         color=dark(GREEN, .15), fontweight="bold", zorder=5)
axB.text(94.5, 1.64, "above the EU limit", ha="right", va="bottom",
         fontsize=8, color=dark(RED, .10), fontweight="bold", zorder=5)
axB.text(10.5, 1.62, "EU 2019/1009 limit, 1.5 mg kg$^{-1}$", ha="left",
         va="bottom", fontsize=8, color="#3A4046", zorder=5)
for i, st in enumerate(STAGE):
    j = list(dd).index(st[0])
    badge(axB, st[0], cd[j] + .34, i + 1, ms=9.0)
axB.set_xlabel("Composting day"); axB.set_ylabel("Cd (" + U + ")")
axB.set_xticks(list(dd))
axB.legend(loc="upper left", bbox_to_anchor=(.005, .99), fontsize=8)
grid_xy(axB)
save(fig, "F7.png")   # conceptual synthesis, cited last

# ═════════════════════════ FIG 6 risk models ════════════════════════════════
fig = plt.figure(figsize=(6.85, 6.30))
gs6 = fig.add_gridspec(2, 2, height_ratios=[1.0, 1.05], wspace=.40, hspace=.46)

ax = fig.add_subplot(gs6[0, 0])
g = np.array(M["mc"]["grid"]) * 100; pe = np.array(M["mc"]["p_exceed"])
ax.set_xlim(0, 100); ax.set_ylim(0, 1.02)
ax.plot(g, pe, color=RED, lw=2.2, zorder=4)
for tgt, c in [(.05, GREEN), (.25, "#D98B3A"), (.50, RED)]:
    # the share at which the risk curve crosses the tolerance, interpolated
    # so that the marker sits exactly on the curve rather than on the
    # nearest grid node, whose probability differs from the target
    xa = float(np.interp(tgt, pe, g))
    ax.plot([0, xa], [tgt, tgt], color=c, lw=1.0, ls=(0, (2, 2)), zorder=3)
    ax.plot([xa, xa], [0, tgt], color=c, lw=1.0, ls=(0, (2, 2)), zorder=3)
    ax.scatter([xa], [tgt], s=30, color=c, ec=WHITE, lw=.8, zorder=5)
    ax.text(xa + 3.0, tgt + .035, f"{xa:.0f}%", fontsize=8, color=c,
            fontweight="bold", zorder=6,
            bbox=dict(fc=WHITE, ec="none", pad=1.0))
ax.set_xlabel("PL in feedstock (%, w/w)")
ax.set_ylabel("P (Cd exceeds 1.5 mg kg$^{-1}$)")
panel(ax, "a", "Monte-Carlo compliance risk")
grid_xy(ax)

ax = fig.add_subplot(gs6[0, 1])
yrs = np.arange(0, 121)
SOIL, BG, LIM = 2600., .41, 1.0
ax.set_xlim(0, 120); ax.set_ylim(.14, 1.30)
for i in range(5):
    ax.plot(yrs, BG + 10 * MEANS["Cd"][i][5] / SOIL * yrs, color=CT[i], lw=1.9,
            label=LB[i], zorder=3)
ax.axhline(LIM, color=RED, lw=1.4, ls=(0, (5, 2)), zorder=2)
ax.text(118, LIM - .030, "EU soil limit, 1.0 mg kg$^{-1}$", fontsize=8,
        color=RED, ha="right", va="top", zorder=6,
        bbox=dict(fc=WHITE, ec="none", pad=1.0))
ax.axhline(BG, color=GRAY, lw=.9, ls=(0, (1, 2)), zorder=2)
ax.text(4, BG - .045, "background", fontsize=8, color=GRAY, ha="left",
        va="top")
for i in range(3):
    ya = (LIM - BG) / (10 * MEANS["Cd"][i][5] / SOIL)
    ax.scatter([ya], [LIM], s=30, color=CT[i], ec=WHITE, lw=.8, zorder=5)
    ax.text(ya, LIM + [.16, .095, .035][i], f"{ya:.0f} yr", fontsize=8,
            color=CT[i], ha="center", fontweight="bold", zorder=6,
            bbox=dict(fc=WHITE, ec="none", pad=1.0))
ax.set_xlabel("Years of annual application")
ax.set_ylabel("Plow-layer Cd (mg kg$^{-1}$)")
panel(ax, "b", "Soil accumulation")
ax.legend(loc="lower right", ncol=1, fontsize=8, title="PL: VW",
          title_fontsize=8, bbox_to_anchor=(1.0, .02), labelspacing=.28,
          handlelength=1.4, borderpad=.35, frameon=True, framealpha=1.0,
          facecolor=WHITE, edgecolor="#C7CDD2")
grid_xy(ax)

ax = fig.add_subplot(gs6[1, :])
mets = ["Zn", "Cr", "Cu", "Ni", "Pb", "Cd"]
TRf = {"Zn": 1, "Cr": 2, "Cu": 5, "Ni": 5, "Pb": 5, "Cd": 30}
cols = ["#8FB8D0", "#7FAE94", "#DCC08A", "#C2A0C2", "#A3AEB8", RED]
bot = np.zeros(5)
for m, c in zip(mets, cols):
    v = np.array([TRf[m] * M["indices"][t]["CF"][m] for t in TRT])
    ax.bar(range(5), v, .66, bottom=bot, color=c, ec=WHITE, lw=.7, label=m,
           zorder=3)
    bot += v
for q, t in enumerate(TRT):
    ax.text(q, bot[q] + 10, f"{M['indices'][t]['RI']:.0f}", ha="center",
            fontsize=8, fontweight="bold", zorder=6, bbox=dict(fc=WHITE, ec="none", pad=1.2))
ax.axhline(150, color="#D98B3A", lw=1.3, ls=(0, (5, 2)), zorder=2)
ax.text(3.45, 150, "moderate risk", fontsize=8, color="#D98B3A", ha="center",
        va="center", zorder=6, bbox=dict(fc=WHITE, ec="none", pad=1.6))
ax.set_xticks(range(5))
ax.set_xticklabels(LB, fontsize=8)
ax.set_xlabel("PL: VW")
ax.set_ylabel("Potential ecological risk index,  RI")
panel(ax, "c", "Ecological risk contributions")
ax.legend(ncol=6, loc="upper right", fontsize=8, columnspacing=1.0,
          handlelength=1.2, handletextpad=.5)
ax.set_ylim(0, 400)
grid_y(ax)
save(fig, "F6.png")

for f in sorted(os.listdir(".")):
    if f.startswith("F") and f.endswith(".png") and len(f) <= 7:
        print(f"  {f}  {os.path.getsize(f)/1024:.0f} kB")


# ═════════════════════════ FIG 5 ordination and model support ═══════════════
A = json.load(open("addstats.json"))
fig, axes = plt.subplots(1, 3, figsize=(6.85, 3.30),
                         gridspec_kw={"width_ratios": [1, 1, 1.15], "wspace": .34})

# (a) sample scores -------------------------------------------------------
ax = axes[0]
sc = np.array(A["pca"]["scores"])
ex = A["pca"]["explained"]
ax.axhline(0, color="#C3CAD0", lw=.8, zorder=1)
ax.axvline(0, color="#C3CAD0", lw=.8, zorder=1)
ax.plot(sc[:, 0], sc[:, 1], "-", color="#C3CAD0", lw=1.0, zorder=2)
for i2 in range(5):
    ax.scatter(sc[i2, 0], sc[i2, 1], s=70, color=CT[i2], ec=WHITE, lw=1.0,
               marker=MK[i2], zorder=5, label=LB[i2])
    ax.annotate(f"T{i2+1}", (sc[i2, 0], sc[i2, 1]), textcoords="offset points",
                xytext=(0, 9), ha="center", fontsize=8, fontweight="bold",
                color=CT[i2], zorder=6, bbox=dict(fc=WHITE, ec="none", pad=1.2))
ax.set_xlabel(f"PC1 ({ex[0]:.1f}% of variance)")
ax.set_ylabel(f"PC2 ({ex[1]:.1f}%)")
ax.set_ylim(-1.55, 2.05)
panel(ax, "a", "Sample ordination")
ax.text(.97, .04, f"PC1 vs PL fraction\n$r$ = {A['pca']['r_pc1']:+.3f}, "
        f"$P$ = {A['pca']['p_pc1']:.4f}", transform=ax.transAxes, va="bottom",
        ha="right", fontsize=8,
        bbox=dict(fc=WHITE, ec="#C7CDD2", lw=.6, pad=2.4), zorder=7)
grid_xy(ax)   # treatments are annotated directly; no legend needed

# (b) variable loadings ---------------------------------------------------
ax = axes[1]
load = np.array(A["pca"]["loadings"])
y = np.arange(len(ALL))
ax.barh(y - .19, load[:, 0], height=.36, color=[BLUE if m in MATRIX else RED
        for m in ALL], ec=WHITE, lw=.6, zorder=3, label="PC1")
ax.barh(y + .19, load[:, 1], height=.36,
        color=[tint(BLUE if m in MATRIX else RED, .55) for m in ALL],
        ec=WHITE, lw=.6, zorder=3, label="PC2")
ax.axvline(0, color=INK, lw=.9, zorder=4)
ax.set_yticks(y); ax.set_yticklabels(ALL)
ax.invert_yaxis()
ax.set_xlabel("Loading")
panel(ax, "b", "Variable loadings")
ax.legend(loc="lower right", fontsize=8, bbox_to_anchor=(1.0, .02))
grid_y(ax)

# (c) model support -------------------------------------------------------
ax = axes[2]
mc = A["model_comparison"]
order = ["Zn", "Cu", "Mn", "Ni", "Fe", "Pb", "Cr", "Cd"]
for j, m in enumerate(order):
    for i2 in range(3):
        v = mc[m][TRT[i2]]["dAICc"]
        ax.bar(j + (i2 - 1) * .25, v, width=.23,
               color=tint(RED if m in TOXIC else BLUE, [.34, .17, .0][i2]),
               ec=WHITE, lw=.5, zorder=3)
ax.axhline(0, color=INK, lw=.9, zorder=4)
ax.axhline(-4, color=GREEN, lw=1.2, ls=(0, (5, 2)), zorder=4)
ax.text(-.45, -5.6, "strong support for liberation", fontsize=8, color=GREEN,
        ha="left", va="top", style="italic", zorder=6, bbox=dict(fc=WHITE, ec="none", pad=1.2))
ax.set_xticks(range(len(order))); ax.set_xticklabels(order)
ax.set_xlim(-.6, 7.6)
ax.set_ylim(-30, 25)
ax.set_ylabel("$\\Delta$AICc  (liberation − null)")
panel(ax, "c", "Support for the liberation model")
ax.text(.03, .96, "conservative set: null favored", transform=ax.transAxes,
        fontsize=8, color=BLUE, ha="left", va="top", style="italic")
grid_y(ax)
save(fig, "F5.png")   # ordination and model support, cited in Sect. 3.7
