ml_meta_perf.report

   1"""Turning a finished study into a written report, without a human in the loop.
   2
   3The whole argument of this project is that the fitted equation can be *read*. That claim
   4would be worth very little if reading it still required an expert to sit down with the
   5weights and narrate what they mean -- the narration would then be the product, and the
   6equation merely its raw material. So every statement this module makes is derived from
   7the fitted object by arithmetic that is written down here, and the same equation always
   8produces the same report.
   9
  10Two ideas carry the analysis.
  11
  12*Standardised weights rank the terms.* The target is centred but never scaled during
  13fitting, so a standardised weight is already in MCC units: it is how far predicted MCC
  14moves when that term moves by one standard deviation of itself. That makes ``beta``
  15directly comparable across terms whose raw units have nothing to do with each other, and
  16it is the number used to decide which terms are major.
  17
  18*Effects state what the terms are worth on this data.* A large weight on a term that
  19barely varies is not important, so ``effect`` -- the swing in predicted MCC across the
  20middle 80% of a term's observed range -- is reported next to every weight. The two
  21disagree often enough to be worth printing together.
  22
  23Feature-level guidance is not re-derived here. It comes from `ml_meta_perf.practices`, which
  24measures direction empirically rather than reading signs off weights, because a feature
  25appearing in two terms or inside a denominator has no single sign to read.
  26
  27Study chapter: [6. Best practices against the equation][study-chapter] -- the rationale, in
  28prose, with the figures.
  29
  30[study-chapter]: https://github.com/mariolpantunes/ml-meta-perf/blob/main/assets/docs/06-practices.md
  31"""
  32
  33from __future__ import annotations
  34
  35from dataclasses import fields
  36from pathlib import Path
  37
  38import numpy as np
  39import polars as pl
  40
  41from ml_meta_perf.attribution import classify, contributions
  42from ml_meta_perf.config import Configuration, SelectionConfig, load_config
  43from ml_meta_perf.data import (
  44    ALL_FEATURES,
  45    DATASET_COLUMN,
  46    DATASET_FEATURES,
  47    FEATURE_GLOSSARY,
  48    MODEL_COLUMN,
  49    MODEL_FEATURES,
  50    corpus_summary,
  51    missing_cells,
  52    target_summary,
  53)
  54from ml_meta_perf.experiment import EquationReport, Report
  55from ml_meta_perf.guidance import as_table as as_guidance_table
  56from ml_meta_perf.guidance import assess, equation_coverage, equation_evidence
  57from ml_meta_perf.guidance import render as render_guidance
  58from ml_meta_perf.model import Equation, direction
  59from ml_meta_perf.practices import render as render_practices
  60from ml_meta_perf.selection import floor_curve, knee_lengths, plateau_start, smoothed
  61from ml_meta_perf.stats import spearman
  62from ml_meta_perf.terms import TRANSFORMS, Atom, Term
  63
  64#: A term is *major* when it falls inside the smallest group of terms whose standardised
  65#: weights account for this share of the total weight mass. Set at four fifths because
  66#: that is where the tail of near-zero weights begins on this data, not because 80% has
  67#: any special status -- ``term_importance`` reports the cumulative share for every term
  68#: so a reader can draw the line somewhere else.
  69MAJOR_MASS = 0.8
  70
  71#: How many terms get written out as sentences. A flat equation can flag most of its
  72#: terms as major, and twenty near-identical paragraphs is not a readable analysis --
  73#: the table carries the rest.
  74SENTENCE_LIMIT = 8
  75
  76#: The grammar's operations, in the order ``operation_usage`` reports them.
  77OPERATIONS = ("atom", "ratio", "product", "sum_ratio", "ratio_of_sums")
  78
  79#: How many raw features each operation names. Used to tell an operation the search
  80#: *declined* from one the arity cap never offered it.
  81OPERATION_ARITY = {"atom": 1, "ratio": 2, "product": 2, "sum_ratio": 3, "ratio_of_sums": 4}
  82
  83
  84def term_importance(
  85    equation: Equation,
  86    columns: dict[str, np.ndarray],
  87    dataset_features: tuple[str, ...],
  88    model_features: tuple[str, ...],
  89    stability: pl.DataFrame | None = None,
  90) -> pl.DataFrame:
  91    """One row per term, ordered by standardised weight, with the major ones flagged.
  92
  93    ``beta`` is the standardised weight: the MCC the term contributes per standard
  94    deviation of itself. ``effect`` is its 10th-to-90th-percentile contribution swing on
  95    the real data. ``share`` is ``|beta|`` as a fraction of the total, and ``cumulative``
  96    its running sum down the ranking, so ``major`` marks the leading terms that together
  97    reach `MAJOR_MASS`.
  98
  99    Ranking on ``beta`` rather than on the raw weight is not cosmetic. Raw weights carry
 100    the units of whatever the term computes, so a term over instance counts and a term
 101    over class entropy have weights that cannot be compared at all; sorting on them would
 102    rank the terms by the size of their units.
 103    """
 104    if not equation.terms:
 105        return pl.DataFrame(
 106            schema={
 107                "rank": pl.Int64,
 108                "term": pl.String,
 109                "group": pl.String,
 110                "features": pl.String,
 111                "weight": pl.Float64,
 112                "beta": pl.Float64,
 113                "effect": pl.Float64,
 114                "share": pl.Float64,
 115                "cumulative": pl.Float64,
 116                "major": pl.Boolean,
 117                "stability": pl.Float64,
 118                "direction": pl.String,
 119            }
 120        )
 121
 122    matrix = contributions(equation, columns)
 123    frequency = (
 124        dict(zip(stability["term"].to_list(), stability["frequency"].to_list(), strict=True))
 125        if stability is not None and stability.height
 126        else {}
 127    )
 128
 129    betas = np.abs(np.asarray(equation.standardized_weights, dtype=np.float64))
 130    total = float(betas.sum())
 131    order = np.argsort(-betas)
 132
 133    rows: list[dict[str, object]] = []
 134    running = 0.0
 135    for position, index in enumerate(order, start=1):
 136        beta = float(equation.standardized_weights[index])
 137        share = float(betas[index] / total) if total > 0.0 else 0.0
 138        running += share
 139        low, high = np.percentile(matrix[:, index], [10.0, 90.0])
 140        term = equation.terms[index]
 141        rows.append(
 142            {
 143                "rank": position,
 144                "term": term.name,
 145                "group": classify(term.features, dataset_features, model_features),
 146                "features": ", ".join(term.features),
 147                "weight": float(equation.weights[index]),
 148                "beta": beta,
 149                "effect": float(high - low),
 150                "share": share,
 151                "cumulative": running,
 152                # The term that crosses the threshold is kept, so the flagged group always
 153                # accounts for at least MAJOR_MASS rather than just short of it.
 154                "major": bool(running - share < MAJOR_MASS),
 155                "stability": float(frequency.get(term.name, float("nan"))),
 156                "direction": direction(beta),
 157            }
 158        )
 159    return pl.DataFrame(rows)
 160
 161
 162def term_sentences(importance: pl.DataFrame, limit: int | None = None) -> list[str]:
 163    """A written line per major term, generated from its row of `term_importance`.
 164
 165    The sentence says three things, in the order a reader needs them: what the term does
 166    to MCC, how much of the equation's weight mass it carries, and how often the folds
 167    agreed it belonged. Nothing is inferred beyond the numbers in the row -- in
 168    particular no claim is made about the *features* inside the term, because a feature's
 169    direction depends on the other terms it appears in and is `ml_meta_perf.practices`' job.
 170    """
 171    if importance.height == 0:
 172        return []
 173    table = importance.filter(pl.col("major")) if "major" in importance.columns else importance
 174    if limit is not None:
 175        table = table.head(limit)
 176
 177    lines: list[str] = []
 178    for row in table.iter_rows(named=True):
 179        beta = float(row["beta"])
 180        stability = float(row["stability"])
 181        agreement = (
 182            "selected in every fold"
 183            if stability >= 0.999
 184            else f"selected in {stability:.0%} of folds"
 185            if not np.isnan(stability)
 186            else "fold agreement not measured"
 187        )
 188        lines.append(
 189            f"`{row['term']}` ({row['group']}) {row['direction']}: one standard deviation of "
 190            f"this term is worth {beta:+.3f} MCC, it moves predicted MCC by {float(row['effect']):.3f} "
 191            f"across the middle 80% of its observed range, it carries {float(row['share']):.1%} of the "
 192            f"equation's weight mass, and it was {agreement}."
 193        )
 194    return lines
 195
 196
 197def group_sentences(shares: pl.DataFrame) -> list[str]:
 198    """A written line per feature group, generated from the share table."""
 199    readings = {
 200        "dataset": "terms over dataset meta-features alone (how hard is this data)",
 201        "model": "terms over model meta-features alone (how capable is this model)",
 202        "mixed": "terms mixing dataset and model features (which model suits which data)",
 203    }
 204    lines: list[str] = []
 205    for row in shares.iter_rows(named=True):
 206        group = str(row["group"])
 207        count = int(row["n_terms"])
 208        lines.append(
 209            f"{readings.get(group, group)}: {count} term{'' if count == 1 else 's'}, "
 210            f"{float(row['share']):.1%} of the equation's output variance."
 211        )
 212    return lines
 213
 214
 215def coverage(importance: pl.DataFrame) -> dict[str, float | int]:
 216    """How concentrated the equation is: how few terms carry how much of it.
 217
 218    ``effective_terms`` is the inverse Simpson index of the weight shares,
 219    $1 / \\sum_i s_i^2$. It answers "how many equally-weighted terms would behave like
 220    this equation": it equals the term count when every term carries the same weight, and
 221    falls toward 1 as one term takes over. It is reported because the headline count and
 222    the major count both depend on where a threshold is drawn, and this does not.
 223    """
 224    if importance.height == 0:
 225        return {
 226            "n_terms": 0,
 227            "n_major": 0,
 228            "major_share": 0.0,
 229            "top_share": 0.0,
 230            "effective_terms": 0.0,
 231        }
 232    major = importance.filter(pl.col("major"))
 233    shares = importance["share"].to_numpy()
 234    return {
 235        "n_terms": importance.height,
 236        "n_major": major.height,
 237        "major_share": float(major["share"].sum()),
 238        "top_share": float(importance["share"][0]),
 239        "effective_terms": float(1.0 / np.sum(shares**2)),
 240    }
 241
 242
 243def unstable_majors(importance: pl.DataFrame, threshold: float = 0.5) -> pl.DataFrame:
 244    """Major terms the folds did not agree on, which are the ones not to build advice on.
 245
 246    A large standardised weight and a low selection frequency together mean the term is
 247    doing a lot of work for *this* training set and would be replaced by something else
 248    on another. That combination is invisible in the equation as printed -- it looks like
 249    any other large coefficient -- so it is extracted rather than left to be noticed.
 250    """
 251    if importance.height == 0 or "stability" not in importance.columns:
 252        return importance
 253    return importance.filter(
 254        pl.col("major") & pl.col("stability").is_not_nan() & (pl.col("stability") < threshold)
 255    ).select("rank", "term", "group", "beta", "share", "stability")
 256
 257
 258def _atoms(term: Term | Atom) -> list[Atom]:
 259    """Every `Atom` inside a term, however deeply nested."""
 260    if isinstance(term, Atom):
 261        return [term]
 262    collected: list[Atom] = []
 263    for operand in term.operands:
 264        collected.extend(_atoms(operand))
 265    return collected
 266
 267
 268def feature_usage(
 269    equation: Equation,
 270    importance: pl.DataFrame,
 271    available: tuple[str, ...],
 272) -> pl.DataFrame:
 273    """Which raw features the equation reached for, how often, and under what transforms.
 274
 275    A second way to read a flat equation. When no single term dominates, the question
 276    "which term matters" has no useful answer, but "which of the features on offer did the
 277    search actually use, and how did it have to bend them" still does -- and it is asked of
 278    the vocabulary rather than of the weights, so a spread of weights does not blunt it.
 279
 280    ``share`` sums the standardised-weight mass of every term a feature appears in. A
 281    feature in two terms is credited both, so the column does not sum to 1: it answers
 282    "how much of the equation touches this feature", not "how much does this feature own".
 283    ``unused`` features are reported too, since a feature the search declined to use after
 284    seeing every transform of it is a finding about the meta-data.
 285    """
 286    mass = dict(zip(importance["term"].to_list(), importance["share"].to_list(), strict=True))
 287    rows: list[dict[str, object]] = []
 288    for feature in available:
 289        indices = [index for index, term in enumerate(equation.terms) if feature in term.features]
 290        transforms = sorted(
 291            {atom.transform for index in indices for atom in _atoms(equation.terms[index]) if atom.feature == feature}
 292        )
 293        operations = sorted({equation.terms[index].operation for index in indices})
 294        rows.append(
 295            {
 296                "feature": feature,
 297                "meaning": FEATURE_GLOSSARY.get(feature, feature),
 298                "n_terms": len(indices),
 299                "share": float(sum(mass.get(equation.terms[index].name, 0.0) for index in indices)),
 300                "transforms": ", ".join(transforms),
 301                "operations": ", ".join(operations),
 302            }
 303        )
 304    return pl.DataFrame(rows).sort("share", descending=True)
 305
 306
 307def operation_usage(
 308    equation: Equation,
 309    importance: pl.DataFrame,
 310    max_arity: int | None = None,
 311) -> pl.DataFrame:
 312    """Which grammar operations and transforms earned their place in the equation.
 313
 314    The vocabulary offers five operations and five transforms; the search is free to
 315    ignore any of them, and what it declined is as much a result as what it chose --
 316    every unused entry is a shape the data turned out not to need.
 317
 318    That reading only holds for entries the search could actually have used, so
 319    ``max_arity`` marks the rest ``offered = no``. Reporting ``ratio_of_sums`` as
 320    "declined" when a three-feature cap kept it out of the library would be reading a
 321    configuration choice as a finding about the data.
 322    """
 323    mass = dict(zip(importance["term"].to_list(), importance["share"].to_list(), strict=True))
 324    rows: list[dict[str, object]] = []
 325    for kind, names in (("operation", OPERATIONS), ("transform", TRANSFORMS)):
 326        for name in names:
 327            if kind == "operation":
 328                indices = [index for index, term in enumerate(equation.terms) if term.operation == name]
 329                offered = max_arity is None or OPERATION_ARITY[name] <= max_arity
 330            else:
 331                indices = [
 332                    index
 333                    for index, term in enumerate(equation.terms)
 334                    if any(atom.transform == name for atom in _atoms(term))
 335                ]
 336                offered = True
 337            rows.append(
 338                {
 339                    "kind": kind,
 340                    "name": name,
 341                    "offered": offered,
 342                    "n_terms": len(indices),
 343                    "share": float(sum(mass.get(equation.terms[index].name, 0.0) for index in indices)),
 344                }
 345            )
 346    return pl.DataFrame(rows)
 347
 348
 349def _shared_features(equation: Equation, members: list[int]) -> list[str]:
 350    """Features appearing in at least half of a block's terms, commonest first.
 351
 352    What names a block. The terms were grouped on how their contributions move, not on
 353    what they contain, so this is a genuine question rather than a restatement of the
 354    grouping -- and when a block has no common feature, the empty answer is worth seeing
 355    too.
 356    """
 357    counts: dict[str, int] = {}
 358    for index in members:
 359        for feature in set(equation.terms[index].features):
 360            counts[feature] = counts.get(feature, 0) + 1
 361    # Strict majority. Half would let a two-term block claim a feature only one of its
 362    # terms names, which is the opposite of shared.
 363    threshold = len(members) // 2 + 1
 364    return [
 365        feature for feature, count in sorted(counts.items(), key=lambda item: (-item[1], item[0])) if count >= threshold
 366    ]
 367
 368
 369def term_groups(
 370    equation: Equation,
 371    columns: dict[str, np.ndarray],
 372    importance: pl.DataFrame,
 373    *,
 374    threshold: float = 0.5,
 375) -> pl.DataFrame:
 376    """Terms whose contributions move together, read as one block.
 377
 378    An additive equation invites reading term by term, and that works when one or two
 379    weights dominate. When they do not -- and here they do not -- the honest unit of
 380    explanation is larger than a term and smaller than the whole equation. Terms are
 381    grouped by the correlation of their per-row contributions, so a group is a set of terms
 382    that rise and fall together across the meta-dataset and can be described in one
 383    sentence.
 384
 385    Grouping on contributions rather than on shared features is deliberate. Two terms can
 386    share no feature and still track each other, and two terms over the same feature can
 387    move independently once their transforms differ; what a reader needs to know is whether
 388    they say the same thing about a row.
 389
 390    Correlation is signed and taken on the contribution, so a term entering with a negative
 391    weight joins the group it actually agrees with rather than the one it merely shares
 392    features with.
 393    """
 394    if not equation.terms:
 395        return pl.DataFrame(
 396            schema={
 397                "group": pl.Int64,
 398                "n_terms": pl.Int64,
 399                "share": pl.Float64,
 400                "effect": pl.Float64,
 401                "direction": pl.String,
 402                "shared": pl.String,
 403                "terms": pl.String,
 404            }
 405        )
 406
 407    matrix = contributions(equation, columns)
 408    centered = matrix - matrix.mean(axis=0)
 409    scale = np.sqrt((centered**2).sum(axis=0))
 410    scale = np.where(scale < 1e-12, 1.0, scale)
 411    correlation = (centered / scale).T @ (centered / scale)
 412
 413    mass = dict(zip(importance["term"].to_list(), importance["share"].to_list(), strict=True))
 414    # Single-link agglomeration: a term joins a group when it tracks *any* member, which is
 415    # the right rule here because a chain of terms that each track the next is one story
 416    # told in several pieces, not several stories.
 417    assignment = list(range(len(equation.terms)))
 418    for i in range(len(equation.terms)):
 419        for j in range(i + 1, len(equation.terms)):
 420            if correlation[i, j] >= threshold:
 421                merged, absorbed = sorted((assignment[i], assignment[j]))
 422                assignment = [merged if label == absorbed else label for label in assignment]
 423
 424    rows: list[dict[str, object]] = []
 425    for label in sorted(set(assignment)):
 426        members = [index for index, value in enumerate(assignment) if value == label]
 427        total = matrix[:, members].sum(axis=1)
 428        low, high = np.percentile(total, [10.0, 90.0])
 429        swing = float(high - low)
 430        # The group's own sign: whether its combined contribution rises with itself, which
 431        # is what a one-sentence description of the block has to state.
 432        signed = float(np.sum([equation.standardized_weights[index] for index in members]))
 433        rows.append(
 434            {
 435                "group": len(rows) + 1,
 436                "n_terms": len(members),
 437                "share": float(sum(mass.get(equation.terms[index].name, 0.0) for index in members)),
 438                "effect": swing,
 439                "direction": direction(signed),
 440                "shared": ", ".join(_shared_features(equation, members)),
 441                "terms": " ; ".join(equation.terms[index].name for index in members),
 442            }
 443        )
 444    return pl.DataFrame(rows).sort("share", descending=True).with_columns(pl.int_range(1, pl.len() + 1).alias("group"))
 445
 446
 447def marginal_versus_conditional(
 448    practices: pl.DataFrame,
 449    columns: dict[str, np.ndarray],
 450    truth: np.ndarray,
 451) -> pl.DataFrame:
 452    """Each practice's direction against the feature's own correlation with MCC.
 453
 454    A practice states what the *equation* does as a feature rises, with every other term
 455    present. A marginal correlation states what the feature does alone. The two are
 456    different quantities and on this data they mostly disagree, which is what conditioning
 457    does rather than a defect: a marginal correlation mixes a feature's effect with
 458    everything it travels with, while inside the equation the terms carrying those
 459    companions are already there.
 460
 461    It is worth printing because the marginal view is the one a reader will accidentally
 462    check against -- plot the feature against MCC and be puzzled -- so the disagreement is
 463    better stated than discovered.
 464    """
 465    if practices.height == 0 or "direction" not in practices.columns:
 466        return pl.DataFrame(
 467            schema={
 468                "feature": pl.String,
 469                "meaning": pl.String,
 470                "marginal": pl.Float64,
 471                "conditional": pl.Float64,
 472                "practice_says": pl.String,
 473                "agrees": pl.Boolean,
 474            }
 475        )
 476    rows: list[dict[str, object]] = []
 477    for row in practices.iter_rows(named=True):
 478        feature = str(row["feature"])
 479        marginal = spearman(columns[feature], truth)
 480        conditional = float(row["direction"])
 481        rows.append(
 482            {
 483                "feature": feature,
 484                "meaning": str(row.get("meaning", feature)),
 485                "marginal": marginal,
 486                "conditional": conditional,
 487                "practice_says": "higher" if conditional > 0.0 else "lower",
 488                "agrees": bool(np.sign(marginal) == np.sign(conditional)),
 489            }
 490        )
 491    return pl.DataFrame(rows)
 492
 493
 494_PROTOCOL_LABELS = (
 495    ("loo_cell", "DHO"),
 496    ("loo_dataset", "LODO"),
 497    ("loo_model", "LOMO"),
 498    ("in_sample", "IS"),
 499    ("loo-cell: both held out", "DHO"),
 500    ("loo-dataset", "LODO"),
 501    ("loo-model", "LOMO"),
 502    ("doubly-held-out", "DHO"),
 503    ("doubly held out", "DHO"),
 504    ("leave-one-dataset-out", "LODO"),
 505    ("leave-one-model-out", "LOMO"),
 506    ("in-sample", "IS"),
 507)
 508
 509
 510def _published_protocol_label(value: object) -> str:
 511    """Translate stable internal protocol names at the publication boundary."""
 512    text = str(value)
 513    for internal, published in _PROTOCOL_LABELS:
 514        text = text.replace(internal, published)
 515    return text
 516
 517
 518def _table(frame: pl.DataFrame, float_format: str = "{:.4f}") -> str:
 519    """A polars frame as a GitHub-flavoured markdown table.
 520
 521    Result files keep the stable internal protocol keys. Published tables use the same
 522    concise labels as the figures: IS, LODO, LOMO, and DHO.
 523    """
 524    if frame.height == 0:
 525        return "_(empty)_"
 526
 527    def cell(value: object) -> str:
 528        # "no" rather than an empty cell: a blank reads as a missing measurement, which
 529        # is a different claim from a measured disagreement.
 530        if isinstance(value, bool):
 531            return "yes" if value else "no"
 532        if value is None:
 533            return ""
 534        if isinstance(value, float):
 535            return "" if np.isnan(value) else float_format.format(value)
 536        return _published_protocol_label(value).replace("|", r"\|")
 537
 538    header = "| " + " | ".join(_published_protocol_label(column) for column in frame.columns) + " |"
 539    rule = "|" + "|".join("---" for _ in frame.columns) + "|"
 540    body = ["| " + " | ".join(cell(value) for value in row) + " |" for row in frame.iter_rows()]
 541    return "\n".join([header, rule, *body])
 542
 543
 544def single_prediction(
 545    equation: Equation,
 546    columns: dict[str, np.ndarray],
 547    truth: np.ndarray,
 548    frame: pl.DataFrame,
 549    *,
 550    limit: int = 6,
 551) -> str:
 552    """One row's prediction broken into its terms, as a printed block.
 553
 554    The interpretability payoff in its most direct form: a reader can add the column up and
 555    get the prediction back, which is a thing no opaque regressor lets anyone do.
 556
 557    **The row is chosen by a rule, not picked.** The equation's median absolute error is a
 558    typical prediction rather than its best case, and choosing it arithmetically is what keeps
 559    this example from being a flattering one somebody selected once. The chapter this lands in
 560    is the chapter arguing that the analysis is generated rather than authored; it carried a
 561    hand-written breakdown until 2026-09-07, and every figure in it had gone stale -- the
 562    intercept it printed was 1.4787 against the equation's actual 1.3586.
 563
 564    Ties are broken by taking the first such row in corpus order, so the choice is a function
 565    of the data and the equation alone.
 566    """
 567    predicted = equation.predict(columns)
 568    error = np.abs(truth - predicted)
 569    # The median error itself, not the row at the middle of the sorted order: with an even
 570    # row count those differ, and `argmin` against the median value picks an actual row.
 571    row = int(np.argmin(np.abs(error - float(np.median(error)))))
 572
 573    shares = contributions(equation, columns)[row]
 574    ordered = sorted(zip(shares, equation.terms, strict=True), key=lambda pair: -abs(pair[0]))
 575    total = equation.intercept + float(shares.sum())
 576
 577    lines = [
 578        f"dataset  : {frame[DATASET_COLUMN][row]}",
 579        f"model    : {frame[MODEL_COLUMN][row]}",
 580        f"actual   : {truth[row]:+.4f}",
 581        f"predicted: {predicted[row]:+.4f}",
 582        "",
 583        "contribution breakdown:",
 584        f"    {equation.intercept:+.4f}   intercept",
 585    ]
 586    for value, term in ordered[:limit]:
 587        lines.append(f"    {value:+.4f}   {term.name}")
 588    if len(ordered) > limit:
 589        remainder = sum(value for value, _ in ordered[limit:])
 590        lines.append(f"    {remainder:+.4f}   the remaining {len(ordered) - limit} terms")
 591    lines.append(f"  = {total:+.4f}   sum")
 592    # `Equation.predict` clips to the range MCC can take at all, so on a row where the linear
 593    # form runs past the end of that range the column does not add up to the number printed
 594    # above it. Saying so is the honest form: a breakdown a reader cannot add up, with no
 595    # explanation of why, is worse than no breakdown.
 596    if abs(total - predicted[row]) > 5e-5:
 597        lines.append(f"    clipped to {predicted[row]:+.4f}, the range MCC can take")
 598    return "```\n" + "\n".join(lines) + "\n```"
 599
 600
 601#: Only E1 is constant within its identity group. E2's model-mean reference is descriptive
 602#: because Processing Units Number can vary within a model, so no ceiling ratio is shown.
 603CEILING_ROWS = {"E1": "E1 reference: true dataset means"}
 604
 605
 606def _headline(report: Report) -> pl.DataFrame:
 607    """The equations, with a ceiling ratio only for dataset-constant E1.
 608
 609    The one table the study can be summarised by. It exists because the summary page kept a
 610    hand-written copy of these numbers and had no way of noticing when they moved -- the same
 611    failure the generated chapter sections were introduced to remove, reappearing one level
 612    up. **The ``reached`` column is the point of the table**, not the R2 column: E1's
 613    structural maximum is a per-dataset constant, so its score and E3's score are not
 614    comparable as achievements, and printing the fraction of each equation's own ceiling
 615    beside them is what stops a reader making that comparison anyway.
 616    """
 617    ceilings = dict(zip(report.comparison["equation"].to_list(), report.comparison["r2"].to_list(), strict=True))
 618    rows: list[dict[str, object]] = []
 619    listed = [("E1", report.e1), ("E2", report.e2), ("E3-Valid", report.e3)]
 620    if report.e3_capability is not report.e3:
 621        listed.append(("E3-MAX", report.e3_capability))
 622    for name, equation in listed:
 623        ceiling = ceilings.get(CEILING_ROWS.get(name, ""))
 624        in_sample = float(equation.in_sample["r2"])
 625        rows.append(
 626            {
 627                "equation": name,
 628                "features": {"E1": "dataset", "E2": "model"}.get(name, "both"),
 629                "terms": equation.equation.n_terms,
 630                "IS R2": in_sample,
 631                "LODO R2": float(equation.cross_validated["loo_dataset"]["r2"]),
 632                "LOMO R2": float(equation.cross_validated["loo_model"]["r2"]),
 633                "DHO R2": float(equation.cross_validated["loo_cell"]["r2"]),
 634                "own ceiling": float(ceiling) if ceiling is not None else float("nan"),
 635                "reached": in_sample / float(ceiling) if ceiling else float("nan"),
 636            }
 637        )
 638    return pl.DataFrame(rows)
 639
 640
 641def _scores(label: str, scores: dict[str, float | int]) -> str:
 642    return (
 643        f"| {label} | {float(scores['r2']):.4f} | {float(scores['mae']):.4f} | "
 644        f"{float(scores['rmse']):.4f} | {int(scores['n'])} |"
 645    )
 646
 647
 648def _crater_note(report: Report, selection: SelectionConfig) -> str:
 649    """The deepest single-length crater on the worst-protocol curve, and what smoothing does to it.
 650
 651    The length rule reads the floor through a running median, and the argument for that is
 652    only convincing next to a length where one fold dragged the curve down. Which length that
 653    is moves whenever the configuration does, so it is found here rather than named in prose.
 654
 655    Deepest is measured against the neighbouring lengths, not against the curve's own mean: a
 656    crater is a local collapse, and a length at the end of a declining run is not one.
 657    """
 658    curve = report.e3.curve
 659    if curve.height < 3:
 660        return ""
 661    floor = floor_curve(curve)
 662    lengths = curve["n_terms"].to_numpy()
 663    neighbours = (floor[:-2] + floor[2:]) / 2.0
 664    index = int(np.argmax(neighbours - floor[1:-1])) + 1
 665    depth = float(neighbours[index - 1] - floor[index])
 666    if depth <= 0.0:
 667        return ""
 668    smooth = smoothed(floor, selection.smoothing)
 669    return (
 670        f"**Why the curve is smoothed before a length is chosen.** The deepest crater on this "
 671        f"curve is at **{int(lengths[index])} terms**, where the worst protocol reads "
 672        f"{float(floor[index]):.3f}, {depth:.3f} below the neighbouring lengths. The running median "
 673        f"of {selection.smoothing} lengths reads {float(smooth[index]):.3f} there. A crater like this "
 674        "is not a property of the length -- it is one held-out dataset sitting outside the convex "
 675        "hull of the other nineteen in term space, where a linear equation extrapolates without "
 676        "limit and `validate._clip_to_training` pins the fold to its training floor -- and the "
 677        "same holds for a single lucky length. Neither should choose the published length.\n"
 678    )
 679
 680
 681def _shared_terms_note(report: Report) -> str:
 682    """How many of each control's terms E3-Valid uses exactly, and how many of E3's are mixed."""
 683    published = {term.name for term in report.e3.equation.terms}
 684    rows = []
 685    for label, control in (("E2 → E3", report.e2), ("E1 → E3", report.e1)):
 686        names = [term.name for term in control.equation.terms]
 687        shared = [name for name in names if name in published]
 688        rows.append(
 689            f"| {label} | **{len(shared)} of {len(names)}** | "
 690            + ("; ".join(f"`{name}`" for name in shared) or "—")
 691            + " |"
 692        )
 693    mixed = sum(
 694        1 for term in report.e3.equation.terms if classify(term.features, DATASET_FEATURES, MODEL_FEATURES) == "mixed"
 695    )
 696    return (
 697        "How much of each control survives into E3-Valid, counting terms shared *exactly*:\n\n"
 698        "| | shared with E3 | which |\n|---|---|---|\n"
 699        + "\n".join(rows)
 700        + f"\n\n{mixed} of E3-Valid's {report.e3.equation.n_terms} terms mix dataset and model "
 701        "features, and a mixed term is available to neither control by construction.\n"
 702    )
 703
 704
 705def _optimism_note(report: Report) -> str:
 706    """The fixed-form LODO against LODO with the form re-chosen per fold, and the gap."""
 707    if report.optimism.height < 2:
 708        return ""
 709    reported, nested = (float(value) for value in report.optimism["r2"].to_list()[:2])
 710    return (
 711        "Every reported held-out number fixes the equation's form -- chosen once, on all observed "
 712        "rows -- and refits only its weights in each fold. The form was therefore chosen with "
 713        "the held-out dataset in view. The second row repeats the *selection* inside every "
 714        f"LODO fold, at E3-Valid's {report.e3.n_terms} terms and configuration, and scores those "
 715        "predictions instead:\n\n"
 716        + _table(report.optimism.select("LODO", "r2", "mae", "spearman"))
 717        + f"\n\n**LODO R² differs by {reported - nested:.3f}** between the globally selected form "
 718        f"({reported:.3f}) and within-fold term reselection ({nested:.3f}). This is a partial "
 719        "selection-sensitivity diagnostic: global configuration, requested length, and library "
 720        "construction remain shared. Neither row independently assesses the complete discovery "
 721        "procedure. The fixed-form scores are conditional on selection.\n"
 722    )
 723
 724
 725def _saturated_note(report: Report) -> str:
 726    """Read the saturated fit against the published equation, in one generated sentence.
 727
 728    The comparison is the point, and stating it in prose beside the table is how it stops
 729    being two numbers a reader has to subtract. Both sides come from this run, so the
 730    sentence cannot drift from the table above it.
 731    """
 732    if report.saturated.height == 0:
 733        return ""
 734    row = report.saturated.to_dicts()[0]
 735    published = float(report.e3.cross_validated["loo_dataset"]["r2"])
 736    return (
 737        f"**The solver is not the hard part; the sample size is.** All "
 738        f"{int(row['terms'])} terms at once fit better under IS than the published equation "
 739        f"({float(row['r2_in_sample']):.4f} against "
 740        f"{float(report.e3.in_sample['r2']):.4f}) and transfer at "
 741        f"{float(row['r2_loo_dataset_clipped']):.4f} under LODO, against the "
 742        f"published equation's {published:.4f}. The unclipped figure — "
 743        f"{float(row['r2_loo_dataset_unclipped']):.1f} — is what the fit does when a held-out "
 744        "dataset falls outside the convex hull of the other nineteen and nothing bounds the "
 745        "extrapolation. A design this much wider than 20 held-out groups can support has "
 746        "nothing to constrain it, which is what selection is for.\n"
 747    )
 748
 749
 750def _coverage_note(report: Report, columns: dict[str, np.ndarray]) -> str:
 751    """How far the equation reaches into the catalogue, and how to read a disagreement.
 752
 753    Two counts a reader will otherwise conflate: the verdict tally above is drawn from corpus
 754    statistics that any study with this data could compute, while this is about the terms. And
 755    one reading that has to be spelled out, because the table looks self-contradictory without
 756    it -- a feature in a denominator enters inverted, so the same feature can carry opposite
 757    signs in two terms and mean one thing.
 758    """
 759    counts = equation_coverage(report, columns)
 760    if not counts["pairs"]:
 761        return ""
 762    lines = [
 763        f"**Read this against the verdict tally above, not as part of it.** Of the "
 764        f"{counts['practices']} practices, {counts['with_feature_claims']} make a claim about a "
 765        "quantity the equation contains; the rest are about a protocol, a metric, or a family "
 766        "of learners, and pairing one of those with a coefficient would be inventing a "
 767        f"connection. Those are carried by **{counts['carrying_terms']} of the equation's "
 768        f"{counts['terms']} terms**, giving {counts['pairs']} directed (practice, term) "
 769        f"pairings: **{counts['agree']} come out the way the practice predicts and "
 770        f"{counts['disagree']} do not**"
 771        + (
 772            f", with {counts['undirected']} pairing too weak to state a direction for and "
 773            f"{counts['unselected']} claim resting on a feature the search never took.\n"
 774            if counts["undirected"] or counts["unselected"]
 775            else ".\n"
 776        ),
 777    ]
 778    if counts["disagree"] and counts["disagree_in_denominator"] == counts["disagree"]:
 779        lines.append(
 780            f"\n**Every one of the {counts['disagree']} disagreements is the same feature "
 781            "entering as a *denominator*, and that is arithmetic rather than conflict.** A "
 782            "negatively-weighted ratio contributes more as its denominator grows, so a term "
 783            "of the form `dataset property / capacity` must rise with capacity. Read the "
 784            "`position` column across those rows and the equation is saying one coherent "
 785            "thing: it has **no marginal statement about capacity at all**, only statements "
 786            "about capacity *relative to* something the dataset demands. Where capacity is a "
 787            "numerator instead — measured against the class count rather than against a "
 788            "difficulty — it lowers predicted MCC, which is the practice's own claim.\n"
 789            "\n**So the mismatch is in how the expectation was written down, not in what the "
 790            "equation says.** A practice recommending that capacity be *matched to the "
 791            "problem* is a conditional claim; encoding it as `capacity lowers MCC` is a "
 792            "marginal one, and a marginal expectation cannot match a term that only ever "
 793            "speaks conditionally. This is the clearest case in the study of why the equation "
 794            "is worth reading term by term: no per-feature summary, and no opaque model, can "
 795            "distinguish 'more capacity is better' from 'more capacity per unit of difficulty "
 796            "is better'.\n"
 797        )
 798    else:
 799        lines.append(
 800            "\n**Read the `position` column before calling a row a disagreement.** A feature "
 801            "in a denominator enters the term inverted, so a negatively-weighted ratio "
 802            "contributes more as that feature rises. The same feature carrying opposite signs "
 803            "in two terms is usually one statement about a ratio rather than two conflicting "
 804            "ones about a quantity.\n"
 805        )
 806    lines.append(
 807        "\n``beta`` is the strength and ``effect`` is what the term is worth on this data; "
 808        "**agreement in sign with a negligible effect is agreement without evidence**, which "
 809        "is why the two are printed together. ``rho`` is the measured association the "
 810        "direction comes from, and a pairing below the same floor the per-feature statements "
 811        "use is reported as having no direction rather than being given a sign it cannot "
 812        "support.\n"
 813    )
 814    return "\n".join(lines)
 815
 816
 817def _identity_note(report: Report) -> str:
 818    """What the identity ceiling is worth, and whether the study's own test calls it real.
 819
 820    A difference of two pooled R2 values over twenty folds is not a measurement -- that has
 821    produced three wrong conclusions on this project -- so the rungs are paired against the
 822    uncorrected equation and the verdict comes from the test rather than from the size of the
 823    gap. The two halves of the paired test are reported separately where they disagree,
 824    because they answer different questions: the sign test asks whether the correction wins
 825    *consistently*, the bootstrap whether the mean gain survives a different draw of datasets.
 826    """
 827    if report.identity.height < 3:
 828        return ""
 829    rows = report.identity.to_dicts()
 830    base, level, slope = rows[0], rows[1], rows[2]
 831    total = float(slope["r2_loo_dataset"]) - float(base["r2_loo_dataset"])
 832    step = float(level["r2_loo_dataset"]) - float(base["r2_loo_dataset"])
 833
 834    lines = [
 835        f"**The gap is {total:.3f} of LODO R²**, of which a per-model *level* "
 836        f"recovers {step:.3f} and the level-plus-slope form the remaining {total - step:.3f}.\n",
 837        "**Whether that is real is a paired question**, so each rung is compared with the "
 838        "uncorrected equation dataset by dataset, on absolute error, over the twenty held-out "
 839        "folds. The two rungs come back differently, and the difference is the finding:\n",
 840    ]
 841    if level["verdict"] == "tie":
 842        lines.append(
 843            f"* a per-model **level** is a **tie** -- it wins on {int(level['wins'])} of the "
 844            "twenty folds and its interval spans zero. A constant shift per model, which is "
 845            "what a level is, adds nothing the equation does not already have.\n"
 846        )
 847    if slope["verdict"] == "real":
 848        lines.append(
 849            f"* a per-model **slope** is **not** a tie: the bootstrap interval "
 850            f"[{float(slope['ci_low']):+.4f}, {float(slope['ci_high']):+.4f}] lies entirely "
 851            f"above zero, on {int(slope['wins'])} winning folds of twenty. Read it with the "
 852            f"sign test beside it, which at p = {float(slope['sign_p']):.3f} does **not** reach "
 853            "significance -- so the gain is carried by its size on the folds it wins rather "
 854            "than by winning nearly all of them. That is weaker evidence than the interval "
 855            "alone suggests, and stronger than a tie.\n"
 856        )
 857        lines.append(
 858            "**So the question this chapter was written to close is not closed.** The half of "
 859            "the correction that survives is the one that lets a model's advantage depend on "
 860            "the data -- exactly what a *capability* descriptor would have to do, and exactly "
 861            "what none of the descriptors this corpus records does. The mixed terms were "
 862            "supposed to absorb that interaction and have absorbed only part of it.\n"
 863        )
 864    else:
 865        lines.append(
 866            "* neither rung survives pairing, so on this corpus perfect model identity adds "
 867            "nothing measurable to the published equation and the model side is as well "
 868            "described as free per-model numbers could make it.\n"
 869        )
 870    lines.append(
 871        "This correction tests information carried by model histories. Its outcome does not "
 872        "establish which alternative descriptors or categorical encodings would transfer to "
 873        "unseen learners. Encoding comparisons remain conditional on form selection.\n"
 874    )
 875    return "\n".join(lines)
 876
 877
 878def _opaque_note(report: Report) -> str:
 879    """Read the opaque table across its protocols, and against the equation on each.
 880
 881    The columns are ordered by how much the estimator was allowed to see, and reading them
 882    left to right is the argument. The last one is the comparison that matters: under
 883    Under DHO neither side has the dataset or the model, so it is the only protocol on
 884    which an opaque regressor and the published equation are denied the same things.
 885    """
 886    if report.opaque.height == 0 or "r2_loo_cell" not in report.opaque.columns:
 887        return ""
 888    rows = report.opaque.to_dicts()
 889    best_fit = max(rows, key=lambda row: float(row["r2_in_sample"]))
 890    best_cell = max(rows, key=lambda row: float(row["r2_loo_cell"]))
 891    return (
 892        f"**Read the {best_fit['model']} row across.** It fits this meta-data at R² "
 893        f"{float(best_fit['r2_in_sample']):.4f}; holding out a whole model leaves it at "
 894        f"{float(best_fit['r2_loo_model']):.4f}; holding out a whole dataset drops it to "
 895        f"{float(best_fit['r2_loo_dataset']):.4f}; and with **both** held out it reaches "
 896        f"{float(best_fit['r2_loo_cell']):.4f}. The published equation is at "
 897        f"{float(report.e3.in_sample['r2']):.4f} and "
 898        f"{float(report.e3.cross_validated['loo_dataset']['r2']):.4f} in the corresponding "
 899        "IS and LODO settings.\n"
 900        "\nThese are fixed-setting comparisons across grouped coefficient-fit partitions. "
 901        "Lower scores after dataset exclusion show sensitivity to the available task histories; "
 902        "they do not by themselves identify the mechanism of a regressor's errors.\n"
 903        f"\n**Under joint exclusion the best reference R² is "
 904        f"{float(best_cell['r2_loo_cell']):.4f}** ({best_cell['model']}). Per-model historical "
 905        "baselines are unavailable under DHO because the target learner has no fitting rows. "
 906        "Regression, ranking, and threshold outcomes must be read separately.\n"
 907        "\nRidgeCV selects its penalty within fitting partitions. The tree ensembles retain "
 908        "documented study-development settings without fresh tuning on the revised corpus. "
 909        "Their raw features and prior selection budgets differ from the globally selected "
 910        "symbolic form; these comparisons do not independently establish superiority of "
 911        "the complete discovery procedure.\n"
 912    )
 913
 914
 915def _length_note(report: Report, selection: SelectionConfig) -> str:
 916    """Explain the length rule with this run's knees, and the paired sensitivity table."""
 917    if report.length_choice.height == 0:
 918        return ""
 919    knees = knee_lengths(report.e3.curve, selection.smoothing)
 920    listed = ", ".join(str(knee) for knee in knees) if knees else "none"
 921    _, how = plateau_start(
 922        report.e3.curve, delta=selection.delta, window=selection.window, smoothing=selection.smoothing
 923    )
 924    reached = (
 925        "It is a knee: the curve stops improving right after it."
 926        if how == "knee"
 927        else "None of the knees is followed by a plateau -- the curve keeps gaining past each -- so "
 928        "the plateau starts at the smoothed maximum."
 929    )
 930    lines = [
 931        f"E3-Valid has **{report.e3.n_terms} terms**. The rule, `selection.plateau_knee`, is "
 932        "the same for E1, E2 and E3-Valid. It reads the worst R² over IS, LODO, LOMO and DHO at "
 933        f"every length, smoothed by a running median of {selection.smoothing} lengths; "
 934        "multi-Kneedle (`kneeliverse`) proposes the lengths where that curve's Pareto front "
 935        f"bends -- here **{listed}** -- and the rule takes the first of them, or the smoothed "
 936        f"maximum, after which the smoothed curve gains at most {selection.delta:g} over the next "
 937        f"{selection.window} lengths: the shortest equation that has stopped improving. {reached}\n",
 938        _crater_note(report, selection),
 939        "The following paired analysis compares every length against E3-Valid, dataset by "
 940        "dataset; it is a sensitivity analysis rather than an additional selector:\n",
 941        _table(report.length_choice) + "\n",
 942    ]
 943    return "\n".join(lines)
 944
 945
 946def _form_stability(report: EquationReport) -> float:
 947    """Mean share of LODO folds that re-select the equation's own terms."""
 948    if report.stability is None or not report.equation.terms:
 949        return float("nan")
 950    frequency = dict(zip(report.stability["term"].to_list(), report.stability["frequency"].to_list(), strict=True))
 951    return float(np.mean([frequency.get(term.name, 0.0) for term in report.equation.terms]))
 952
 953
 954def _capability_note(report: Report) -> str:
 955    """The configured capability/length comparison at the raw floor maximum."""
 956    capability = report.e3_capability
 957    if not capability.equation.terms:
 958        return ""
 959    published = float(report.e3.in_sample["r2"])
 960    reached = float(capability.in_sample["r2"])
 961    lines = [
 962        f"The published equation uses arity {report.e3.arity}. The same configuration under arity "
 963        f"{capability.arity}, at the raw maximum of its worst-protocol curve (E3-MAX), reaches "
 964        f"{len(capability.equation.terms)} terms at R² {reached:.4f} under IS:\n",
 965        "| | terms | IS | LODO | LOMO | DHO | form stability |",
 966        "|---|---|---|---|---|---|---|",
 967        f"| E3-Valid (arity {report.e3.arity}) | {len(report.e3.equation.terms)} | {published:.4f} | "
 968        f"{float(report.e3.cross_validated['loo_dataset']['r2']):.4f} | "
 969        f"{float(report.e3.cross_validated['loo_model']['r2']):.4f} | "
 970        f"{float(report.e3.cross_validated['loo_cell']['r2']):.4f} | {_form_stability(report.e3):.2f} |",
 971        f"| E3-MAX (arity {capability.arity}) | {len(capability.equation.terms)} | {reached:.4f} | "
 972        f"{float(capability.cross_validated['loo_dataset']['r2']):.4f} | "
 973        f"{float(capability.cross_validated['loo_model']['r2']):.4f} | "
 974        f"{float(capability.cross_validated['loo_cell']['r2']):.4f} | {_form_stability(capability):.2f} |",
 975        "",
 976        "This is a **capability measurement, not a recommendation**. It answers the question "
 977        "the published equation cannot answer about itself — whether the additive form is out "
 978        f"of room or whether this equation is short of it — and the answer is that {reached - published:+.4f} "
 979        "of IS R² separates the equations under the configured arities and length rules. Form stability "
 980        "is the mean share of LODO folds that re-select an equation's own terms when selection "
 981        "is re-run without the held-out dataset.\n",
 982    ]
 983    return "\n".join(lines)
 984
 985
 986def _reach_note(report: Report) -> str:
 987    """The grammar's implied ceiling ladder, and where the equation lands on it.
 988
 989    Everything here is arithmetic over `analysis.grammar_ceiling` and the fitted R2. The
 990    point of the section is that the three rungs are computable *before* any search runs, so
 991    the equation's score can be read against an expectation rather than against nothing.
 992    """
 993    if report.ceiling.height == 0:
 994        return ""
 995    ladder = report.ceiling.to_dicts()[0]
 996    raw = float(ladder["r2_raw_additive"])
 997    best = float(ladder["r2_best_per_feature"])
 998    single = float(ladder["r2_all_single_feature"])
 999    fitted = float(report.e3.in_sample["r2"])
1000    n_single = int(float(ladder["n_single_feature_terms"]))
1001    n_raw = int(report.reach.filter(pl.col("r2_raw").is_finite()).height)
1002    n_terms = len(report.e3.equation.terms)
1003
1004    strongest = report.reach.head(1).to_dicts()[0] if report.reach.height else None
1005    lines = [
1006        "Three levels of what the vocabulary can explain, each a least-squares fit over the "
1007        "library and each computable before the search runs. They bound a *sum of "
1008        "per-feature functions*, which is a different question from the additive mean-based reference "
1009        "above: that one bounds a per-dataset value plus a per-model value.\n",
1010        "| level | terms | R² |",
1011        "|---|---|---|",
1012        f"| every raw term admitted by the grammar | {n_raw} | {raw:.4f} |",
1013        f"| the best single-feature term per feature | {len(ALL_FEATURES)} | {best:.4f} |",
1014        f"| every single-feature term at once | {n_single} | {single:.4f} |",
1015        f"| **the fitted equation (E3-Valid)** | **{n_terms}** | **{fitted:.4f}** |",
1016        "",
1017    ]
1018    if strongest is not None:
1019        lines.append(
1020            f"No individual feature carries much: the strongest is `{strongest['feature']}` at "
1021            f"R² {float(strongest['r2_best']):.3f}, so any accuracy beyond that is combination "
1022            "rather than a single dominant driver. Taking one best admissible term per "
1023            f"feature is worth {best - raw:+.3f} over the admissible raw-term fit.\n"
1024        )
1025    verdict = (
1026        f"E3-Valid reaches {fitted:.4f} with {n_terms} terms, **above** the {single:.4f} that all "
1027        f"{n_single} single-feature terms reach together. An equation cannot pass that level "
1028        "by describing features one at a time, so the excess is what the cross-feature terms "
1029        "buy — the same conclusion the additive mean-based reference reaches, by an independent route."
1030        if fitted > single
1031        else f"E3-Valid reaches {fitted:.4f} with {n_terms} terms against the {single:.4f} available "
1032        f"from all {n_single} single-feature terms, so on this configuration its accuracy is "
1033        "still within what per-feature description alone could explain."
1034    )
1035    lines.append(verdict + "\n")
1036    return "\n".join(lines)
1037
1038
1039def _protocol_note(report: Report) -> str:
1040    """What the equation is worth as more of the problem becomes unfamiliar.
1041
1042    Built from the two decision tables, which now carry the protocol of every row. The point
1043    of the section is that the strictest protocol -- both the dataset and the model out of the
1044    training set -- is the one the practitioner's question actually needs, and that the trivial
1045    baselines cannot be computed under it at all.
1046    """
1047    ranking = report.ranking_baselines
1048    decision = report.decision_baselines
1049    if ranking.height == 0 or decision.height == 0:
1050        return ""
1051
1052    order = [
1053        ("in-sample", "IS", "nothing held out"),
1054        ("loo-dataset", "LODO", "the dataset unseen, the model known"),
1055        ("loo-model", "LOMO", "the model unseen, the dataset known"),
1056        ("loo-cell", "DHO", "**both unseen**"),
1057    ]
1058    rank_rows = {str(row["predictor"]): row for row in ranking.to_dicts()}
1059    mid = sorted({float(value) for value in decision["threshold"]})
1060    threshold = mid[len(mid) // 2]
1061    dec_rows = {str(row["predictor"]): row for row in decision.filter(pl.col("threshold") == threshold).to_dicts()}
1062
1063    lines = [
1064        "The ranking and the go/no-go decision are reported with **both the dataset and the "
1065        "model of every cell held out of the fit**. Neither single-group protocol answers the "
1066        "question those tasks pose: LODO has seen the learner on the other "
1067        "nineteen problems, and LOMO has seen the dataset. A recommendation is "
1068        "asked about a pair that has not been run.\n",
1069        f"| Evaluation protocol | AP | MRR | hit@1 | regret | F1 @ {threshold:g} | MCC @ {threshold:g} |",
1070        "|---|---|---|---|---|---|---|",
1071    ]
1072    for key, label, shown in order:
1073        name = next((n for n in rank_rows if key in n), None)
1074        dname = next((n for n in dec_rows if key in n), None)
1075        if name is None or dname is None:
1076            continue
1077        r, d = rank_rows[name], dec_rows[dname]
1078        emphasis = "**" if key == "loo-cell" else ""
1079        lines.append(
1080            f"| {emphasis}{label}{emphasis} — {shown} | {float(r['ap']):.3f} | {float(r['mrr']):.3f} | "
1081            f"{float(r['hit_at_1']):.2f} | {float(r['regret']):.3f} | "
1082            f"{float(d['f1']):.3f} | {float(d['mcc']):.3f} |"
1083        )
1084    lines.append("")
1085
1086    baselines = [str(name) for name in rank_rows if "per-model" in name]
1087    if baselines:
1088        best = max(baselines, key=lambda n: float(rank_rows[n]["ap"]))
1089        lines.append(
1090            "The per-model mean and median predictors are reported below under "
1091            "LODO, the only held-out protocol in which the test model still "
1092            "has training rows. **Under the strictest protocol they cannot be computed at "
1093            "all**: a model held out of every fold has no rows to average, so "
1094            f'"how well does this model usually do" has no value. The best of them reaches AP '
1095            f"{float(rank_rows[best]['ap']):.3f} and F1 {float(dec_rows[best]['f1']):.3f} while "
1096            "being shown the model identity the strictest row of the equation is denied.\n"
1097        )
1098    return "\n".join(lines)
1099
1100
1101def _ranking_verdict(baselines: pl.DataFrame) -> str:
1102    """State the ranking comparison from the paired test, never from the means.
1103
1104    The means over twenty datasets are the trap this study has fallen into three times.
1105    hit@1 moves only in steps of 0.05 on twenty folds, so one dataset flipping its top pick
1106    shifts it further than the gaps in the table; and average precision differences of this
1107    size sit well inside the bootstrap interval. The sentence is therefore built from
1108    `ap_vs_e3_significant`, not from which row has the larger mean.
1109    """
1110    if "ap_vs_e3_significant" not in baselines.columns:
1111        return ""
1112    rivals = [row for row in baselines.to_dicts() if row.get("ap_vs_e3_significant") is not None]
1113    if not rivals:
1114        return ""
1115    beaten = [row for row in rivals if row["ap_vs_e3_significant"]]
1116    names = ", ".join(_published_protocol_label(row["predictor"]) for row in rivals)
1117    if beaten:
1118        return (
1119            "Paired over the datasets, the equation differs significantly from: "
1120            + ", ".join(_published_protocol_label(row["predictor"]) for row in beaten)
1121            + ". The remaining comparisons are ties.\n"
1122        )
1123    return (
1124        f"**None of these differences survives a paired test.** Against {names} the sign test "
1125        "and the bootstrap interval over datasets both include zero, so on ranking the equation "
1126        "is indistinguishable from ordering the models by how well they usually do. Read the "
1127        "means in the table above as ties, not as a ranking of the predictors — including where "
1128        "a baseline's mean is the larger one.\n"
1129    )
1130
1131
1132def _baseline_centre_note(baselines: pl.DataFrame) -> str:
1133    """Which centre is the harder baseline, per metric, derived rather than asserted.
1134
1135    The mean minimises squared error and the median minimises absolute error, so the fair
1136    opponent differs by metric. Rather than state that as a claim, the sentence is built from
1137    the table: for each metric, the best mean row is compared with the best median row.
1138    """
1139    rows = {str(row["baseline"]): row for row in baselines.to_dicts()}
1140    pairs = [(name, name.replace(" mean ", " median ")) for name in rows if " mean " in name]
1141    if not pairs:
1142        return ""
1143
1144    def best(names: list[str], metric: str, lower_is_better: bool) -> tuple[str, float]:
1145        scored = [(name, float(rows[name][metric])) for name in names]
1146        return min(scored, key=lambda item: item[1]) if lower_is_better else max(scored, key=lambda item: item[1])
1147
1148    means = [mean for mean, _ in pairs]
1149    medians = [median for _, median in pairs if median in rows]
1150    if not medians:
1151        return ""
1152
1153    lines = [
1154        "Every trivial predictor appears at its mean and at its median, because the metrics "
1155        "disagree about which is the honest opponent: the mean minimises squared error and "
1156        "the median minimises absolute error, so an MAE quoted against a mean baseline is "
1157        "quoted against a predictor that is not minimising the metric it is judged on. "
1158        "Reading the strongest baseline of each kind, per metric:\n",
1159        "| metric | strongest mean baseline | strongest median baseline | harder |",
1160        "|---|---|---|---|",
1161    ]
1162    for metric, lower_is_better in (("r2", False), ("mae", True), ("smape", True)):
1163        mean_name, mean_value = best(means, metric, lower_is_better)
1164        median_name, median_value = best(medians, metric, lower_is_better)
1165        median_wins = median_value < mean_value if lower_is_better else median_value > mean_value
1166        lines.append(
1167            f"| {metric.upper()} | {_published_protocol_label(mean_name)} ({mean_value:.4f}) | "
1168            f"{_published_protocol_label(median_name)} ({median_value:.4f}) | "
1169            f"**{'median' if median_wins else 'mean'}** |"
1170        )
1171    return "\n".join(lines) + "\n"
1172
1173
1174def _configuration(config: Configuration, selection: SelectionConfig | None) -> str:
1175    rows = [{"setting": f"search.{field.name}", "value": str(getattr(config, field.name))} for field in fields(config)]
1176    if selection is not None:
1177        rows += [
1178            {"setting": f"selection.{field.name}", "value": str(getattr(selection, field.name))}
1179            for field in fields(selection)
1180        ]
1181    return _table(pl.DataFrame(rows))
1182
1183
1184def render(
1185    report: Report,
1186    columns: dict[str, np.ndarray],
1187    truth: np.ndarray,
1188    dataset_features: tuple[str, ...],
1189    model_features: tuple[str, ...],
1190    *,
1191    frame: pl.DataFrame | None = None,
1192    config: Configuration | None = None,
1193    selection: SelectionConfig | None = None,
1194    source: str | None = None,
1195) -> str:
1196    """The full study as a markdown report.
1197
1198    Every table and every sentence is computed from ``report``. The prose that is fixed
1199    -- section headings and the explanations of what each table means -- is fixed because
1200    it describes the *method*, which does not change between runs; anything describing
1201    the *result* is generated.
1202    """
1203    rule = selection or load_config().selection
1204    importance = term_importance(report.e3.equation, columns, dataset_features, model_features, report.e3.stability)
1205    concentration = coverage(importance)
1206    equation = report.e3.equation
1207    # Counted from the corpus rather than written into the prose below it. A generated
1208    # sentence carrying a literal "20 datasets" is the same staleness the generated sections
1209    # exist to prevent, one level further in.
1210    datasets = int(frame[DATASET_COLUMN].n_unique()) if frame is not None else 0
1211
1212    parts: list[str] = []
1213    parts.append("# ml-meta-perf — fitted equation and extracted practices\n")
1214    parts.append(
1215        "Generated by `ml_meta_perf.report.render`. Every number, table and sentence below is "
1216        "computed from the fitted equation and the validation run; nothing is written by "
1217        "hand at report time.\n"
1218    )
1219
1220    parts.append("## 1. The configuration this run used\n")
1221    if source:
1222        parts.append(f"Meta-dataset: `{source}`\n")
1223    if config is not None:
1224        parts.append("Every hyperparameter below is read from `config/study.json`; E1, E2 and E3 share all of them.\n")
1225        parts.append(_configuration(config, selection) + "\n")
1226
1227    parts.append("## 1b. The corpus\n")
1228    parts.append(
1229        "What the meta-dataset is, computed from the file the run was fitted on rather than "
1230        "transcribed into prose beside it.\n"
1231    )
1232    if frame is not None:
1233        parts.append(_table(corpus_summary(frame)) + "\n")
1234        parts.append(
1235            "**The absent cells are not missing at random.** Every dataset short of models is "
1236            "one of the smallest in the corpus, which is what the instance counts show:\n"
1237        )
1238        parts.append(_table(missing_cells(frame), float_format="{:.0f}") + "\n")
1239        parts.append("And how MCC is distributed over those rows:\n")
1240        parts.append(_table(target_summary(frame)) + "\n")
1241        parts.append(
1242            "**The table counts observations exactly at 0 or 1, and those below "
1243            "zero.** That is "
1244            "what makes MAE rather than SMAPE the reported error: SMAPE divides by "
1245            "`|truth| + |prediction|`, so every row at exactly zero contributes the full 200% "
1246            "unless the prediction is exactly zero too, and the metric ends up dominated by "
1247            "the rows the equation is already known to handle worst.\n"
1248        )
1249        parts.append(
1250            "Two things these tables cannot say, both of which bound every number in the study. "
1251            "Each row is the **best of five seeds**, not their mean, so the target is "
1252            "optimistic and subject to seed variation that is not modelled separately; and "
1253            "tuning uses a sample **capped at 100,000 rows**, while later performance evaluation "
1254            "uses available source data with dataset/model-specific limits. `nr_inst` is the "
1255            "source dataset's size rather than necessarily the training set's. These are properties of the "
1256            "corpus builder upstream, and are audited against it in the prose above.\n"
1257        )
1258    else:
1259        parts.append("_(not summarised: the meta-dataset was not supplied to the renderer)_\n")
1260
1261    parts.append("## 1c. The headline\n")
1262    parts.append(
1263        "The one table the study is summarised by, so that the summary cannot drift from "
1264        f"the chapters. Every row is scored on the same {int(report.e3.in_sample['n'])} rows "
1265        "under the same protocols. E1, E2, and E3-Valid share the retained base configuration; "
1266        "E3-MAX takes the raw floor maximum at the configured capability arity; the current "
1267        "primary and capability arities both equal three.\n"
1268    )
1269    parts.append(_table(_headline(report)) + "\n")
1270    parts.append(
1271        "**Do not read these R² values as achievements against each other.** They share a "
1272        "scale but not a ceiling: E1 sees only dataset features, every row of a dataset "
1273        "shares one feature vector, and so E1 can predict nothing but a per-dataset "
1274        "constant. The `true dataset means` row gives E1's structural maximum. "
1275        "The last column gives only E1's fraction of that "
1276        "bound. E2's model-mean reference is descriptive because Processing Units Number "
1277        "can vary with dataset shape; no ceiling fraction is reported for E2.\n"
1278    )
1279    parts.append(
1280        "**E2, E3-Valid, and E3-MAX have blank ceiling cells.** Their recorded predictors are "
1281        "not constant within the corresponding identity groups. The comparison table in "
1282        "chapter 5 also shows the additive mean-based reference. On the incomplete grid that "
1283        "reference is not a jointly fitted least-squares additive bound; mixed terms can "
1284        "represent interactions beyond marginal group effects.\n"
1285    )
1286
1287    parts.append("## 2. The equation\n")
1288    parts.append(
1289        f"E3-Valid uses **{equation.n_terms} terms** over dataset and model meta-features, "
1290        "simplified and refitted after pruning, so it evaluates exactly as printed.\n"
1291    )
1292    parts.append("```\n" + str(equation) + "\n```\n")
1293    parts.append("LaTeX:\n")
1294    parts.append("```latex\n" + equation.to_latex() + "\n```\n")
1295
1296    parts.append("## 3. How well it does\n")
1297    parts.append("| protocol | R² | MAE | RMSE | n |")
1298    parts.append("|---|---|---|---|---|")
1299    parts.append(_scores("IS", report.e3.in_sample))
1300    protocol_labels = {"loo_dataset": "LODO", "loo_model": "LOMO", "loo_cell": "DHO"}
1301    for label, scores in report.e3.cross_validated.items():
1302        parts.append(_scores(protocol_labels.get(label, label.replace("_", "-")), scores))
1303    parts.append("")
1304    parts.append(
1305        "Cross-validated rows hold out a whole dataset or a whole model, so the equation "
1306        "is scored on a group it has never seen. That is the number that matters, and it "
1307        "is well below the IS value at this sample size.\n"
1308    )
1309    parts.append("Against the baselines and descriptive group-mean references:\n")
1310    parts.append(_table(report.comparison) + "\n")
1311
1312    parts.append("### The trivial predictors, at both centres\n")
1313    parts.append(_table(report.baselines) + "\n")
1314    parts.append(_baseline_centre_note(report.baselines))
1315
1316    parts.append("## 3c. How far the form could reach\n")
1317    parts.append(
1318        "The following library fits and length comparison describe representational reach "
1319        "under the recorded corpus and settings. They are not independent performance "
1320        "estimates or universal ceilings.\n"
1321    )
1322    parts.append("### How far the additive form reaches\n")
1323    parts.append(_capability_note(report))
1324
1325    parts.append("### What the vocabulary could reach, before any search\n")
1326    parts.append(_reach_note(report))
1327
1328    parts.append("### Where the variance is, before any equation\n")
1329    parts.append("Variance of MCC explained by group identity alone, with no equation involved:\n")
1330    parts.append(_table(report.decomposition) + "\n")
1331
1332    parts.append("### How fast interaction pays\n")
1333    parts.append(
1334        "The additive form cannot represent dataset-by-model interaction beyond what its mixed "
1335        "terms reach. The ladder below adds interaction components to an oracle that is "
1336        "handed the true group means, so it measures the ceiling rather than any equation:\n"
1337    )
1338    parts.append(_table(report.oracles) + "\n")
1339    parts.append(
1340        "That is a ceiling, not a score. Whether the equation reaches any of it is a separate "
1341        "question, and the answer is that it reaches some: below, `alignment` is the squared "
1342        "correlation between the equation's own interaction residual and the leading components "
1343        "of the oracle's, over observed cells. `leading_share` is how much of the interaction "
1344        "variance those components carry, and `interaction_share` how much of MCC's variance is "
1345        "interaction at all.\n"
1346    )
1347    parts.append(_table(report.interaction) + "\n")
1348
1349    parts.append("## 3b. Why a subset rather than every term\n")
1350    parts.append(
1351        "The control for the whole selection stage. If handing every candidate term to "
1352        "unpenalised least squares in one go transferred well, the beam search and the length "
1353        "rule would be machinery in search of a problem.\n"
1354    )
1355    parts.append(_table(report.saturated) + "\n")
1356    parts.append(_saturated_note(report))
1357
1358    parts.append("## 4. Equation analysis\n")
1359    parts.append(
1360        f"The equation has {concentration['n_terms']} terms, of which "
1361        f"**{concentration['n_major']}** carry {float(concentration['major_share']):.0%} of the "
1362        f"standardised weight mass; the single largest carries "
1363        f"{float(concentration['top_share']):.1%}, and the weights behave like "
1364        f"**{float(concentration['effective_terms']):.1f} equally-weighted terms** "
1365        "(inverse Simpson index of the shares).\n"
1366    )
1367    flatness = float(concentration["effective_terms"]) / max(int(concentration["n_terms"]), 1)
1368    parts.append(
1369        "That last number is the one to read for concentration, because it does not depend "
1370        "on where a threshold is drawn. At "
1371        f"{flatness:.0%} of the term count the equation is "
1372        + (
1373            "**flat**: no single term dominates. That is a statement about the *unit of "
1374            "explanation*, not about the quality of the equation — MCC here is inferred by "
1375            "a set of terms acting together rather than by one or two that could be quoted "
1376            "on their own. Three readings follow, and the sections below give each one: "
1377            "read the terms in the blocks that move together, read which features the "
1378            "search reached for, and read which operations it needed to apply to them.\n"
1379            if flatness > 0.6
1380            else "**concentrated**: a minority of terms does most of the work, and reading "
1381            "those few is close to reading the whole equation.\n"
1382        )
1383    )
1384    parts.append(
1385        "`beta` is the standardised weight — the MCC contributed per standard deviation of "
1386        "the term, which is what makes terms in unrelated units comparable. `effect` is the "
1387        "swing in predicted MCC across the middle 80% of the term's observed range. "
1388        "`stability` is the fraction of LODO folds that selected the term.\n"
1389    )
1390    parts.append(_table(importance) + "\n")
1391
1392    shown = min(SENTENCE_LIMIT, int(concentration["n_major"]))
1393    parts.append(f"### The {shown} largest terms, in words\n")
1394    if int(concentration["n_major"]) > shown:
1395        parts.append(
1396            f"({int(concentration['n_major'])} terms are flagged major; the leading {shown} are "
1397            "written out, and the table above carries the rest.)\n"
1398        )
1399    for index, sentence in enumerate(term_sentences(importance, limit=SENTENCE_LIMIT), start=1):
1400        parts.append(f"{index}. {sentence}")
1401    parts.append("")
1402
1403    unstable = unstable_majors(importance)
1404    parts.append("### Large terms the folds disagreed on\n")
1405    if unstable.height == 0:
1406        parts.append("None: every major term was selected by at least half of the LODO folds.\n")
1407    else:
1408        parts.append(
1409            f"{unstable.height} of the major terms were selected by fewer than half of the "
1410            "folds. A large weight and a low selection frequency together mean the term is "
1411            "doing its work for *this* training set and would be replaced by something else "
1412            "on another, which the equation as printed does not show. **Do not build "
1413            "guidance on these.**\n"
1414        )
1415        parts.append(_table(unstable) + "\n")
1416
1417    blocks = term_groups(report.e3.equation, columns, importance)
1418    parts.append("### Reading the terms in blocks\n")
1419    parts.append(
1420        "An additive form invites reading one term at a time, and that works when one or "
1421        "two weights dominate. When they do not, the honest unit is larger than a term and "
1422        "smaller than the equation: terms whose per-row contributions move together say the "
1423        "same thing about a row and can be read as one block. Grouping is on the "
1424        "contributions rather than on shared features, because two terms can share no "
1425        "feature and still track each other.\n"
1426    )
1427    parts.append(
1428        f"{blocks.height} blocks over {equation.n_terms} terms"
1429        + (
1430            f", the largest holding {int(blocks['n_terms'][0])} terms and "
1431            f"{float(blocks['share'][0]):.0%} of the weight mass.\n"
1432            if blocks.height
1433            else ".\n"
1434        )
1435    )
1436    parts.append(_table(blocks) + "\n")
1437
1438    parts.append("### Which features the search reached for\n")
1439    usage = feature_usage(report.e3.equation, importance, dataset_features + model_features)
1440    used = usage.filter(pl.col("n_terms") > 0)
1441    parts.append(
1442        f"**{used.height} of {usage.height}** available meta-features appear in the "
1443        "equation. `share` sums the weight mass of every term a feature appears in, so a "
1444        "feature in two terms is credited both and the column does not sum to 1 — it "
1445        "answers how much of the equation touches this feature, not how much it owns. A "
1446        "feature the search declined to use after seeing every transform of it is itself a "
1447        "result.\n"
1448    )
1449    parts.append(_table(usage) + "\n")
1450
1451    parts.append("### Which operations the equation needed\n")
1452    operations = operation_usage(report.e3.equation, importance, config.max_arity if config is not None else None)
1453    parts.append(
1454        "The vocabulary offers five operations and five transforms and the search is free "
1455        "to ignore any of them, so a row that was offered and went unused is a shape this "
1456        "data turned out not to need. Rows marked `offered = no` were kept out of the "
1457        "library by the arity cap and say nothing about the data:\n"
1458    )
1459    parts.append(_table(operations) + "\n")
1460
1461    parts.append("### Where the equation's variance comes from\n")
1462    for sentence in group_sentences(report.shares):
1463        parts.append(f"- {sentence}")
1464    parts.append("")
1465    parts.append(_table(report.shares) + "\n")
1466
1467    parts.append("## 4b. The ceiling on model descriptors\n")
1468    parts.append(
1469        "Under LODO every model appears in every training fold, so the "
1470        "equation's residual can be averaged per model on the training rows and applied to "
1471        "the held-out dataset using fitting histories. These model-identity corrections are "
1472        "conditional diagnostic references, not universal ceilings on better descriptors. "
1473        "They are unavailable for a learner whose entire history is excluded.\n"
1474    )
1475    parts.append(_table(report.identity) + "\n")
1476    parts.append(_identity_note(report))
1477
1478    parts.append("## 5. Best practices\n")
1479    parts.append(
1480        "A best practice is general, transferable advice that already circulates in the "
1481        "field — not a property of this equation. So the practices below are taken from the "
1482        "literature and this study is used to *weigh* them: each verdict, and the numbers "
1483        "inside it, are computed from this run by `ml_meta_perf.guidance`, against a stated "
1484        "threshold, so other data can overturn any of them.\n"
1485    )
1486    if frame is not None:
1487        verdicts = assess(frame, report)
1488        parts.append(render_guidance(verdicts) + "\n")
1489        parts.append("At a glance:\n")
1490        parts.append(_table(as_guidance_table(verdicts).select("practice", "verdict", "magnitude")) + "\n")
1491    else:
1492        parts.append("_(not assessed: the meta-dataset was not supplied to the renderer)_\n")
1493
1494    parts.append("### Each practice against the equation's own terms\n")
1495    parts.append(
1496        "The verdicts above are drawn from corpus averages -- family means, variance shares, "
1497        "paired tests -- which any study with this corpus could compute. This table asks the "
1498        "stronger question, and the one an interpretability-first study is uniquely able to "
1499        "ask: **which of the equation's terms carries this practice, with what strength and "
1500        "which sign?**\n"
1501    )
1502    if frame is not None:
1503        parts.append(_table(equation_evidence(report, columns)) + "\n")
1504        parts.append(
1505            "**One row per (practice, term) pair**, because the equation is a sum of terms and "
1506            "a term is the unit a practice can be held against. `expected` is what the "
1507            "practice predicts as the feature rises; `direction` is what *that term's own "
1508            "contribution* does, measured on the data rather than read off the weight sign -- "
1509            "which would be wrong the moment the feature sits in a denominator, and several "
1510            "here do. `not selected` marks a claim resting on a feature the search never took, "
1511            "so the equation is silent on it rather than supporting it.\n"
1512        )
1513        parts.append(_coverage_note(report, columns))
1514    else:
1515        parts.append("_(not checked: the meta-dataset was not supplied to the renderer)_\n")
1516
1517    parts.append("## 5b. The measurements underneath\n")
1518    parts.append(
1519        "What the fitted equation says about each raw feature it uses, kept only when the "
1520        "feature moves predicted MCC enough to matter, does so monotonically enough for a "
1521        "sentence to be true of it, and does so through terms that survived most folds. "
1522        "Directions are measured on the data rather than read off weight signs, because a "
1523        "feature can appear in several terms and inside denominators.\n"
1524    )
1525    parts.append(
1526        f"**These are associations across {datasets} datasets, not causal claims, and not practices "
1527        "on their own** — a statement about a meta-feature column is a measurement. Section "
1528        "5 is where they become advice, by supporting or failing to support something a "
1529        "practitioner could already have been told.\n"
1530    )
1531    parts.append(render_practices(report.practices) + "\n")
1532    parts.append("Evidence:\n")
1533    parts.append(
1534        _table(
1535            report.practices.select("feature", "meaning", "n_terms", "direction", "effect", "stability", "confidence")
1536            if report.practices.height
1537            else report.practices
1538        )
1539        + "\n"
1540    )
1541
1542    parts.append("### These are conditional statements, not marginal ones\n")
1543    conditional = marginal_versus_conditional(report.practices, columns, truth)
1544    agree = int(conditional["agrees"].sum()) if conditional.height else 0
1545    parts.append(
1546        "A practice states what the *equation* does as a feature rises, with every other "
1547        "term present. A marginal correlation states what the feature does alone. They are "
1548        f"different quantities, and here **{agree} of {conditional.height} agree** on the "
1549        "sign:\n"
1550    )
1551    parts.append(_table(conditional) + "\n")
1552    parts.append(
1553        "Disagreement is what conditioning does, not a defect. A marginal correlation mixes "
1554        "a feature's effect with everything it travels with; inside the equation the terms "
1555        "carrying those companions are already present, so what is left for this feature is "
1556        "what it adds beyond them. The practical consequence: **these statements describe "
1557        "what to expect once the other factors are accounted for, not what a scatter plot of "
1558        "that one feature will show** — and the scatter plot is what a reader will "
1559        "accidentally check against.\n"
1560    )
1561
1562    parts.append("## 5c. Reading a single prediction\n")
1563    parts.append(
1564        "The same equation on one row, term by term. The column adds up, and a reader can "
1565        "check that it does, which is the interpretability payoff in its most direct form "
1566        "-- and the row is "
1567        "chosen arithmetically, at the equation's **median absolute error**, so the example "
1568        "is a typical prediction rather than a flattering one somebody picked.\n"
1569    )
1570    if frame is not None:
1571        parts.append(single_prediction(equation, columns, truth, frame) + "\n")
1572    else:
1573        parts.append("_(not shown: the meta-dataset was not supplied to the renderer)_\n")
1574
1575    parts.append("## 6. Acting on it\n")
1576    parts.append(
1577        "R² is the wrong question for a practitioner, who asks whether a model will work on "
1578        "some data rather than what its MCC will be to three decimals. Thresholding both "
1579        "the prediction and the truth turns the equation into a go/no-go rule, scored here "
1580        "on held-out datasets:\n"
1581    )
1582    parts.append(_table(report.decision) + "\n")
1583    parts.append(
1584        "`majority` is the floor any such rule has to clear. The harder comparison is a "
1585        'predictor that answers "how well does this model usually do", thresholded the same '
1586        "way — at both centres, for the reason the error metrics report both:\n"
1587    )
1588    parts.append(_table(report.decision_baselines) + "\n")
1589
1590    parts.append(_protocol_note(report))
1591
1592    parts.append("Ranking models within a held-out dataset:\n")
1593    regret = float(np.mean(report.selection["regret"].to_numpy()))
1594    parts.append(
1595        f"- mean top-1 regret **{regret:.3f}** MCC — what you give up by taking the model the equation ranks first\n"
1596    )
1597    parts.append(_table(report.ranking_baselines) + "\n")
1598    parts.append(_ranking_verdict(report.ranking_baselines))
1599
1600    parts.append("## 6b. What an opaque model reaches, and does not\n")
1601    parts.append(
1602        "The other side of the trade, priced. Three standard regressors on all eighteen "
1603        "raw corpus columns, under the same protocols, with the same clip to the training "
1604        "fold's range that every reported number uses.\n"
1605    )
1606    parts.append(_table(report.opaque) + "\n")
1607    parts.append(_opaque_note(report))
1608
1609    parts.append("## 7. Why a random split is not a protocol\n")
1610    parts.append(
1611        "The same equation under three splits. A random k-fold puts rows of one dataset on "
1612        "both sides of the fold, so it does not test transfer to an unseen dataset. The "
1613        "difference between that row and the grouped splits measures how much this potential "
1614        "leakage changes the score; here random k-fold and LODO are essentially "
1615        "identical:\n"
1616    )
1617    parts.append(_table(report.leakage) + "\n")
1618
1619    parts.append("## 8. Equation length\n")
1620    parts.append(_length_note(report, rule))
1621    parts.append("The full curve reports all four protocols at every length; the length rule reads their minimum:\n")
1622    parts.append(_table(report.e3.curve) + "\n")
1623
1624    parts.append("## 8b. How much of the transfer is the form's selection\n")
1625    parts.append(_optimism_note(report))
1626    parts.append("## 9. The dataset-only and model-only controls\n")
1627    parts.append(
1628        "E1 sees dataset meta-features only, so it can predict just one value per dataset; "
1629        "E2 sees model meta-features only. Together they show how much of MCC each half of "
1630        "the meta-data explains on its own.\n"
1631    )
1632    parts.append(f"**E1** ({report.e1.equation.n_terms} terms):\n")
1633    parts.append("```\n" + str(report.e1.equation) + "\n```\n")
1634    parts.append(f"**E2** ({report.e2.equation.n_terms} terms):\n")
1635    parts.append("```\n" + str(report.e2.equation) + "\n```\n")
1636    parts.append(_shared_terms_note(report))
1637
1638    if report.e3_capability is report.e3:
1639        parts.append("### E3-Valid and E3-MAX coincide\n")
1640        parts.append(
1641            f"The length rule and the capability rule select the same arity-{report.e3.arity}, "
1642            f"{report.e3.equation.n_terms}-term equation. The full "
1643            "equation above therefore serves both roles; repeating it here would create a "
1644            "duplicate result.\n"
1645        )
1646    else:
1647        parts.append("### The capability variant, in full\n")
1648        parts.append(
1649            "The same feature sets under the most capable searched grammar. This variant is "
1650            "reported so that the published equation's accuracy can be read against what the "
1651            "wider additive form can do, rather than only against oracles and baselines.\n"
1652        )
1653        capability_scores = report.e3_capability.cross_validated
1654        parts.append(
1655            f"It reaches **{float(report.e3_capability.in_sample['r2']):.4f}** under IS against "
1656            f"the published equation's {float(report.e3.in_sample['r2']):.4f}, and "
1657            f"**{float(capability_scores['loo_dataset']['r2']):.4f}** under LODO "
1658            f"against {float(report.e3.cross_validated['loo_dataset']['r2']):.4f}.\n"
1659        )
1660        parts.append(f"**E3-MAX** ({report.e3_capability.equation.n_terms} terms):\n")
1661        parts.append("```\n" + str(report.e3_capability.equation) + "\n```\n")
1662        parts.append("LaTeX:\n")
1663        parts.append("```latex\n" + report.e3_capability.equation.to_latex() + "\n```\n")
1664
1665    return "\n".join(parts)
1666
1667
1668#: Which generated section belongs in which chapter, by the ``## N. Title`` heading `render`
1669#: gives it. The study has six chapters and its results belong *in* them: an earlier layout
1670#: put every generated table in a separate chapter 10 and had the written chapters point at
1671#: it, which kept the numbers from going stale and left the results dispersed across two
1672#: places a reader had to hold at once.
1673CHAPTER_SECTIONS: dict[str, tuple[str, ...]] = {
1674    "index.md": ("1c. The headline",),
1675    "01-dataset.md": ("1b. The corpus",),
1676    # Ownership, one topic to one chapter. The length rule lives with the selection procedure
1677    # that applies it, the ceilings with the equation they bound, the protocols with the
1678    # evaluation. Each generated section has one chapter owner.
1679    "03-term-selection.md": (
1680        "1. The configuration this run used",
1681        "3b. Why a subset rather than every term",
1682        "8. Equation length",
1683        "8b. How much of the transfer is the form's selection",
1684    ),
1685    "04-equation.md": (
1686        "2. The equation",
1687        "3c. How far the form could reach",
1688        "4. Equation analysis",
1689        "4b. The ceiling on model descriptors",
1690        "9. The dataset-only and model-only controls",
1691    ),
1692    "05-evaluation.md": (
1693        "3. How well it does",
1694        "6. Acting on it",
1695        "6b. What an opaque model reaches, and does not",
1696        "7. Why a random split is not a protocol",
1697    ),
1698    "06-practices.md": ("5. Best practices", "5b. The measurements underneath", "5c. Reading a single prediction"),
1699}
1700
1701#: The markers a generated block sits between. Everything between them is replaced on every
1702#: run; everything outside them is written by hand and never touched.
1703BEGIN, END = "<!-- generated: do not edit below -->", "<!-- end generated -->"
1704#: The generated section the README carries between its markers.
1705README_SECTION = "1c. The headline"
1706
1707
1708def sections(text: str) -> dict[str, str]:
1709    """Split a rendered report into its ``## N. Title`` sections, keyed by that title."""
1710    found: dict[str, str] = {}
1711    title, body = None, []
1712    for line in text.split("\n"):
1713        if line.startswith("## "):
1714            if title is not None:
1715                found[title] = "\n".join(body).strip()
1716            title, body = line[3:].strip(), []
1717        elif title is not None:
1718            body.append(line)
1719    if title is not None:
1720        found[title] = "\n".join(body).strip()
1721    return found
1722
1723
1724def _demote(body: str) -> str:
1725    """Push a section's own subheadings down one level, so it nests under a chapter heading."""
1726    return "\n".join(("#" + line) if line.startswith("###") else line for line in body.split("\n"))
1727
1728
1729def splice(page: str, block: str) -> str:
1730    """Replace the generated block in ``page``, or insert one if it has none.
1731
1732    A new block goes *before* the chapter's closing limitations section, so a chapter reads
1733    method, then results, then caveats -- the order a paper is read in. Replacing an existing
1734    block leaves it wherever it already is, so the position survives a hand edit.
1735    """
1736    stamped = f"{BEGIN}\n\n{block.strip()}\n\n{END}"
1737    if BEGIN in page and END in page:
1738        head, rest = page.split(BEGIN, 1)
1739        _, tail = rest.split(END, 1)
1740        return head + stamped + tail
1741    for line in page.split("\n"):
1742        if line.startswith("## Limitations"):
1743            head, tail = page.split(line, 1)
1744            return head.rstrip() + "\n\n" + stamped + "\n\n" + line + tail
1745    return page.rstrip() + "\n\n" + stamped + "\n"
1746
1747
1748def write_into_chapters(
1749    report: Report,
1750    columns: dict[str, np.ndarray],
1751    truth: np.ndarray,
1752    dataset_features: tuple[str, ...],
1753    model_features: tuple[str, ...],
1754    docs: str | Path,
1755    *,
1756    frame: pl.DataFrame | None = None,
1757    config: Configuration | None = None,
1758    selection: SelectionConfig | None = None,
1759    source: str | None = None,
1760    readme: str | Path | None = None,
1761) -> list[Path]:
1762    """Write each generated section into the chapter it belongs to.
1763
1764    The results live in the chapters that discuss them rather than in a separate generated
1765    chapter, and they are still generated rather than narrated: everything between `BEGIN` and
1766    `END` is replaced on every run, and everything outside is hand-written prose that the
1767    pipeline never touches. A chapter therefore cannot carry a stale table, which is the
1768    property the separate report chapter existed to guarantee.
1769    """
1770    text = render(
1771        report,
1772        columns,
1773        truth,
1774        dataset_features,
1775        model_features,
1776        frame=frame,
1777        config=config,
1778        selection=selection,
1779        source=source,
1780    )
1781    available = sections(text)
1782    folder = Path(docs)
1783    written: list[Path] = []
1784    for filename, titles in CHAPTER_SECTIONS.items():
1785        page = folder / filename
1786        if not page.is_file():
1787            continue
1788        blocks: list[str] = []
1789        for title in titles:
1790            body = available.get(title)
1791            if body:
1792                # Drop the "N. " ordinal: the chapter supplies the position, and a generated
1793                # heading numbered against the report would contradict the chapter it sits in.
1794                heading = title.split(". ", 1)[-1]
1795                blocks.append(f"## {heading}\n\n{_demote(body)}")
1796        if blocks:
1797            content = page.read_text(encoding="utf-8")
1798            page.write_text(
1799                splice(content, "\n\n".join(blocks)),
1800                encoding="utf-8",
1801            )
1802            written.append(page)
1803    # The README carries the headline too, so the first numbers a reader sees are generated
1804    # like every other number. It must already hold the markers: a README without them is
1805    # left alone rather than given a block in a place nobody chose.
1806    if readme is not None and Path(readme).is_file():
1807        page = Path(readme)
1808        content = page.read_text(encoding="utf-8")
1809        body = available.get(README_SECTION)
1810        if body and BEGIN in content and END in content:
1811            heading = README_SECTION.split(". ", 1)[-1]
1812            page.write_text(splice(content, f"## {heading}\n\n{_demote(body)}"), encoding="utf-8")
1813            written.append(page)
1814    return written
1815
1816
1817def glossary() -> pl.DataFrame:
1818    """The feature glossary as a table, for appending to a report or a paper.
1819
1820    Restricted to `ALL_FEATURES`, which is the generated corpus schema and exactly the set
1821    of columns the equations may draw on.
1822    """
1823    return pl.DataFrame([{"feature": name, "meaning": FEATURE_GLOSSARY[name]} for name in ALL_FEATURES])