ml_meta_perf.experiment

   1"""The end-to-end experiment: fit E1, E2 and E3, validate all three, and report.
   2
   3The three equations differ only in which features they may draw on -- E1 sees dataset
   4meta-features, E2 sees model meta-features, E3 sees both -- so the gaps between them
   5measure what each half of the meta-data is worth. Everything else is shared: **one
   6configuration** (`config.Configuration`, read from ``config/study.json``), one search, one
   7set of protocols, and **one length rule**, `selection.plateau_knee`. Each equation searches
   8every feature it is allowed; the search, not the configuration, decides which it keeps.
   9
  10E3-Valid is the E3 the study publishes. E3-MAX is the same configuration at the configured
  11capability arity (``selection.capability_arity``) at the raw maximum of its worst-protocol curve: a
  12capability bound on the additive form, not an equation put forward to be read.
  13
  14Study chapter: [4. The equation][study-chapter] -- the rationale, in
  15prose, with the figures.
  16
  17[study-chapter]: https://github.com/mariolpantunes/ml-meta-perf/blob/main/assets/docs/04-equation.md
  18"""
  19
  20from __future__ import annotations
  21
  22import dataclasses
  23from dataclasses import dataclass, field
  24
  25import numpy as np
  26import polars as pl
  27
  28from ml_meta_perf.analysis import feature_reach, grammar_ceiling, saturated_fit, screen
  29from ml_meta_perf.attribution import group_shares, term_effects, variance_decomposition
  30from ml_meta_perf.config import Configuration, OpaqueConfig, SelectionConfig, load_config
  31from ml_meta_perf.data import (
  32    DATASET_COLUMN,
  33    DATASET_FEATURES,
  34    MODEL_COLUMN,
  35    MODEL_FEATURES,
  36    columns_as_arrays,
  37    groups,
  38    load,
  39    target,
  40)
  41from ml_meta_perf.identity import correct_out_of_fold
  42from ml_meta_perf.model import Equation
  43from ml_meta_perf.opaque import Builder as OpaqueBuilder
  44from ml_meta_perf.opaque import OpaqueRun
  45from ml_meta_perf.opaque import estimators as opaque_estimators
  46from ml_meta_perf.opaque import evaluate as opaque_evaluate
  47from ml_meta_perf.practices import best_practices, feature_practices
  48from ml_meta_perf.search import prune, search
  49from ml_meta_perf.selection import complexity, floor_argmax, floor_curve, plateau_knee, protocol_spread, smoothed
  50from ml_meta_perf.stats import mae, r2_score
  51from ml_meta_perf.terms import Library, build_library
  52from ml_meta_perf.validate import (
  53    CrossValidation,
  54    Scores,
  55    additive_mean_reference,
  56    baseline_group_centre,
  57    cross_validate_doubly_held_out,
  58    cross_validate_fixed_form,
  59    decision_report,
  60    fold_selections,
  61    interaction_capture,
  62    leave_one_group_out,
  63    library_rows,
  64    oracle_ladder,
  65    paired_comparison,
  66    random_kfold_groups,
  67    ranking_report,
  68    score,
  69    term_stability,
  70)
  71
  72#: E1 is fitted on all 484 rows like the other two. An earlier version fitted it on the 20
  73#: aggregated per-dataset means, on the grounds that a predictor constant inside a group can
  74#: only ever predict that group's mean anyway. That was true and still the wrong choice: it put
  75#: E1's R2 on a 20-point denominator, so its headline could not be compared with E3's without a
  76#: paragraph of explanation. See [chapter 4][study-chapter].
  77
  78
  79@dataclass
  80class EquationReport:
  81    """One equation together with everything said about it."""
  82
  83    equation: Equation
  84    #: The selected length: `selection.plateau_knee` for E1, E2 and E3-Valid, the raw floor
  85    #: maximum for E3-MAX. Every output reads the result from here.
  86    n_terms: int
  87    #: The grammar this equation was searched under -- the configuration's for E1, E2 and
  88    #: E3-Valid, ``selection.capability_arity`` for E3-MAX.
  89    arity: int
  90    in_sample: dict[str, float | int]
  91    curve: pl.DataFrame
  92    cross_validated: dict[str, dict[str, float | int]] = field(default_factory=dict)
  93    stability: pl.DataFrame | None = None
  94    #: The full cross-validated path per protocol, kept so later tables can reuse it.
  95    #: ``leakage_demonstration`` and ``decision_quality`` ask for the same folds at the
  96    #: same settings, and recomputing them was a third of the study's runtime.
  97    paths: dict[str, dict[int, CrossValidation]] = field(default_factory=dict)
  98
  99
 100def _curve(
 101    sizes: tuple[int, ...],
 102    in_sample: dict[int, Scores],
 103    paths: dict[str, dict[int, CrossValidation]],
 104    truth: np.ndarray,
 105) -> pl.DataFrame:
 106    """One row per equation length, carrying every metric under every protocol.
 107
 108    IS gets the same three metrics as the cross-validated columns, not just R²:
 109    an error curve that omits the IS line cannot show how far apart fit and
 110    transfer run, which is the whole point of plotting them together.
 111    """
 112    rows: list[dict[str, object]] = []
 113    for size in sizes:
 114        if size not in in_sample:
 115            continue
 116        row: dict[str, object] = {
 117            "n_terms": size,
 118            "r2_in_sample": in_sample[size].r2,
 119            "mae_in_sample": in_sample[size].mae,
 120            "smape_in_sample": in_sample[size].smape,
 121        }
 122        for label, path in paths.items():
 123            scores = path[size].scores(truth)
 124            row[f"r2_{label}"] = scores.r2
 125            row[f"mae_{label}"] = scores.mae
 126            row[f"smape_{label}"] = scores.smape
 127        rows.append(row)
 128    return pl.DataFrame(rows)
 129
 130
 131def _fixed_form_path(
 132    columns: dict[str, np.ndarray],
 133    truth: np.ndarray,
 134    labels: np.ndarray,
 135    config: Configuration,
 136    library: Library | None = None,
 137    model_features: tuple[str, ...] = MODEL_FEATURES,
 138) -> tuple[dict[int, CrossValidation], Library]:
 139    """Fit once, then cross-validate the resulting forms with only their weights refit.
 140
 141    The fallback used when a caller has no precomputed path to reuse. It exists so that every
 142    cross-validated number in this module comes from the same protocol: one equation, its
 143    terms fixed, its constants recalibrated per fold. Re-running selection inside the folds
 144    would answer a different question and is never done for a reported figure.
 145    """
 146    if library is None:
 147        library = build_library(
 148            DATASET_FEATURES,
 149            model_features,
 150            columns,
 151            max_arity=config.max_arity,
 152            max_abs_zscore=config.max_abs_zscore,
 153        )
 154    result = search(
 155        library,
 156        truth,
 157        max_terms=config.max_terms,
 158        penalty=config.penalty,
 159        pool_size=config.pool_size,
 160        beam_width=config.beam_width,
 161    )
 162    path = cross_validate_fixed_form(library, columns, truth, labels, result.equations, penalty=config.penalty)
 163    return path, library
 164
 165
 166@dataclass(frozen=True)
 167class FittedPath:
 168    """Everything `run_equation` computes before a length is chosen.
 169
 170    Searching and cross-validating every length is the expensive half of an equation; choosing
 171    one length and finalising it is cheap. Keeping the two apart means a length rule only ever
 172    reads a curve that was computed once, and lets `ml_meta_perf.sweep` score hundreds of
 173    configurations without finalising any of them.
 174    """
 175
 176    library: Library
 177    columns: dict[str, np.ndarray]
 178    truth: np.ndarray
 179    datasets: np.ndarray
 180    equations: dict[int, Equation]
 181    paths: dict[str, dict[int, CrossValidation]]
 182    curve: pl.DataFrame
 183    config: Configuration
 184
 185
 186def fit_path(
 187    frame: pl.DataFrame,
 188    dataset_features: tuple[str, ...],
 189    model_features: tuple[str, ...],
 190    config: Configuration,
 191    name: str,
 192) -> FittedPath:
 193    """Search one grammar and cross-validate every length on it, choosing none of them."""
 194    columns = columns_as_arrays(frame, dataset_features + model_features)
 195    truth = target(frame)
 196    datasets = groups(frame, DATASET_COLUMN)
 197    models = groups(frame, MODEL_COLUMN)
 198
 199    library = build_library(
 200        dataset_features,
 201        model_features,
 202        columns,
 203        max_arity=config.max_arity,
 204        max_abs_zscore=config.max_abs_zscore,
 205    )
 206    result = search(
 207        library,
 208        truth,
 209        max_terms=config.max_terms,
 210        penalty=config.penalty,
 211        pool_size=config.pool_size,
 212        beam_width=config.beam_width,
 213        name=name,
 214    )
 215    # Fixed form: the terms are chosen once, here, and only the weights are refit in each
 216    # fold. See `validate.cross_validate_fixed_form` for why that is the reported protocol.
 217    paths = {
 218        label: cross_validate_fixed_form(library, columns, truth, labels, result.equations, penalty=config.penalty)
 219        for label, labels in (("loo_dataset", datasets), ("loo_model", models))
 220    }
 221    # DHO at **every** length, not just the selected one: every selection rule reads the
 222    # complete curve, and E3-MAX selects on this protocol through `selection.floor_curve`.
 223    paths["loo_cell"] = cross_validate_doubly_held_out(
 224        library, columns, truth, datasets, models, result.equations, penalty=config.penalty
 225    )
 226    in_sample = {k: score(truth, eq.predict(columns)) for k, eq in result.equations.items()}
 227    curve = _curve(tuple(sorted(result.equations)), in_sample, paths, truth)
 228    return FittedPath(library, columns, truth, datasets, result.equations, paths, curve, config)
 229
 230
 231def finalise(fitted: FittedPath, size: int) -> EquationReport:
 232    """Publish one length of a fitted path: prune it, and measure how stable its form is."""
 233    if size not in fitted.equations:
 234        raise ValueError(f"selected equation length {size} is unavailable")
 235    config, truth = fitted.config, fitted.truth
 236    # The same folds with selection re-run inside them, kept only for `stability`: the share
 237    # of the equation's terms that survive when one dataset is removed. That is what
 238    # licenses fixing the form, and it is not part of any reported score.
 239    reselected = fold_selections(
 240        fitted.library,
 241        truth,
 242        fitted.datasets,
 243        n_terms=size,
 244        penalty=config.penalty,
 245        pool_size=config.pool_size,
 246        beam_width=config.beam_width,
 247    )
 248    equation = prune(fitted.equations[size], fitted.columns, truth, penalty=config.penalty)
 249
 250    return EquationReport(
 251        equation=equation,
 252        n_terms=size,
 253        arity=config.max_arity,
 254        in_sample=score(truth, equation.predict(fitted.columns)).as_dict(),
 255        curve=fitted.curve,
 256        cross_validated={
 257            label: path[size].scores(truth).as_dict() | path[size].dispersion() for label, path in fitted.paths.items()
 258        },
 259        stability=term_stability(reselected, tuple(term.name for term in equation.terms)),
 260        paths=fitted.paths,
 261    )
 262
 263
 264def _resolve(config: Configuration | None, selection: SelectionConfig | None) -> tuple[Configuration, SelectionConfig]:
 265    """Fill whichever of the two was not given from ``config/study.json``."""
 266    if config is None or selection is None:
 267        study = load_config()
 268        return config or study.search, selection or study.selection
 269    return config, selection
 270
 271
 272def select_length(curve: pl.DataFrame, selection: SelectionConfig) -> int:
 273    """The length rule every published equation shares: `selection.plateau_knee`."""
 274    return plateau_knee(curve, delta=selection.delta, window=selection.window, smoothing=selection.smoothing)
 275
 276
 277def run_equation(
 278    frame: pl.DataFrame,
 279    dataset_features: tuple[str, ...],
 280    model_features: tuple[str, ...],
 281    config: Configuration,
 282    selection: SelectionConfig,
 283    name: str,
 284) -> EquationReport:
 285    """Fit one equation and validate it. The three equations differ **only** in the
 286    features they may draw on, and this is the single code path that says so.
 287
 288    Every one of them is fitted on all 484 rows and scored on all 484 rows, under all four
 289    protocols, and has its length chosen by the same rule. That uniformity is the point: the
 290    gaps between E1, E2 and E3 are only evidence about what each half of the meta-data is
 291    worth if nothing else differs between them -- not the fitting scale, not the protocol,
 292    not the denominator of the R2, not the way the length was picked.
 293
 294    **The length is derived, not asserted**: `select_length` reads it off the complete
 295    cross-validated curve.
 296    """
 297    fitted = fit_path(frame, dataset_features, model_features, config, name)
 298    return finalise(fitted, select_length(fitted.curve, selection))
 299
 300
 301def run_e1(
 302    frame: pl.DataFrame, config: Configuration | None = None, selection: SelectionConfig | None = None
 303) -> EquationReport:
 304    """Dataset features only -- how much of MCC the data alone explains.
 305
 306    Every model on a given dataset shares one feature vector, so this equation can only
 307    ever predict a per-dataset constant. That is not a flaw to be corrected, it is the
 308    control: whatever E3 reaches beyond this is what knowing the model buys.
 309
 310    Least squares finds that per-dataset constant on its own, so aggregating the rows to
 311    twenty dataset means first would be unnecessary as well as harmful to the comparison.
 312    """
 313    config, selection = _resolve(config, selection)
 314    return run_equation(frame, DATASET_FEATURES, (), config, selection, "E1")
 315
 316
 317def run_e2(
 318    frame: pl.DataFrame, config: Configuration | None = None, selection: SelectionConfig | None = None
 319) -> EquationReport:
 320    """Model features only -- the mirror image of E1.
 321
 322    There are six model features and five of them are constant per model, so the library
 323    is tiny and the equation is short by necessity rather than by choice. That is itself
 324    the finding: the meta-data describes datasets far better than it describes models.
 325
 326    Aggregating this one to 25 per-model means -- the mirror of aggregating E1 -- was
 327    measured and is worse. `Processing Units Number` is intentionally allowed to vary with
 328    the dataset because it is a capacity calculation over that dataset's shape. Averaging
 329    it away discards real variation, and the optimum collapses to a single term.
 330    """
 331    config, selection = _resolve(config, selection)
 332    return run_equation(frame, (), MODEL_FEATURES, config, selection, "E2")
 333
 334
 335def run_e3(
 336    frame: pl.DataFrame, config: Configuration | None = None, selection: SelectionConfig | None = None
 337) -> EquationReport:
 338    """Dataset and model features: E3-Valid, the equation the study publishes.
 339
 340    Every dataset feature and every model feature is in the term pool; the configuration
 341    restricts the search only through its hyperparameters.
 342    """
 343    config, selection = _resolve(config, selection)
 344    return run_equation(frame, DATASET_FEATURES, MODEL_FEATURES, config, selection, "E3")
 345
 346
 347def run_capability(
 348    frame: pl.DataFrame, config: Configuration | None = None, selection: SelectionConfig | None = None
 349) -> EquationReport:
 350    """E3-MAX: E3's configuration at the configured capability arity, at the raw maximum of its floor.
 351
 352    **A capability measurement, not a recommendation.** It answers the question E3-Valid cannot
 353    answer about itself -- whether the additive form is out of reach or merely out of the
 354    readable range -- so it takes `selection.floor_argmax` rather than the length rule, with
 355    no smoothing and no complexity penalty. Only the arity differs from E3-Valid.
 356    """
 357    config, selection = _resolve(config, selection)
 358    wider = dataclasses.replace(config, max_arity=selection.capability_arity)
 359    fitted = fit_path(frame, DATASET_FEATURES, MODEL_FEATURES, wider, "E3-MAX")
 360    return finalise(fitted, floor_argmax(fitted.curve))
 361
 362
 363def selection_summary(
 364    e3: EquationReport, capability: EquationReport, config: Configuration, selection: SelectionConfig
 365) -> pl.DataFrame:
 366    """E3-Valid and E3-MAX side by side: the rule each used and where it landed.
 367
 368    ``smoothed_floor`` is the quantity the length rule reads; ``floor`` is the raw worst
 369    protocol at that length; ``at_horizon`` flags a bound that stopped at ``max_terms`` rather
 370    than at a maximum of the data.
 371    """
 372    rows: list[dict[str, object]] = []
 373    for role, report, rule in (
 374        ("E3-Valid", e3, "first plateau knee"),
 375        ("E3-MAX", capability, "raw floor maximum"),
 376    ):
 377        curve = report.curve
 378        position = list(curve["n_terms"]).index(report.n_terms)
 379        rows.append(
 380            {
 381                "role": role,
 382                "rule": rule,
 383                "arity": report.arity,
 384                "n_terms": report.n_terms,
 385                "complexity": complexity(report.arity, report.n_terms),
 386                "r2_in_sample": float(curve["r2_in_sample"][position]),
 387                "r2_loo_dataset": float(curve["r2_loo_dataset"][position]),
 388                "r2_loo_model": float(curve["r2_loo_model"][position]),
 389                "r2_loo_cell": float(curve["r2_loo_cell"][position]),
 390                "floor": float(floor_curve(curve)[position]),
 391                "smoothed_floor": float(smoothed(floor_curve(curve), selection.smoothing)[position]),
 392                "spread": float(protocol_spread(curve)[position]),
 393                "at_horizon": report.n_terms == config.max_terms,
 394            }
 395        )
 396    return pl.DataFrame(rows)
 397
 398
 399def reach_analysis(frame: pl.DataFrame, config: Configuration) -> tuple[pl.DataFrame, pl.DataFrame]:
 400    """What each raw feature is worth, and how far the grammar reaches before any search.
 401
 402    Returns the per-feature table and the three-level ceiling ladder as a one-row frame.
 403    Both are computed from the library alone, so they are available *before* an equation
 404    exists and can be read as an expectation the fitted equation is then held against.
 405    """
 406    columns = columns_as_arrays(frame, DATASET_FEATURES + MODEL_FEATURES)
 407    library = build_library(
 408        DATASET_FEATURES,
 409        MODEL_FEATURES,
 410        columns,
 411        max_arity=config.max_arity,
 412        max_abs_zscore=config.max_abs_zscore,
 413    )
 414    truth = target(frame)
 415    features = DATASET_FEATURES + MODEL_FEATURES
 416    ladder = grammar_ceiling(library, truth, features)
 417    return feature_reach(library, truth, features), pl.DataFrame([ladder])
 418
 419
 420def saturated_analysis(frame: pl.DataFrame, config: Configuration) -> dict[str, float]:
 421    """`analysis.saturated_fit` over the library E3 actually searches.
 422
 423    The library rather than the screened pool, and the E3 grammar rather than the capability
 424    one, because the claim it supports is about the procedure the study publishes: handing
 425    *these* candidates to least squares in one go is what selection is being compared with.
 426    """
 427    columns = columns_as_arrays(frame, DATASET_FEATURES + MODEL_FEATURES)
 428    library = build_library(
 429        DATASET_FEATURES,
 430        MODEL_FEATURES,
 431        columns,
 432        max_arity=config.max_arity,
 433        max_abs_zscore=config.max_abs_zscore,
 434    )
 435    return saturated_fit(library, target(frame), groups(frame, DATASET_COLUMN))
 436
 437
 438def correlation_analysis(frame: pl.DataFrame, config: Configuration, top: int = 15) -> pl.DataFrame:
 439    """Rank candidate terms by how they relate to MCC, linearly and monotonically.
 440
 441    Both correlations are reported because they say different things and the difference
 442    is what guides term design. A term whose Pearson and Spearman correlations agree is
 443    linear in MCC and needs no transform. A term whose Spearman correlation is clearly
 444    the larger is monotone but curved, which is the signal that a log, an inverse or a
 445    ratio will pay for itself. The ``within_pearson`` column is the strict one: both the
 446    term and the target are centred inside each dataset first, so a term is credited
 447    only for variation that dataset identity does not already account for.
 448    """
 449    columns = columns_as_arrays(frame, DATASET_FEATURES + MODEL_FEATURES)
 450    library = build_library(
 451        DATASET_FEATURES,
 452        MODEL_FEATURES,
 453        columns,
 454        max_arity=config.max_arity,
 455        max_abs_zscore=config.max_abs_zscore,
 456    )
 457    table = screen(library, target(frame), groups(frame, DATASET_COLUMN))
 458    return table.with_columns((pl.col("spearman").abs() - pl.col("pearson").abs()).alias("monotone_gap")).head(top)
 459
 460
 461def baselines(frame: pl.DataFrame) -> pl.DataFrame:
 462    """What the equations have to beat, and the ceiling neither of them can pass.
 463
 464    Every trivial predictor is reported at both its mean and its median, because the metrics
 465    disagree about which is the honest opponent. R2 and RMSE are squared-error metrics and the
 466    mean minimises squared error; MAE is minimised by the median, and SMAPE is an
 467    absolute-error ratio that behaves the same way. Quoting the equation's MAE against a
 468    *mean* baseline compares it with a predictor that is not minimising the metric it is being
 469    judged on, which flatters the equation. Both rows are here so each metric can be read
 470    against whichever centre is hardest to beat on it.
 471    """
 472    truth = target(frame)
 473    datasets = groups(frame, DATASET_COLUMN)
 474    models = groups(frame, MODEL_COLUMN)
 475
 476    rows: list[dict[str, object]] = []
 477    for centre in ("mean", "median"):
 478        for label, outer, inner in (
 479            (f"global {centre} (loo-dataset)", datasets, None),
 480            (f"per-model {centre} (loo-dataset)", datasets, models),
 481            (f"global {centre} (loo-model)", models, None),
 482            (f"per-dataset {centre} (loo-model)", models, datasets),
 483        ):
 484            prediction = baseline_group_centre(truth, outer, inner, centre=centre)
 485            rows.append({"baseline": label, **score(truth, prediction).as_dict()})
 486    rows.append(
 487        {
 488            "baseline": "additive mean-based reference (in-sample)",
 489            **score(truth, additive_mean_reference(truth, datasets, models)).as_dict(),
 490        }
 491    )
 492    return pl.DataFrame(rows)
 493
 494
 495def leakage_demonstration(
 496    frame: pl.DataFrame,
 497    n_terms: int,
 498    config: Configuration,
 499    known: dict[str, dict[int, CrossValidation]] | None = None,
 500) -> pl.DataFrame:
 501    """The same equation scored under a random split and under a grouped split.
 502
 503    Dataset features are constant across a dataset's rows, so a random k-fold split puts
 504    the same dataset on both sides of the fold. The equation then recognises the dataset
 505    rather than generalising to it. This table exists so the difference is visible in
 506    numbers rather than asserted in prose.
 507
 508    ``known`` supplies grouped paths that `run_e3` has already computed. They are the same
 509    folds over the same library at the same settings, and beam search records its best
 510    subset at every size as it goes, so a path built to ``max_terms`` contains the shorter
 511    entry this table needs and does not have to be recomputed.
 512    """
 513    columns = columns_as_arrays(frame, DATASET_FEATURES + MODEL_FEATURES)
 514    truth = target(frame)
 515    protocols = {
 516        "random 10-fold (leaky)": ("random", random_kfold_groups(truth.shape[0])),
 517        "leave-one-dataset-out": ("loo_dataset", groups(frame, DATASET_COLUMN)),
 518        "leave-one-model-out": ("loo_model", groups(frame, MODEL_COLUMN)),
 519    }
 520    library: Library | None = None
 521    rows: list[dict[str, object]] = []
 522    for label, (key, labels) in protocols.items():
 523        path = (known or {}).get(key)
 524        if path is None or n_terms not in path:
 525            path, library = _fixed_form_path(columns, truth, labels, config, library)
 526        size = n_terms if n_terms in path else max(path)
 527        rows.append({"protocol": label, **path[size].scores(truth).as_dict()})
 528    return pl.DataFrame(rows)
 529
 530
 531def selection_optimism(frame: pl.DataFrame, e3: EquationReport, config: Configuration) -> pl.DataFrame:
 532    """LODO with the terms **re-chosen inside every fold**, beside the reported fixed-form LODO.
 533
 534    The reported protocol fixes the form -- chosen once, on all rows -- and refits only the
 535    weights per fold (`validate.cross_validate_fixed_form`). That is a position about what the
 536    equation is, and it has a price: the form was chosen with the held-out dataset in view. This
 537    measures the price, by repeating the selection itself without the held-out dataset, at
 538    E3-Valid's length and configuration, and scoring those predictions.
 539
 540    **A diagnostic, not a score.** Twenty folds fit twenty different equations, so the nested
 541    row describes the discovery procedure rather than the published equation; it is reported
 542    so the fixed-form transfer numbers are read as the upper estimate they are.
 543    """
 544    columns = columns_as_arrays(frame, DATASET_FEATURES + MODEL_FEATURES)
 545    truth, datasets = target(frame), groups(frame, DATASET_COLUMN)
 546    library = build_library(
 547        DATASET_FEATURES, MODEL_FEATURES, columns, max_arity=config.max_arity, max_abs_zscore=config.max_abs_zscore
 548    )
 549    position = {name: index for index, name in enumerate(library.names)}
 550    nested = np.zeros_like(truth)
 551    for _, train, test in leave_one_group_out(datasets):
 552        result = search(
 553            library_rows(library, train),
 554            truth[train],
 555            max_terms=e3.n_terms,
 556            penalty=config.penalty,
 557            pool_size=config.pool_size,
 558            beam_width=config.beam_width,
 559        )
 560        equation = result.equations[max(size for size in result.equations if size <= e3.n_terms)]
 561        held = library.matrix[test][:, [position[term.name] for term in equation.terms]]
 562        nested[test] = np.clip(
 563            equation.intercept + held @ np.asarray(equation.weights), truth[train].min(), truth[train].max()
 564        )
 565    fixed = e3.paths["loo_dataset"][e3.n_terms].predictions
 566    return pl.DataFrame(
 567        [
 568            {"LODO": "form chosen once, on all rows (reported)", **score(truth, fixed).as_dict()},
 569            {"LODO": "form re-chosen inside every fold", **score(truth, nested).as_dict()},
 570        ]
 571    )
 572
 573
 574def identity_ceiling(frame: pl.DataFrame, e3: EquationReport) -> pl.DataFrame:
 575    """The upper bound on what any model descriptor could add, measured rather than argued.
 576
 577    Under LODO **every one of the 25 models appears in every training fold**,
 578    so the equation's residual can be averaged per model on the training rows and applied to
 579    the held-out dataset with no leak. That replaces the model descriptors with the best
 580    possible substitute -- the model's *identity*, fitted freely -- and what it adds is
 581    therefore a ceiling on what any descriptor set could reach by telling these classifiers
 582    apart. It is the number the study's central negative result is measured against: every
 583    rejected model-side encoding is rejected for failing to recover it.
 584
 585    Two rungs, because the ceiling depends on what the identity is allowed to carry. A level
 586    per model shifts every one of its rows equally; a level and a slope on a dataset feature
 587    lets a model's advantage depend on the data, which is what a *capability* descriptor would
 588    have to do. `identity.correct_out_of_fold` reuses the per-fold equations the reported path
 589    already recorded, so the protocol is unchanged and no search is repeated.
 590
 591    Generated because the chapter's hand-written copy of this table had all three rows wrong
 592    and listed its two rungs as identical, which no run of this function can produce.
 593    """
 594    columns = columns_as_arrays(frame, DATASET_FEATURES + MODEL_FEATURES)
 595    truth = target(frame)
 596    datasets, models = groups(frame, DATASET_COLUMN), groups(frame, MODEL_COLUMN)
 597    path = e3.paths.get("loo_dataset")
 598    size = len(e3.equation.terms)
 599    if path is None or size not in path:
 600        return pl.DataFrame(schema={"correction": pl.String, "r2_loo_dataset": pl.Float64, "mae": pl.Float64})
 601    fold = path[size]
 602    rungs = {
 603        f"none (E3-Valid, {size} terms)": fold.predictions,
 604        "per-model level": correct_out_of_fold(fold, columns, truth, datasets, models),
 605        "per-model level and slope": correct_out_of_fold(
 606            fold, columns, truth, datasets, models, features=DATASET_FEATURES
 607        ),
 608    }
 609
 610    def per_dataset_mae(prediction: np.ndarray) -> np.ndarray:
 611        return np.array([mae(truth[test], prediction[test]) for _, _, test in leave_one_group_out(datasets)])
 612
 613    # **Whether the gap is real is a paired question, not a difference of two pooled numbers.**
 614    # That distinction has produced three wrong conclusions on this project, and this gap is
 615    # the one the study's central negative result is measured against -- so each rung is paired
 616    # against the uncorrected equation over the twenty held-out datasets, on per-dataset MAE.
 617    baseline = per_dataset_mae(fold.predictions)
 618    rows: list[dict[str, object]] = []
 619    for label, values in rungs.items():
 620        row: dict[str, object] = {
 621            "correction": label,
 622            "r2_loo_dataset": r2_score(truth, values),
 623            "mae": mae(truth, values),
 624        }
 625        if label.startswith("none"):
 626            row |= {
 627                "gain": 0.0,
 628                "ci_low": float("nan"),
 629                "ci_high": float("nan"),
 630                "sign_p": float("nan"),
 631                "wins": 0,
 632                "verdict": "baseline",
 633            }
 634        else:
 635            paired = paired_comparison(per_dataset_mae(values), baseline, lower_is_better=True)
 636            row |= {
 637                "gain": paired.mean,
 638                "ci_low": paired.low,
 639                "ci_high": paired.high,
 640                "sign_p": paired.p_value,
 641                "wins": paired.wins,
 642                "verdict": "real" if paired.significant else "tie",
 643            }
 644        rows.append(row)
 645    return pl.DataFrame(rows)
 646
 647
 648def interaction_reached(frame: pl.DataFrame, e3: EquationReport) -> pl.DataFrame:
 649    """How much of the leading interaction pattern the fitted equation actually reaches.
 650
 651    The oracle ladder prices interaction components; this says whether E3 gets any of
 652    them. Reported under both the fit and the held-out protocol, because an equation can
 653    align with a pattern under IS and lose it out of fold -- and that difference is the
 654    whole question for a component nobody can predict from meta-features.
 655    """
 656    columns = columns_as_arrays(frame, DATASET_FEATURES + MODEL_FEATURES)
 657    truth = target(frame)
 658    datasets, models = groups(frame, DATASET_COLUMN), groups(frame, MODEL_COLUMN)
 659    size = len(e3.equation.terms)
 660    predictions = {"in-sample": e3.equation.predict(columns)}
 661    path = e3.paths.get("loo_dataset")
 662    if path is not None and size in path:
 663        predictions["leave-one-dataset-out"] = path[size].predictions
 664    rows = [
 665        {"protocol": label, "rank": rank, **interaction_capture(truth, prediction, datasets, models, rank=rank)}
 666        for label, prediction in predictions.items()
 667        for rank in (1, 2)
 668    ]
 669    return pl.DataFrame(rows)
 670
 671
 672def comparison(
 673    frame: pl.DataFrame,
 674    e1: EquationReport,
 675    e3: EquationReport,
 676    e2: EquationReport | None = None,
 677    e3_capability: EquationReport | None = None,
 678) -> pl.DataFrame:
 679    """E1 against E3 on the one scale where they are comparable: all rows.
 680
 681    All three equations are already fitted and scored on the same 484 rows, so this table
 682    is not correcting a denominator -- an earlier version of E1 was fitted on 20 aggregated
 683    dataset means and it was. What it adds is the two *ceilings* beside the equations: the
 684    true per-dataset and per-model means, which are the most any equation restricted to
 685    that half of the meta-data could reach. Without them a reader can compare the raw E1 and
 686    E3 scores and conclude the dataset side is weak, even when E1 is already close to its
 687    own structural ceiling.
 688    """
 689    columns = columns_as_arrays(frame, DATASET_FEATURES + MODEL_FEATURES)
 690    truth = target(frame)
 691    datasets = groups(frame, DATASET_COLUMN)
 692    models = groups(frame, MODEL_COLUMN)
 693
 694    dataset_ceiling = np.zeros_like(truth)
 695    for label in np.unique(datasets):
 696        mask = datasets == label
 697        dataset_ceiling[mask] = truth[mask].mean()
 698    model_ceiling = np.zeros_like(truth)
 699    for label in np.unique(models):
 700        mask = models == label
 701        model_ceiling[mask] = truth[mask].mean()
 702    if e2 is None:
 703        e2 = run_e2(frame)
 704
 705    return pl.DataFrame(
 706        [
 707            {
 708                "equation": f"E1, dataset only ({e1.equation.n_terms} terms)",
 709                "n_terms": e1.equation.n_terms,
 710                **score(truth, e1.equation.predict(columns)).as_dict(),
 711            },
 712            {
 713                "equation": "E1 reference: true dataset means",
 714                "n_terms": None,
 715                **score(truth, dataset_ceiling).as_dict(),
 716            },
 717            {
 718                "equation": f"E2, model only ({e2.equation.n_terms} terms)",
 719                "n_terms": e2.equation.n_terms,
 720                **score(truth, e2.equation.predict(columns)).as_dict(),
 721            },
 722            {
 723                "equation": "E2 reference: true model means",
 724                "n_terms": None,
 725                **score(truth, model_ceiling).as_dict(),
 726            },
 727            {
 728                "equation": f"E3-Valid, dataset + model ({e3.equation.n_terms} terms)",
 729                "n_terms": e3.equation.n_terms,
 730                **score(truth, e3.equation.predict(columns)).as_dict(),
 731            },
 732            *(
 733                [
 734                    {
 735                        "equation": (
 736                            f"E3-MAX, arity "
 737                            f"{max((len(term.features) for term in e3_capability.equation.terms), default=0)}"
 738                            f" ({len(e3_capability.equation.terms)} terms)"
 739                        ),
 740                        "n_terms": len(e3_capability.equation.terms),
 741                        **score(truth, e3_capability.equation.predict(columns)).as_dict(),
 742                    }
 743                ]
 744                if e3_capability is not None and e3_capability is not e3
 745                else []
 746            ),
 747            {
 748                "equation": "reference: additive mean-based reference",
 749                "n_terms": None,
 750                **score(truth, additive_mean_reference(truth, datasets, models)).as_dict(),
 751            },
 752        ]
 753    )
 754
 755
 756def decision_quality(
 757    frame: pl.DataFrame,
 758    n_terms: int,
 759    config: Configuration,
 760    known: dict[int, CrossValidation] | None = None,
 761) -> pl.DataFrame:
 762    """Go/no-go decision quality, with both the dataset and the model held out.
 763
 764    Same requirement as `model_selection` and for the same reason: the question is asked about
 765    a pair nobody has run, so neither half of it may be in the training set. ``known`` is
 766    `run_e3`'s LODO path and is the fallback when the DHO fit
 767    cannot be built.
 768    """
 769    truth = target(frame)
 770    datasets = groups(frame, DATASET_COLUMN)
 771    doubly = doubly_held_out_predictions(frame, n_terms, config)
 772    if doubly is not None:
 773        return decision_report(truth, doubly, datasets)
 774    columns = columns_as_arrays(frame, DATASET_FEATURES + MODEL_FEATURES)
 775    path = known
 776    if path is None or n_terms not in path:
 777        path, _ = _fixed_form_path(columns, truth, datasets, config)
 778    size = n_terms if n_terms in path else max(path)
 779    return decision_report(truth, path[size].predictions, datasets)
 780
 781
 782def length_comparison(frame: pl.DataFrame, e3: EquationReport, config: Configuration) -> pl.DataFrame:
 783    """Every equation length paired against the E3-Valid selection, per held-out dataset.
 784
 785    E3-Valid is selected before this diagnostic is computed. The table measures how each
 786    length differs from that retained equation; it does not define another selection rule.
 787
 788    Per-dataset **MAE** rather than per-dataset R2, because three datasets here have almost no
 789    within-dataset variance -- `5G_Slicing`'s 25 models all score about 0.986 -- and an R2
 790    over a near-constant target is dominated by its denominator. See
 791    `validate.CrossValidation.dispersion`.
 792
 793    A row marked ``worse`` differs significantly from E3-Valid under this paired comparison.
 794    """
 795    truth = target(frame)
 796    datasets = groups(frame, DATASET_COLUMN)
 797    labels = np.unique(datasets)
 798    path = e3.paths.get("loo_dataset", {})
 799    if not path:
 800        return pl.DataFrame()
 801    # Referenced to what the rule chose, not to a length passed in: the table has to be able
 802    # to say that the rule's own pick is the right one, which it cannot do if the pick is the
 803    # thing being assumed.
 804    chosen = e3.n_terms
 805    if chosen not in path:
 806        chosen = max(path)
 807    if chosen not in path:
 808        return pl.DataFrame()
 809    published = chosen
 810
 811    def per_fold(size: int) -> np.ndarray:
 812        prediction = path[size].predictions
 813        return np.array([mae(truth[datasets == label], prediction[datasets == label]) for label in labels])
 814
 815    reference = per_fold(published)
 816    rows: list[dict[str, object]] = []
 817    for size in sorted(path):
 818        result = paired_comparison(reference, per_fold(size), lower_is_better=True)
 819        # `mean` is positive when the first argument -- the published equation -- is better.
 820        verdict = (
 821            "selected"
 822            if size == published
 823            else ("worse" if result.significant and result.mean > 0 else ("better" if result.significant else "tie"))
 824        )
 825        rows.append(
 826            {
 827                "n_terms": size,
 828                "r2_loo_dataset": path[size].scores(truth).r2,
 829                "mae_loo_dataset": float(per_fold(size).mean()),
 830                "mean_difference": result.mean,
 831                "p_value": result.p_value,
 832                "ci_low": result.low,
 833                "ci_high": result.high,
 834                "verdict": verdict,
 835            }
 836        )
 837    return pl.DataFrame(rows)
 838
 839
 840def doubly_held_out_predictions(frame: pl.DataFrame, n_terms: int, config: Configuration) -> np.ndarray | None:
 841    """E3's predictions with **both** the dataset and the model of each cell held out.
 842
 843    Used for the ranking and the threshold decision and for nothing else. Those two are the
 844    questions a practitioner actually asks -- *which model should I run on this data*, and
 845    *will it clear my bar* -- and both are asked about a pair that has not been run. Neither
 846    single-group protocol answers that: LODO has seen the learner on nineteen
 847    other problems, LOMO has seen the dataset.
 848
 849    The R2 curve and the length rule stay on the single-group protocols, which is the right
 850    scope for them: they are about how the equation degrades as one axis becomes unfamiliar,
 851    and a per-cell refit would answer a different question at 32 times the cost.
 852    """
 853    columns = columns_as_arrays(frame, DATASET_FEATURES + MODEL_FEATURES)
 854    truth = target(frame)
 855    datasets = groups(frame, DATASET_COLUMN)
 856    models = groups(frame, MODEL_COLUMN)
 857    library = build_library(
 858        DATASET_FEATURES,
 859        MODEL_FEATURES,
 860        columns,
 861        max_arity=config.max_arity,
 862        max_abs_zscore=config.max_abs_zscore,
 863    )
 864    result = search(
 865        library,
 866        truth,
 867        max_terms=config.max_terms,
 868        penalty=config.penalty,
 869        pool_size=config.pool_size,
 870        beam_width=config.beam_width,
 871    )
 872    if n_terms not in result.equations:
 873        return None
 874    size = n_terms
 875    path = cross_validate_doubly_held_out(
 876        library,
 877        columns,
 878        truth,
 879        datasets,
 880        models,
 881        {size: result.equations[size]},
 882        penalty=config.penalty,
 883    )
 884    return path[size].predictions if size in path else None
 885
 886
 887def _e3_predictions(frame: pl.DataFrame, e3: EquationReport, protocol: str, config: Configuration) -> np.ndarray | None:
 888    """E3's predictions under one protocol.
 889
 890    ``in_sample``, ``loo_dataset``, ``loo_model``, or ``loo_cell`` -- the last holding out both
 891    the dataset and the model of every cell, which is the protocol the ranking and the
 892    threshold decision are reported under.
 893    """
 894    if protocol == "loo_cell":
 895        return doubly_held_out_predictions(frame, e3.n_terms, config)
 896    if protocol == "in_sample":
 897        columns = columns_as_arrays(frame, DATASET_FEATURES + MODEL_FEATURES)
 898        return e3.equation.predict(columns)
 899    path = e3.paths.get(protocol)
 900    if not path:
 901        return None
 902    size = e3.n_terms if e3.n_terms in path else max(path)
 903    return path[size].predictions if size in path else None
 904
 905
 906def model_selection(
 907    frame: pl.DataFrame,
 908    e3: EquationReport,
 909    config: Configuration,
 910    protocol: str = "loo_cell",
 911) -> pl.DataFrame:
 912    """Can the equation pick a good model for a dataset it has never run it on?
 913
 914    **Reported under ``loo_cell``: both the dataset and the model of every cell are out of the
 915    training set.** Anything weaker answers a different question. LODO has
 916    seen the learner on nineteen other problems; LOMO has seen the dataset. A
 917    recommendation is asked about a pair that has not been run, so both have to go.
 918
 919    The cost of that is small, and structurally so rather than by luck: a ranking depends only
 920    on the *within-dataset* ordering of models, and the dataset-side terms shift every model
 921    on a dataset by the same amount, so holding the dataset out moves the level and not the
 922    order. `ranking_baselines` reports every protocol side by side so that stays visible.
 923    """
 924    truth = target(frame)
 925    datasets = groups(frame, DATASET_COLUMN)
 926    prediction = _e3_predictions(frame, e3, protocol, config)
 927    if prediction is None:
 928        return pl.DataFrame()
 929    return ranking_report(truth, prediction, datasets)
 930
 931
 932#: How an opaque estimator's protocol keys are labelled in the comparison tables. Every row
 933#: names its protocol, because a comparison is only a comparison if both sides were scored
 934#: under the same one -- which this project has shipped wrong twice.
 935OPAQUE_PROTOCOLS = {
 936    "in_sample": "in-sample",
 937    "loo_dataset": "loo-dataset",
 938    "loo_model": "loo-model",
 939    "loo_cell": "loo-cell: both held out",
 940}
 941
 942
 943def _opaque_candidates(opaque: OpaqueRun) -> dict[str, np.ndarray]:
 944    """Every opaque estimator under every protocol, keyed the way the equation's rows are."""
 945    return {
 946        f"{label} ({OPAQUE_PROTOCOLS[protocol]})": values
 947        for label, held in opaque.predictions.items()
 948        for protocol, values in held.items()
 949        if protocol in OPAQUE_PROTOCOLS
 950    }
 951
 952
 953def ranking_baselines(
 954    frame: pl.DataFrame,
 955    e3: EquationReport,
 956    config: Configuration,
 957    opaque: OpaqueRun | None = None,
 958) -> pl.DataFrame:
 959    """The equation's ranking against the trivial rankings, averaged over datasets.
 960
 961    **Every predictor is scored under the same protocol**, which is what makes the comparison
 962    a comparison. The equation appears under IS, LODO, LOMO, and DHO, so the cost of
 963    generalisation on this task is visible rather than assumed. The trivial predictors appear
 964    under LODO, which is the only way a per-model
 965    centre can be computed honestly.
 966
 967    Both centres are reported for the same reason the error metrics report both: the per-model
 968    mean and the per-model median are different orderings, and which is harder to beat is a
 969    question rather than an assumption.
 970    """
 971    truth = target(frame)
 972    datasets = groups(frame, DATASET_COLUMN)
 973    models = groups(frame, MODEL_COLUMN)
 974
 975    candidates: dict[str, np.ndarray] = {}
 976    for label, protocol in (
 977        ("equation (in-sample)", "in_sample"),
 978        ("equation (loo-dataset)", "loo_dataset"),
 979        ("equation (loo-model)", "loo_model"),
 980    ):
 981        prediction = _e3_predictions(frame, e3, protocol, config)
 982        if prediction is not None:
 983            candidates[label] = prediction
 984    doubly = doubly_held_out_predictions(frame, e3.n_terms, config)
 985    if doubly is not None:
 986        candidates["equation (loo-cell: both held out)"] = doubly
 987    # The trivial predictors cannot be computed under DHO (internally ``loo_cell``) at all: a model
 988    # held out of every fold has no rows to average, so "how well does this model usually do"
 989    # has no value. They appear under LODO, which is the only way they exist -- and
 990    # that hands them the model identity denied to the DHO row of the equation.
 991    candidates["per-model mean (loo-dataset)"] = baseline_group_centre(truth, datasets, models, centre="mean")
 992    candidates["per-model median (loo-dataset)"] = baseline_group_centre(truth, datasets, models, centre="median")
 993    # The opaque opponent under every protocol it has, so each of its rows can be read against
 994    # the equation row that was allowed to see the same things. The DHO rows are the
 995    # comparison that matters: there neither side has the dataset or the model, and the trivial
 996    # baselines cannot be computed at all.
 997    if opaque is not None:
 998        candidates |= _opaque_candidates(opaque)
 999
1000    tables = {
1001        label: ranking_report(truth, prediction, datasets).sort("group") for label, prediction in candidates.items()
1002    }
1003    reference = tables["equation (loo-dataset)"] if "equation (loo-dataset)" in tables else next(iter(tables.values()))
1004
1005    rows: list[dict[str, object]] = []
1006    for label, table in tables.items():
1007        row: dict[str, object] = {
1008            "predictor": label,
1009            **{
1010                column: float(table[column].to_numpy().mean())
1011                for column in ("ap", "mrr", "hit_at_1", "regret")
1012                if column in table.columns
1013            },
1014            "datasets": table.height,
1015        }
1016        # The mean over twenty datasets is not the comparison. hit@1 moves only in steps of
1017        # 0.05 on twenty folds, so a single dataset flipping its top pick shifts it by more
1018        # than the gaps being read. The paired test against the reported equation is what
1019        # licenses any claim from this table.
1020        if label != "equation (loo-dataset)" and "ap" in table.columns:
1021            paired = paired_comparison(reference["ap"].to_numpy(), table["ap"].to_numpy())
1022            row["ap_vs_e3_p"] = paired.p_value
1023            row["ap_vs_e3_significant"] = paired.significant
1024        rows.append(row)
1025    return pl.DataFrame(rows)
1026
1027
1028def decision_baselines(
1029    frame: pl.DataFrame,
1030    config: Configuration,
1031    known: dict[int, CrossValidation] | None = None,
1032    e3: EquationReport | None = None,
1033    opaque: OpaqueRun | None = None,
1034) -> pl.DataFrame:
1035    """The above-or-below-threshold decision, for the equation and both trivial centres.
1036
1037    Every predictor under the same protocol, as in `ranking_baselines`, and the equation under
1038    all four so the cost of generalisation on the decision is measured rather than assumed.
1039    `decision_report` already carries a majority-class column, which is the floor any rule has
1040    to clear; these are the harder comparison.
1041    """
1042    columns = columns_as_arrays(frame, DATASET_FEATURES + MODEL_FEATURES)
1043    truth = target(frame)
1044    datasets = groups(frame, DATASET_COLUMN)
1045    models = groups(frame, MODEL_COLUMN)
1046
1047    candidates: dict[str, np.ndarray] = {}
1048    if e3 is not None:
1049        for label, protocol in (
1050            ("equation (in-sample)", "in_sample"),
1051            ("equation (loo-dataset)", "loo_dataset"),
1052            ("equation (loo-model)", "loo_model"),
1053        ):
1054            prediction = _e3_predictions(frame, e3, protocol, config)
1055            if prediction is not None:
1056                candidates[label] = prediction
1057        doubly = doubly_held_out_predictions(frame, e3.n_terms, config)
1058        if doubly is not None:
1059            candidates["equation (loo-cell: both held out)"] = doubly
1060    else:
1061        path = known
1062        if path is None or (e3 is not None and e3.n_terms not in path):
1063            path, _ = _fixed_form_path(columns, truth, datasets, config)
1064        candidates["equation (loo-dataset)"] = path[e3.n_terms if e3 is not None else max(path)].predictions
1065    candidates["per-model mean (loo-dataset)"] = baseline_group_centre(truth, datasets, models, centre="mean")
1066    candidates["per-model median (loo-dataset)"] = baseline_group_centre(truth, datasets, models, centre="median")
1067    if opaque is not None:
1068        candidates |= _opaque_candidates(opaque)
1069
1070    rows: list[dict[str, object]] = []
1071    for label, prediction in candidates.items():
1072        for row in decision_report(truth, prediction, datasets).to_dicts():
1073            rows.append({"predictor": label, **row})
1074    return pl.DataFrame(rows)
1075
1076
1077@dataclass
1078class Report:
1079    """Everything the experiment produces."""
1080
1081    e1: EquationReport
1082    e3: EquationReport
1083    e2: EquationReport
1084    #: E3-MAX: the same configuration at the configured capability arity at the raw floor maximum.
1085    #: Reported to show how far the additive form reaches, not as the study's
1086    #: recommendation. See `run_capability`.
1087    e3_capability: EquationReport
1088    #: E3-Valid and E3-MAX side by side, with the rule each used. See `selection_summary`.
1089    grammars: pl.DataFrame
1090    #: Feature-level effects and stability before evidence thresholds are applied.
1091    feature_practices: pl.DataFrame
1092    practices: pl.DataFrame
1093    effects: pl.DataFrame
1094    shares: pl.DataFrame
1095    decomposition: pl.DataFrame
1096    correlations: pl.DataFrame
1097    #: What each raw feature carries on its own, and the three-level ceiling the grammar
1098    #: implies before any search runs. See `analysis.grammar_ceiling`.
1099    reach: pl.DataFrame
1100    ceiling: pl.DataFrame
1101    baselines: pl.DataFrame
1102    comparison: pl.DataFrame
1103    leakage: pl.DataFrame
1104    selection: pl.DataFrame
1105    #: Every length paired against E3-Valid as a sensitivity analysis.
1106    length_choice: pl.DataFrame
1107    oracles: pl.DataFrame
1108    #: How far the equation's own interactions lie along the leading component the
1109    #: oracle ladder finds. Without it the ladder measures a ceiling and says nothing
1110    #: about whether the equation reaches any of it.
1111    interaction: pl.DataFrame
1112    decision: pl.DataFrame
1113    #: The ranking and threshold decisions against the trivial predictors, at both
1114    #: centres. An evaluation without them says how well the equation does and not
1115    #: whether it beats ordering the models by how well they usually do.
1116    ranking_baselines: pl.DataFrame
1117    decision_baselines: pl.DataFrame
1118    #: Standard opaque regressors on the same rows under the same protocols. The priced
1119    #: other side of the study's trade: see `ml_meta_perf.opaque`.
1120    opaque: pl.DataFrame
1121    #: What per-model *identity* adds to the reported path, which is the ceiling on what any
1122    #: model descriptor could reach. See `identity_ceiling`.
1123    identity: pl.DataFrame
1124    #: What an unpenalised least-squares fit over the *whole* library reaches under IS and
1125    #: held out. The control for the selection stage: if this transferred, the beam search
1126    #: and the length rule would be machinery in search of a problem. See
1127    #: `analysis.saturated_fit`.
1128    saturated: pl.DataFrame
1129    #: LODO with E3-Valid's terms re-chosen in every fold, beside the reported fixed form: how
1130    #: much of the transfer number belongs to choosing the form on all rows. See
1131    #: `selection_optimism`.
1132    optimism: pl.DataFrame
1133
1134
1135def run(
1136    path: str | None = None,
1137    *,
1138    config: Configuration | None = None,
1139    selection: SelectionConfig | None = None,
1140    opaque: OpaqueConfig | None = None,
1141    opaque_models: tuple[tuple[str, OpaqueBuilder], ...] | None = None,
1142) -> Report:
1143    """Run the whole study.
1144
1145    ``config``, ``selection`` and ``opaque`` default to ``config/study.json``. **One
1146    configuration serves all three equations**: every table that mentions an equation is
1147    computed from the same object, so an override changes the study consistently.
1148
1149    ``opaque_models`` replaces the estimators the comparison is run against, in the shape
1150    `opaque.evaluate` takes. The opaque side is scikit-learn's and the DHO refit around it is
1151    this project's, so a caller checking the wiring can exercise the second without paying
1152    for the first -- which at 484 refits per estimator is most of what a run costs.
1153    """
1154    frame = load(path)
1155    config, selection = _resolve(config, selection)
1156    e1 = run_e1(frame, config, selection)
1157    e2 = run_e2(frame, config, selection)
1158    e3 = run_e3(frame, config, selection)
1159    e3_capability = run_capability(frame, config, selection)
1160    config_e3 = config
1161    columns = columns_as_arrays(frame, DATASET_FEATURES + MODEL_FEATURES)
1162    reach, ceiling = reach_analysis(frame, config_e3)
1163    # Fitted once and shared: the folds are the expensive part, and the regression table and
1164    # the two decision comparisons have to be scored from the same predictions or they can
1165    # disagree with each other.
1166    if opaque_models is None:
1167        opaque_models = opaque_estimators(opaque or load_config().opaque)
1168    opaque_run = opaque_evaluate(frame, models=opaque_models)
1169    return Report(
1170        opaque=opaque_run.table,
1171        saturated=pl.DataFrame([saturated_analysis(frame, config_e3)]),
1172        optimism=selection_optimism(frame, e3, config_e3),
1173        e1=e1,
1174        e2=e2,
1175        e3=e3,
1176        e3_capability=e3_capability,
1177        grammars=selection_summary(e3, e3_capability, config, selection),
1178        feature_practices=feature_practices(e3.equation, columns, e3.stability),
1179        practices=best_practices(e3.equation, columns, e3.stability),
1180        effects=term_effects(e3.equation, columns, DATASET_FEATURES, MODEL_FEATURES),
1181        shares=group_shares(e3.equation, columns, DATASET_FEATURES, MODEL_FEATURES),
1182        decomposition=variance_decomposition(target(frame), groups(frame, DATASET_COLUMN), groups(frame, MODEL_COLUMN)),
1183        correlations=correlation_analysis(frame, config_e3),
1184        reach=reach,
1185        ceiling=ceiling,
1186        baselines=baselines(frame),
1187        comparison=comparison(frame, e1, e3, e2, e3_capability),
1188        leakage=leakage_demonstration(frame, e3.n_terms, config_e3, e3.paths),
1189        selection=model_selection(frame, e3, config_e3, "loo_cell"),
1190        decision=decision_quality(frame, e3.n_terms, config_e3, e3.paths.get("loo_dataset")),
1191        ranking_baselines=ranking_baselines(frame, e3, config_e3, opaque_run),
1192        decision_baselines=decision_baselines(frame, config_e3, e3.paths.get("loo_dataset"), e3, opaque_run),
1193        length_choice=length_comparison(frame, e3, config_e3),
1194        oracles=oracle_ladder(target(frame), groups(frame, DATASET_COLUMN), groups(frame, MODEL_COLUMN)),
1195        interaction=interaction_reached(frame, e3),
1196        identity=identity_ceiling(frame, e3),
1197    )