ml_meta_perf.plots

  1"""Figures for the study.
  2
  3matplotlib is imported through the Agg backend so these run headless, in CI and over
  4ssh, without a display. Every function takes a destination path, writes one figure, and
  5returns that path, so the caller decides where output lands and nothing is written as a
  6side effect of importing.
  7
  8**No titles and no annotations.** These figures are written for a LaTeX document where
  9the caption carries the description, and a title rendered into the PNG would duplicate it
 10in a font the document cannot restyle. What stays is the part a caption cannot replace:
 11axis labels, tick labels, and legends identifying the series. Anything a reader would
 12otherwise have to be *told* is instead drawn -- a reference level becomes a line with a
 13legend entry, not a sentence.
 14
 15Each figure carries one message and gets its own caption. The ranking-quality figure uses two
 16aligned axes because its unitless ranking scores and MCC regret must not share a scale.
 17
 18The palette is colour-blind safe and every series is distinguished by marker or line
 19style as well as colour, so the figures survive being printed in greyscale.
 20"""
 21
 22from __future__ import annotations
 23
 24import tempfile
 25import time
 26from pathlib import Path
 27
 28import matplotlib
 29
 30matplotlib.use("Agg")
 31
 32import matplotlib.pyplot as plt
 33import numpy as np
 34import polars as pl
 35from matplotlib.figure import Figure
 36from matplotlib.patches import Rectangle
 37from matplotlib.text import Text
 38
 39IN_SAMPLE = "#1b6ca8"
 40LOO_DATASET = "#d1495b"
 41LOO_MODEL = "#00798c"
 42CEILING = "#6b7280"
 43POSITIVE = "#00798c"
 44NEGATIVE = "#d1495b"
 45
 46FIGURE_DPI = 150
 47
 48# The current paper uses elsarticle's preprint, 12pt layout. In that layout the
 49# standard article class gives \textwidth = 390pt. DejaVu Sans looks larger than
 50# the paper's Latin Modern Roman at the same nominal size; 8.5pt labels align
 51# visually with its body text and captions after \includegraphics scales the PDF.
 52PAPER_FONT_SIZE_PT = 8.5
 53PAPER_TEXT_WIDTH_PT = 390.0
 54
 55
 56#: Written alongside the raster for every figure. A vector copy is what a paper's
 57#: typesetting actually wants, and generating it here costs one extra ``savefig``.
 58VECTOR_SUFFIX = ".pdf"
 59
 60
 61def _save_atomic(figure: Figure, path: Path, *, vector: bool) -> None:
 62    """Render beside ``path`` and replace it, tolerating brief Windows reader locks."""
 63    with tempfile.NamedTemporaryFile(dir=path.parent, prefix=f".{path.name}.", suffix=".tmp", delete=False) as handle:
 64        temporary = Path(handle.name)
 65    try:
 66        if vector:
 67            figure.savefig(
 68                temporary,
 69                format="pdf",
 70                bbox_inches="tight",
 71                transparent=True,
 72                metadata={"CreationDate": None},
 73            )
 74        else:
 75            figure.savefig(
 76                temporary,
 77                format="png",
 78                dpi=FIGURE_DPI,
 79                bbox_inches="tight",
 80                transparent=True,
 81            )
 82        for attempt in range(10):
 83            try:
 84                temporary.replace(path)
 85                return
 86            except OSError:
 87                if attempt == 9:
 88                    raise
 89                time.sleep(0.1)
 90    finally:
 91        temporary.unlink(missing_ok=True)
 92
 93
 94def _match_paper_font(figure: Figure, width_fraction: float, *, tight_layout: bool) -> None:
 95    """Match label size to the body text after scaling to the paper width."""
 96    target_width_pt = PAPER_TEXT_WIDTH_PT * width_fraction
 97    font_size = PAPER_FONT_SIZE_PT
 98    for _ in range(8):
 99        for axes in figure.axes:
100            axes.tick_params(axis="both", labelsize=font_size)
101        for label in figure.findobj(match=Text):
102            label.set_fontsize(font_size)
103        if tight_layout:
104            figure.tight_layout()
105        figure.canvas.draw()
106        bounds = figure.get_tightbbox(figure.canvas.get_renderer())
107        saved_width_pt = 72 * (bounds.width + 2 * plt.rcParams["savefig.pad_inches"])
108        revised = PAPER_FONT_SIZE_PT * saved_width_pt / target_width_pt
109        if abs(revised - font_size) < 0.01:
110            break
111        font_size = revised
112
113
114def _finish(
115    figure: Figure,
116    destination: str | Path,
117    *,
118    paper_width_fraction: float = 0.99,
119    tight_layout: bool = True,
120) -> Path:
121    """Write the figure as PNG and PDF, both on a transparent background.
122
123    Transparent rather than white so a figure sits on whatever the page behind it is,
124    without a rectangle of the wrong shade around it. Text and ticks keep matplotlib's
125    near-black default, which suits the white or near-white page these are written for;
126    on a dark background they would need restyling, and transparency alone would not be
127    enough.
128
129    The returned path is the PNG: it is what the markdown chapters embed. The PDF sits
130    beside it under the same stem, since a paper's typesetting wants the vector copy.
131    """
132    path = Path(destination)
133    path.parent.mkdir(parents=True, exist_ok=True)
134    _match_paper_font(figure, paper_width_fraction, tight_layout=tight_layout)
135    _save_atomic(figure, path, vector=False)
136    # `CreationDate: None` because matplotlib otherwise stamps the PDF with the wall clock,
137    # so every run rewrote seven tracked figures with byte-different, visually identical
138    # files. That makes `git status` dirty after any run and makes "did this change the
139    # output?" -- the check this project verifies optimisations with -- unanswerable for the
140    # vector figures. The PNGs were already deterministic.
141    _save_atomic(figure, path.with_suffix(VECTOR_SUFFIX), vector=True)
142    plt.close(figure)
143    return path
144
145
146#: Model features abbreviated in rendered term labels. The full names are five syllables long
147#: and a full equation figure cannot carry them; the glossary in chapter 1 expands them.
148TERM_ABBREVIATIONS: dict[str, str] = {
149    "Processing Units Number": "PUN",
150    "Model Capability": "MC",
151    "Input Distribution Modelling": "IDM",
152    "Fitting Regime": "FR",
153    "Solution Stochasticity": "SS",
154    "Loss Margin Behaviour": "LMB",
155}
156
157
158#: Division rendered inline, as ``\nicefrac`` would set it, rather than as a built-up
159#: ``\frac``. Mathtext shrinks the two levels of a built-up fraction to fit the line, so a
160#: figure mixing ``a * b`` with ``a / b`` labels drew half its terms at full size and half at
161#: roughly two thirds of it -- see the term-effect figure before 2026-09-07, where
162#: ``log(gravity) / log(PUN)`` was visibly smaller than the product above it. An inline slash
163#: keeps every label on one baseline at one size, which is the whole readability argument.
164DIVIDE = r"\,/\,"
165
166
167def _math_atom(text: str) -> tuple[str, bool]:
168    """One operand of a term as mathtext, and whether it nests without parentheses.
169
170    The flag is what keeps an inline slash unambiguous. ``1/f`` renders as ``1 / f``, so
171    dropping it into ``[1/f] / [g]`` unparenthesised would read ``1 / f / g`` -- a different
172    expression. Atoms that bracket themselves (a parenthesised ``\\log``, a radical, a power,
173    a bare name) need no help; the reciprocal is the only form that does.
174    """
175    text = text.strip()
176    # ``\log(x)``, not ``\log x``: without the parentheses a product reads as
177    # ``log gravity * log MC``, where nothing on the page says how far the first log extends.
178    for pattern, wrap in ((r"log(", r"\log({})"), (r"sqrt(", r"\sqrt{{{}}}")):
179        if text.startswith(pattern) and text.endswith(")"):
180            return wrap.format(_math_atom(text[len(pattern) : -1])[0]), True
181    if text.endswith("^2"):
182        return _math_atom(text[:-2])[0] + "^{2}", True
183    if text.startswith("1/"):
184        return "1" + DIVIDE + _nested(text[2:]), False
185    if all(character.isdigit() or character == "." for character in text):
186        return text, True
187    name = TERM_ABBREVIATIONS.get(text, text)
188    return r"\mathrm{" + name.replace("_", r"\_").replace(" ", r"\ ") + "}", True
189
190
191def _unwrap(text: str) -> str:
192    """Drop the ``[...]`` a term name puts around each operand."""
193    text = text.strip()
194    return text[1:-1] if text.startswith("[") and text.endswith("]") else text
195
196
197def _math_operand(text: str) -> str:
198    """An operand standing alone, where no parentheses are needed."""
199    return _math_atom(_unwrap(text))[0]
200
201
202def _nested(text: str) -> str:
203    """An operand standing inside a larger expression, parenthesised if it has to be."""
204    expression, atomic = _math_atom(_unwrap(text))
205    return expression if atomic else r"\left(" + expression + r"\right)"
206
207
208def _split_top(text: str, operator: str) -> tuple[str, str] | None:
209    """Split on ``operator`` at bracket depth zero, so operators inside `[...]` are ignored."""
210    depth = 0
211    for index, character in enumerate(text):
212        if character in "[(":
213            depth += 1
214        elif character in "])":
215            depth -= 1
216        elif depth == 0 and text.startswith(operator, index):
217            return text[:index], text[index + len(operator) :]
218    return None
219
220
221def term_to_math(name: str) -> str:
222    """A term name as a mathtext expression, set on one line at one size.
223
224    Term names are written for a CSV -- ``[log(gravity)] / [log(Processing Units Number)]`` --
225    and a figure full of those is a wall of brackets. Rendered as mathematics the same
226    term is one expression, which is how it would appear in the paper the equation is for.
227
228    Three conventions, all chosen for a figure whose labels have to be read at a glance
229    beside each other rather than for typographic elegance in isolation:
230
231    * **Division inline, not built up.** ``a / b`` rather than a two-storey ``\\frac``; see
232      ``DIVIDE``. Every label then has the same x-height, which a mixed set does not.
233    * **Logarithms parenthesised.** ``\\log(a) \\times \\log(b)``, never ``\\log a \\log b``,
234      which does not say where the first logarithm stops.
235    * **Multiplication as ``\\times``.** A centre dot is easy to lose next to a decimal point
236      and next to the slash, at the size these labels are set.
237
238    Precedence is explicit rather than assumed: a term name has exactly one operator at the
239    top level, and any operand that is itself a division is parenthesised on the way in, so
240    ``[1/f] / [g]`` sets as ``(1 / f) / g`` and never as ``1 / f / g``.
241    """
242    text = name.strip()
243    ratio = _split_top(text, " / ")
244    if ratio is not None:
245        numerator, denominator = ratio
246        inner = numerator.strip()
247        if inner.startswith("(") and inner.endswith(")"):
248            summed = _split_top(inner[1:-1], " + ")
249            if summed is not None:
250                left, right = summed
251                summed_math = _math_operand(left) + " + " + _math_operand(right)
252                return "$" + r"\left(" + summed_math + r"\right)" + DIVIDE + _nested(denominator) + "$"
253        return "$" + _nested(numerator) + DIVIDE + _nested(denominator) + "$"
254    product = _split_top(text, " * ")
255    if product is not None:
256        left, right = product
257        return "$" + _nested(left) + r" \times " + _nested(right) + "$"
258    return "$" + _math_operand(text) + "$"
259
260
261def _length_ticks(sizes: np.ndarray, limit: int = 16) -> np.ndarray:
262    """Tick positions for a length axis, thinned so the labels stay readable.
263
264    The curve is reported at every length now, which is right for the detector and wrong for
265    the axis: 32 labels collide into a grey band. Every nth length is labelled instead, with
266    the last one always kept so the axis states its own range.
267    """
268    if sizes.shape[0] <= limit:
269        return sizes
270    step = int(np.ceil(sizes.shape[0] / limit))
271    kept = list(sizes[::step])
272    # Append the final length so the axis states its own range -- unless it would sit on top
273    # of the tick before it, which is what produced a "3132" smudge at the right-hand end.
274    if sizes[-1] not in kept:
275        if kept and sizes[-1] - kept[-1] < step:
276            kept[-1] = sizes[-1]
277        else:
278            kept.append(sizes[-1])
279    return np.asarray(kept)
280
281
282def term_count_curve(
283    curve: pl.DataFrame,
284    destination: str | Path,
285    *,
286    reference: float | None = None,
287    marker: int | None = None,
288    marker_label: str | None = None,
289) -> Path:
290    """Accuracy against equation length: the explainability trade.
291
292    The reference line is the point of the figure. Without it a reader sees a curve still
293    climbing and assumes more terms would keep paying, when the whole approach is bounded
294    well below 1. It is drawn as a labelled line rather than described in text.
295
296    **It is labelled "additive mean-based reference", not "additive ceiling".**
297    `validate.additive_mean_reference` implements the descriptive reference that is additive
298    in *dataset effect plus model effect* -- observed group means and nothing else. E3 carries
299    mixed terms, each multiplying a dataset feature by a model one, so it represents
300    interactions the additive reference cannot and can cross it at longer lengths. Calling
301    that line a ceiling would therefore misstate what it bounds.
302    """
303    figure, axes = plt.subplots(figsize=(7.0, 4.4))
304    sizes = curve["n_terms"].to_numpy()
305
306    axes.plot(sizes, curve["r2_in_sample"].to_numpy(), "o-", color=IN_SAMPLE, label="IS", linewidth=2)
307    if "r2_loo_dataset" in curve.columns:
308        axes.plot(sizes, curve["r2_loo_dataset"].to_numpy(), "s--", color=LOO_DATASET, label="LODO")
309    if "r2_loo_model" in curve.columns:
310        axes.plot(sizes, curve["r2_loo_model"].to_numpy(), "^:", color=LOO_MODEL, label="LOMO")
311    if "r2_loo_cell" in curve.columns:
312        axes.plot(sizes, curve["r2_loo_cell"].to_numpy(), "d-.", color="#D1495B", label="DHO")
313    if reference is not None:
314        axes.axhline(
315            reference,
316            color=CEILING,
317            linestyle="-.",
318            linewidth=1.2,
319            label=f"additive mean-based reference ({reference:.3f})",
320        )
321        # Headroom above whichever is higher. Pinning the top to the reference cropped the
322        # IS curve the moment it crossed -- which is exactly when it matters most.
323        highest = max(float(curve["r2_in_sample"].to_numpy().max()), reference)
324        axes.set_ylim(top=highest + 0.05)
325
326    if marker is not None:
327        # This figure is *about* choosing a length, so the chosen one belongs on it. Named by
328        # the caller for the same reason `error_curve`'s is: a line whose label is fixed in
329        # the plotting code cannot be kept in step with what is actually being passed.
330        name = marker_label or "marked length"
331        axes.axvline(marker, color=CEILING, linestyle=":", linewidth=1.4, label=f"{name} ({marker} terms)")
332
333    axes.set_xlabel("number of terms")
334    axes.set_ylabel("$R^2$")
335    axes.set_xticks(_length_ticks(sizes))
336    axes.grid(alpha=0.25, linestyle=":")
337    axes.legend(loc="lower right", facecolor="white", edgecolor="none", framealpha=0.9)
338    return _finish(figure, destination, paper_width_fraction=0.95)
339
340
341def scatter_limits(
342    truth: np.ndarray,
343    predicted: np.ndarray,
344    margin: float = 0.05,
345    floor: float = 0.0,
346) -> tuple[float, float]:
347    """Square axis bounds covering both series, padded by ``margin`` and stopping at ``floor``.
348
349    Deliberately not MCC's theoretical [-1, 1]. Nothing here approaches -1 -- the observed
350    minimum is -0.29, one row -- and the interval below zero is worse than empty, it is
351    uninformative: around MCC = 0 sit 15 degenerate results and one anti-correlated one,
352    carrying no linear structure for the diagonal to be read against. Including them only
353    compresses the range where the relationship lives.
354    """
355    lowest = max(float(min(truth.min(), predicted.min())) - margin, floor)
356    highest = float(max(truth.max(), predicted.max())) + margin
357    return lowest, highest
358
359
360def count_below_floor(truth: np.ndarray, predicted: np.ndarray, floor: float = 0.0) -> int:
361    """How many points ``predicted_versus_actual`` leaves outside its axes.
362
363    Kept separate from the figure so the disclosure goes in the caption rather than being
364    rendered into the image, where it could not be restyled or edited.
365    """
366    return int(np.sum((truth < floor) | (predicted < floor)))
367
368
369def predicted_versus_actual(
370    truth: np.ndarray,
371    predicted: np.ndarray,
372    destination: str | Path,
373    *,
374    margin: float = 0.05,
375    floor: float = 0.0,
376    groups: np.ndarray | None = None,
377) -> Path:
378    """A scatter against the diagonal.
379
380    17% of the meta-dataset sits at exactly MCC = 1.0 and 8% at exactly 0.0. Both pile-ups
381    are the most important thing to see about this target, and a scatter shows them in a
382    way no summary statistic does.
383
384    ``groups`` is accepted and unused. A rug of the target's marginal distribution was drawn
385    along the bottom axis and removed: on an axis that starts at zero it reads as a row of
386    predictions at MCC = 0, which is exactly the region this target genuinely occupies, and
387    no caption reliably undoes that.
388    """
389    limits = scatter_limits(truth, predicted, margin, floor)
390
391    figure, axes = plt.subplots(figsize=(5.6, 5.4))
392    axes.plot(limits, limits, color=CEILING, linewidth=1.0, linestyle="--", label="perfect", zorder=1)
393    axes.scatter(truth, predicted, s=18, alpha=0.55, color=IN_SAMPLE, edgecolor="none", zorder=2)
394    axes.set_xlim(limits)
395    axes.set_ylim(limits)
396    axes.set_aspect("equal")
397    axes.set_xlabel("actual MCC")
398    axes.set_ylabel("predicted MCC")
399    axes.grid(alpha=0.25, linestyle=":")
400    axes.legend(frameon=False, loc="upper left", fontsize=9)
401    return _finish(figure, destination, paper_width_fraction=0.64)
402
403
404def equation_comparison(comparison: pl.DataFrame, destination: str | Path) -> Path:
405    """Every equation and every ceiling on one scale.
406
407    The headline result of the study. Equations are drawn solid and the reference levels
408    they are read against are drawn hatched immediately beside them, so the gap each
409    equation leaves against its own limit is read directly off the figure instead of
410    computed by the reader from a table.
411
412    **The capability equation is deliberately absent.** It appears in one table in chapter 4
413    and nowhere else, because a bar chart of headline equations is exactly the place a reader
414    would take it for a second recommendation. Its purpose is to answer one question about how
415    far the form reaches, and a figure cannot carry that qualification.
416
417    **"Reference level", not "ceiling", because one of the three is not a ceiling for the
418    bar next to it.** The E1 and E2 references are genuine ceilings -- true per-group means
419    are the most a predictor constant within that group can achieve. The additive mean-based reference is
420    a ceiling only for an equation additive in dataset *and* model effects, and E3 is not
421    one: its mixed terms carry interactions that the additive reference cannot represent. The current
422    E3 passes that reference, which is permitted because a mixed equation is not structurally
423    bounded by it.
424    """
425    table = comparison.filter(~pl.col("equation").str.contains("capability"))
426    labels = table["equation"].to_list()
427    values = table["r2"].to_numpy()
428    is_ceiling = [("reference" in label) or ("oracle" in label) for label in labels]
429
430    figure, axes = plt.subplots(figsize=(8.6, 4.2))
431    positions = np.arange(len(labels))
432    axes.bar(
433        positions,
434        values,
435        color=[CEILING if ceiling else IN_SAMPLE for ceiling in is_ceiling],
436        hatch=["//" if ceiling else "" for ceiling in is_ceiling],
437        alpha=0.85,
438        edgecolor="white",
439    )
440    axes.set_xticks(positions)
441    short_labels = []
442    for label in labels:
443        if label.startswith("E1 reference"):
444            short_labels.append("dataset\nmean")
445        elif label.startswith("E2 reference"):
446            short_labels.append("model\nmean")
447        elif label.startswith("reference: additive"):
448            short_labels.append("additive\nmean")
449        elif label.startswith("E3-Valid"):
450            short_labels.append("E3-\nValid")
451        elif label.startswith("E3-MAX"):
452            short_labels.append("E3-\nMAX")
453        else:
454            short_labels.append(label.split(",", 1)[0])
455    axes.set_xticklabels(short_labels)
456    axes.set_ylabel("$R^2$ on observed pairs")
457    axes.grid(axis="y", alpha=0.25, linestyle=":")
458
459    solid = Rectangle((0, 0), 1, 1, facecolor=IN_SAMPLE, alpha=0.85)
460    hatched = Rectangle((0, 0), 1, 1, facecolor=CEILING, alpha=0.85, hatch="//")
461    axes.legend([solid, hatched], ["fitted equation", "reference level"], frameon=False, fontsize=9)
462    return _finish(figure, destination)
463
464
465def term_effects(effects: pl.DataFrame, destination: str | Path, *, top: int | None = None) -> Path:
466    """Per-term effect sizes in MCC units, signed, strongest at the top.
467
468    ``top`` defaults to every term in the equation so the small-effect terms remain visible;
469    those are the terms a reader most needs to see when judging the brevity argument.
470
471    Terms are rendered as mathematics rather than as their CSV names -- a ratio becomes a
472    fraction, a product a centre dot -- because a figure full of bracketed strings is a
473    wall of punctuation, and the equation is written for a paper. `term_to_math` does it and
474    `TERM_ABBREVIATIONS` shortens the five-syllable model features; chapter 1 expands them.
475    """
476    table = (effects if top is None else effects.head(top)).reverse()
477    labels = [term_to_math(name) for name in table["term"].to_list()]
478    values = table["effect"].to_numpy() * np.sign(table["beta"].to_numpy())
479
480    # At 12pt on the page, each of the 18 term labels needs its own readable row.
481    figure_height = max(2.5, 0.29 * len(labels) + 1.6)
482    figure, axes = plt.subplots(figsize=(8.0, figure_height))
483    axes.barh(
484        range(len(labels)),
485        values,
486        height=0.4,
487        color=[POSITIVE if value > 0 else NEGATIVE for value in values],
488        alpha=0.85,
489    )
490    axes.set_yticks(range(len(labels)))
491    axes.set_yticklabels(labels, fontsize=7.5)
492    axes.axvline(0.0, color="black", linewidth=0.8)
493    axes.set_xlabel("effect on predicted MCC")
494    axes.tick_params(axis="x", labelsize=8)
495    axes.grid(axis="x", alpha=0.25, linestyle=":")
496    return _finish(figure, destination, paper_width_fraction=0.99)
497
498
499def practice_effects(practices: pl.DataFrame, destination: str | Path) -> Path:
500    """What the fitted equation says each raw feature does to MCC, and how stable that is.
501
502    One bar per feature the equation uses. Its length is the change in **predicted MCC**
503    between that feature's lowest and highest decile, holding the rest of the equation --
504    so a bar reaching -0.6 means the equation predicts 0.6 less MCC at the top of that
505    feature's range than at the bottom. Sign is measured on the data rather than read off a
506    weight, because a feature can appear in several terms and inside denominators, and then
507    it has no single weight to read.
508
509    Opacity is the confidence its written practice was rated at, which is driven by how often
510    the feature's terms survived reselection across folds. A long bar at low opacity is a
511    large effect the folds disagreed about; both facts are needed and neither is the other.
512
513    This is the evidence layer of chapter 6, and it is not itself advice: a statement about a
514    meta-feature column becomes a practice only when it supports or contradicts something a
515    practitioner could already have been told.
516    """
517    table = practices.reverse()
518    labels = [TERM_ABBREVIATIONS.get(str(name), str(name)) for name in table["feature"].to_list()]
519    values = table["effect"].to_numpy()
520    alphas = _confidence_alpha(table)
521
522    figure, axes = plt.subplots(figsize=(7.8, 0.44 * len(labels) + 1.4))
523    for index, (value, alpha) in enumerate(zip(values, alphas, strict=True)):
524        axes.barh(index, value, color=POSITIVE if value > 0 else NEGATIVE, alpha=alpha)
525    axes.set_yticks(range(len(labels)))
526    axes.set_yticklabels(labels, fontsize=9)
527    axes.axvline(0.0, color="black", linewidth=0.8)
528    axes.set_xlabel("change in predicted MCC")
529    axes.grid(axis="x", alpha=0.25, linestyle=":")
530
531    # Only the levels actually present, and drawn in the palette the bars use. A legend in a
532    # third colour reads as a third series, and one advertising levels the table does not
533    # contain describes a distinction that is not on the plot.
534    present = _confidence_levels(table)
535    if present:
536        handles = [Rectangle((0, 0), 1, 1, facecolor=POSITIVE, alpha=CONFIDENCE_ALPHA[name]) for name in present]
537        figure.legend(handles, present, frameon=False, loc="upper right", title="confidence")
538    figure.subplots_adjust(left=0.18, right=0.98, bottom=0.24, top=0.72)
539    return _finish(figure, destination, tight_layout=False)
540
541
542#: Bar opacity per rated confidence. ``unrated`` is deliberately the faintest rather than
543#: sharing ``moderate``'s value: an unrated practice is one whose term never stabilised, which
544#: is weaker evidence than a moderate rating, not equal to it.
545CONFIDENCE_ALPHA: dict[str, float] = {"strong": 0.95, "moderate": 0.65, "weak": 0.4, "unrated": 0.25}
546
547
548def _confidence_levels(table: pl.DataFrame) -> list[str]:
549    """The rated levels present in ``table``, strongest first."""
550    if "confidence" not in table.columns:
551        return []
552    seen = {str(value) for value in table["confidence"].to_list()}
553    return [name for name in CONFIDENCE_ALPHA if name in seen]
554
555
556def _confidence_alpha(table: pl.DataFrame) -> list[float]:
557    if "confidence" not in table.columns:
558        return [0.85] * table.height
559    return [CONFIDENCE_ALPHA.get(str(value), 0.25) for value in table["confidence"].to_list()]
560
561
562def _shorten(name: str, limit: int = 46) -> str:
563    return name if len(name) <= limit else name[: limit - 1] + "…"
564
565
566def _wrap(label: str, width: int = 18) -> str:
567    words = label.split()
568    lines: list[str] = []
569    current = ""
570    for word in words:
571        candidate = f"{current} {word}".strip()
572        if len(candidate) > width and current:
573            lines.append(current)
574            current = word
575        else:
576            current = candidate
577    if current:
578        lines.append(current)
579    return "\n".join(lines)
580
581
582def decision_quality(decision: pl.DataFrame, destination: str | Path) -> Path:
583    """The go/no-go decision against threshold, as **F1**, one line per protocol.
584
585    One metric, four lines, and what differs between them is only how much the equation was
586    allowed to see. Previously this drew the equation's accuracy, the equation's F1 and a
587    majority-class baseline on one axis -- two metrics and a predictor, so no reader could
588    tell which metric the baseline was being scored on.
589
590    F1 rather than accuracy: the classes are unbalanced at the outer thresholds, where always
591    answering with the larger one reaches 0.51 accuracy at an F1 of exactly zero. Accuracy
592    would draw four flattering lines and hide that. The trivial predictors stay in the table
593    rather than on the plot, because they cannot be computed under the strictest protocol at
594    all -- a model held out of every fold has no rows to average.
595
596    The four lines are the point: the gap between the top and the bottom one is the whole cost
597    of generalisation on this task, and on this corpus it is small. The legend uses IS for
598    in-sample (IS), leave-one-dataset-out (LODO), leave-one-model-out (LOMO), and doubly
599    held-out (DHO), in which both the dataset and model are absent from training.
600    """
601    figure, axes = plt.subplots(figsize=(6.8, 4.2))
602    table = decision.sort("threshold")
603
604    order = ["in-sample", "loo-dataset", "loo-model", "loo-cell"]
605    styles = {
606        "in-sample": (IN_SAMPLE, "o-", 1.8),
607        "loo-dataset": (LOO_DATASET, "s--", 1.8),
608        "loo-model": (LOO_MODEL, "^:", 1.8),
609        "loo-cell": ("#2f2f2f", "D-", 2.2),
610    }
611    protocol_labels = {
612        "in-sample": "IS",
613        "loo-dataset": "LODO",
614        "loo-model": "LOMO",
615        "loo-cell": "DHO",
616    }
617
618    if "predictor" in table.columns:
619        names = [str(name) for name in table["predictor"].unique(maintain_order=True)]
620        equations = [name for name in names if "equation" in name.lower()]
621        equations.sort(key=lambda name: next((i for i, k in enumerate(order) if k in name), len(order)))
622        for name in equations:
623            rows = table.filter(pl.col("predictor") == name).sort("threshold")
624            key = next((k for k in order if k in name), "in-sample")
625            colour, marker, width = styles[key]
626            axes.plot(
627                rows["threshold"],
628                rows["f1"],
629                marker,
630                color=colour,
631                label=protocol_labels[key],
632                linewidth=width,
633            )
634    else:
635        axes.plot(table["threshold"], table["f1"], "o-", color=IN_SAMPLE, label="equation", linewidth=2)
636
637    axes.set_xlabel("MCC threshold for the go/no-go decision")
638    axes.set_ylabel("F1 of the decision")
639    axes.set_xticks(sorted({float(value) for value in table["threshold"]}))
640    axes.grid(alpha=0.25, linewidth=0.6)
641    axes.legend(fontsize=9, loc="lower left", framealpha=0.0, title="Evaluation protocol")
642    return _finish(figure, destination, paper_width_fraction=0.73)
643
644
645def ranking_quality(selection: pl.DataFrame, destination: str | Path) -> Path:
646    """Per-dataset ranking quality beside what a bad ranking actually *costs*, in MCC.
647
648    Ranking metrics alone raise a question they cannot answer. Four datasets here score
649    average precision below 0.85 and two below 0.5, which reads as serious failure -- until
650    one asks what following the bad ranking costs, and the answer is at most 0.13 MCC and on
651    most of them nothing at all. The two belong together: a dataset whose models all score
652    within a hair of each other is *unrankable*, and ranking it badly is not a cost.
653
654    Two panels sharing the dataset axis rather than two scales on one axis. Regret is in MCC
655    and the ranking metrics are unitless, so overlaying them would repeat the mistake this
656    figure exists to avoid -- a reader cannot tell which axis a marker belongs to. Sorted by
657    average precision, so the hard datasets group at the bottom and the cost of each sits
658    directly beside it.
659    """
660    table = selection.sort("ap", descending=True)
661    positions = np.arange(table.height)
662
663    figure, (axes, cost) = plt.subplots(
664        1,
665        2,
666        figsize=(9.0, 0.37 * table.height + 2.0),
667        sharey=True,
668        gridspec_kw={"width_ratios": [2.1, 1.0], "wspace": 0.10},
669    )
670
671    axes.barh(positions - 0.16, table["ap"], height=0.28, color=IN_SAMPLE, label="average precision")
672    axes.barh(positions + 0.16, table["mrr"], height=0.28, color=LOO_MODEL, label="reciprocal rank")
673    hits = [index for index, value in enumerate(table["hit_at_1"]) if value >= 1.0]
674    if hits:
675        axes.scatter(
676            [0.03] * len(hits),
677            hits,
678            marker="*",
679            s=70,
680            color="white",
681            edgecolor="#333333",
682            linewidth=0.7,
683            zorder=5,
684            label="top choice within 0.01 of best",
685        )
686    axes.set_yticks(positions)
687    axes.set_yticklabels([_shorten(name, 24) for name in table["group"]], fontsize=8)
688    axes.set_xlim(0.0, 1.0)
689    axes.set_xticks([0.0, 0.2, 0.4, 0.6, 0.8])
690    axes.set_xlabel("average precision /\nreciprocal rank")
691    axes.invert_yaxis()
692    axes.grid(alpha=0.25, linewidth=0.6, axis="x")
693    handles, labels = axes.get_legend_handles_labels()
694    figure.legend(handles, labels, loc="upper center", bbox_to_anchor=(0.5, 0.995), ncol=2, frameon=False)
695
696    regret = table["regret"].to_numpy()
697    cost.barh(positions, regret, height=0.42, color=[CEILING if v <= 1e-9 else NEGATIVE for v in regret])
698    cost.set_xlabel("top-choice regret\n(MCC)")
699    cost.grid(alpha=0.25, linewidth=0.6, axis="x")
700    cost.tick_params(labelleft=False)
701    # ``tight_layout`` warns on this shared-y pair in Matplotlib even though the axes are
702    # placed correctly. Set the margins explicitly and retain ``bbox_inches='tight'`` in
703    # `_finish`, which still expands the saved canvas around labels and the legend.
704    figure.subplots_adjust(left=0.30, right=0.98, bottom=0.10, top=0.89, wspace=0.10)
705    return _finish(figure, destination, paper_width_fraction=0.99, tight_layout=False)