ml_meta_perf.opaque

  1"""Opaque regressors on the same rows, scored under the same protocols.
  2
  3The study's central claim is a trade: accuracy given up for an equation that can be read.
  4A trade needs both sides priced, and the expensive side is what a flexible model reaches on
  5this meta-data. That number is quoted throughout the chapters and was, for a long time,
  6measured once by hand and reproducible by nobody -- which for a study whose whole argument is
  7that its analysis is generated rather than narrated was the wrong way round. This module
  8computes it on every run.
  9
 10The evaluation protocols are in-sample (IS), leave-one-dataset-out (LODO),
 11leave-one-model-out (LOMO), and doubly held out (DHO).
 12
 13Three regressors, chosen because they are what a reader would reach for rather than because
 14they are the strongest available:
 15
 16* **Ridge** with a cross-validated penalty -- the linear-but-opaque control. It sees the same
 17  eighteen raw columns the grammar is built over, unstandardised structure and all.
 18* **Random forest** -- the model that produces the R2 ~ 0.9 figures reported for opaque
 19  meta-learners, and the one whose transfer collapse is the point.
 20* **Gradient boosting** -- the same family with a different bias, so the forest result cannot
 21  be read as an artefact of one estimator.
 22
 23**Every one is scored under the study's own protocols**, by the same
 24`ml_meta_perf.validate.leave_one_group_out` splits, with the same clip to the training fold's
 25observed range that every reported number uses. A comparison is only a comparison if both
 26sides are scored the same way, and this went wrong twice in this project's history in ways
 27that flattered the equation.
 28
 29Fixed seeds, and no tuning. Tuning the opaque side would strengthen it and change nothing
 30about the conclusion -- the failure is that 20 dataset groups let a flexible model memorise
 31group identity, which more trees cannot fix -- but an untuned opponent should be described as
 32one, so `n_estimators` and depth are stated in the table rather than buried.
 33
 34Study chapter: [5. Evaluation][study-chapter] -- the rationale, in
 35prose, with the figures.
 36
 37[study-chapter]: https://github.com/mariolpantunes/ml-meta-perf/blob/main/assets/docs/05-evaluation.md
 38"""
 39
 40from __future__ import annotations
 41
 42from collections.abc import Callable
 43from dataclasses import dataclass
 44from typing import Protocol
 45
 46import numpy as np
 47import polars as pl
 48from sklearn.ensemble import GradientBoostingRegressor, RandomForestRegressor
 49from sklearn.linear_model import RidgeCV
 50from sklearn.utils.parallel import Parallel, delayed
 51
 52from ml_meta_perf.config import OpaqueConfig, load_config
 53from ml_meta_perf.data import (
 54    ALL_FEATURES,
 55    DATASET_COLUMN,
 56    MODEL_COLUMN,
 57    columns_as_arrays,
 58    groups,
 59    target,
 60)
 61from ml_meta_perf.model import MCC_LOWER, MCC_UPPER
 62from ml_meta_perf.stats import mae, r2_score
 63from ml_meta_perf.validate import leave_one_group_out
 64
 65#: Seed for every estimator that takes one. Fixed so the comparison is reproducible; the
 66#: conclusion does not depend on it: the gap being measured is two orders of magnitude
 67#: larger than the spread a different seed produces.
 68SEED = 0
 69
 70#: The ensemble sizes are hyperparameters and live in ``config/study.json`` (`config.OpaqueConfig`).
 71#: They were chosen on a grid scored under all four protocols: a third of the features per split
 72#: transfers far better than all of them, and a hundred trees reach within noise of three
 73#: hundred at a third of the cost. Neither conclusion the study draws depends on the size: the
 74#: forest's argument is that a flexible model fits this meta-data almost perfectly and reaches
 75#: nothing on an unseen cell, which is a property of the design matrix -- twenty dataset groups,
 76#: features constant within a group -- rather than of how many trees vote on it.
 77
 78
 79class Regressor(Protocol):
 80    """The whole of what this module needs from an estimator: fit it, then predict with it.
 81
 82    Stated as a protocol rather than a union of the three scikit-learn classes, because that
 83    is the actual contract. Anything satisfying it can be handed to `evaluate`, which is what
 84    lets the tests exercise the protocols without paying for real ensembles.
 85    """
 86
 87    # Positional-only: scikit-learn names these ``X`` and ``y``, and a protocol that insisted
 88    # on its own names would match nothing in the library it exists to describe.
 89    def fit(self, design: np.ndarray, truth: np.ndarray, /) -> object: ...
 90
 91    def predict(self, design: np.ndarray, /) -> object: ...
 92
 93
 94#: A factory taking the thread budget for one fit and returning a fresh unfitted regressor.
 95#:
 96#: **A callable, not a name to dispatch on.** The protocols below are the study's contribution
 97#: on this side -- the DHO refit especially -- and the estimators are scikit-learn's
 98#: contribution, which needs no testing here. Injecting the factory separates them: a caller
 99#: that wants to check the *plumbing* passes something trivial and pays nothing, and the tests
100#: do exactly that. It also makes these factories what this module always said they were.
101Builder = Callable[[int], Regressor]
102
103
104def estimators(opaque: OpaqueConfig) -> tuple[tuple[str, Builder], ...]:
105    """The published regressors, as (label, factory).
106
107    Factories rather than instances so each fold gets a fresh unfitted estimator -- refitting
108    a fitted scikit-learn estimator is fine, but sharing one across folds is the kind of thing
109    that silently leaks if an estimator ever caches.
110
111    **The size is in the label**, so a table produced at a reduced size says so rather than
112    claiming a size it never grew.
113    """
114    return (
115        # A wide log-spaced grid rather than a chosen penalty: the point of the control is
116        # that a *well-regularised* linear model on raw columns still transfers badly, so
117        # picking the penalty by hand would leave the result open to that objection.
118        ("RidgeCV (linear)", lambda _threads=-1: RidgeCV(alphas=np.logspace(-3, 3, 25))),
119        (
120            f"RandomForest ({opaque.forest_trees} trees)",
121            lambda threads=-1: RandomForestRegressor(
122                n_estimators=opaque.forest_trees,
123                max_features=opaque.forest_max_features,
124                random_state=SEED,
125                n_jobs=threads,
126            ),
127        ),
128        (
129            f"GradientBoosting ({opaque.boosting_stages} stages)",
130            lambda _threads=-1: GradientBoostingRegressor(
131                n_estimators=opaque.boosting_stages,
132                learning_rate=opaque.boosting_learning_rate,
133                max_depth=opaque.boosting_max_depth,
134                random_state=SEED,
135            ),
136        ),
137    )
138
139
140def _design(frame: pl.DataFrame) -> np.ndarray:
141    """The eighteen raw meta-feature columns, in `ALL_FEATURES` order.
142
143    Raw and untransformed on purpose. The grammar's transforms are part of what this study
144    contributes, so handing them to the opaque side would be comparing the equation against
145    itself; handing the opaque side nothing but the columns is the comparison a reader means
146    by "would a forest have done better".
147    """
148    columns = columns_as_arrays(frame, ALL_FEATURES)
149    return np.column_stack([columns[name] for name in ALL_FEATURES])
150
151
152def _cross_validate(design: np.ndarray, truth: np.ndarray, labels: np.ndarray, build: Builder) -> np.ndarray:
153    """Out-of-fold predictions under one leave-one-group-out split.
154
155    Clipped to the training fold's own observed range, exactly as
156    `validate._clip_to_training` does for every reported number. Without it the ridge control
157    is scored on an unbounded extrapolation the equation is never scored on, which would make
158    the comparison unfair in the equation's favour.
159    """
160    predictions = np.zeros_like(truth)
161    for _, train, test in leave_one_group_out(labels):
162        estimator = build(-1)
163        estimator.fit(design[train], truth[train])
164        held = np.asarray(estimator.predict(design[test]), dtype=np.float64)
165        predictions[test] = np.clip(held, float(truth[train].min()), float(truth[train].max()))
166    return predictions
167
168
169@dataclass(frozen=True)
170class OpaqueRun:
171    """One pass over the opaque regressors: the score table, and the predictions behind it.
172
173    Both together because the folds are the expensive part -- a random forest refits 45 times
174    across the two protocols -- and the ranking and threshold comparisons need the same
175    predictions the regression table was scored from. Computing them twice would double the
176    cost and, worse, leave open the possibility of the two disagreeing.
177    """
178
179    table: pl.DataFrame
180    #: ``predictions[label][protocol]``, for every estimator `estimators` names and all four
181    #: protocols -- ``in_sample``, ``loo_dataset``, ``loo_model`` and ``loo_cell``. Clipped
182    #: exactly as the reported numbers are.
183    predictions: dict[str, dict[str, np.ndarray]]
184
185
186def evaluate(frame: pl.DataFrame, *, models: tuple[tuple[str, Builder], ...] | None = None) -> OpaqueRun:
187    """Fit and score every opaque regressor once, keeping the predictions.
188
189    **Read the IS and DHO columns of the forest row together.** That pair is
190    the study's argument in two numbers: a flexible model fits this meta-data almost perfectly
191    and, with neither the dataset nor the model in training, reaches nothing. With twenty
192    dataset groups and features constant within a group it identifies the dataset and looks the
193    answer up, and identification is worth nothing on a cell nobody has run.
194
195    It is also the likely provenance of the R2 ~ 0.9 figures reported for opaque meta-models
196    elsewhere: a forest evaluated under IS or a random split reproduces them exactly.
197    """
198    design = _design(frame)
199    truth = target(frame)
200    datasets, model_labels = groups(frame, DATASET_COLUMN), groups(frame, MODEL_COLUMN)
201    splits = {"loo_dataset": datasets, "loo_model": model_labels}
202
203    rows: list[dict[str, object]] = []
204    predictions: dict[str, dict[str, np.ndarray]] = {}
205    for label, build in models or estimators(load_config().opaque):
206        fitted = build(-1)
207        fitted.fit(design, truth)
208        held: dict[str, np.ndarray] = {
209            "in_sample": np.clip(np.asarray(fitted.predict(design), dtype=np.float64), MCC_LOWER, MCC_UPPER)
210        }
211        for name, labels in splits.items():
212            held[name] = _cross_validate(design, truth, labels, build)
213        held["loo_cell"] = _doubly_held_out(design, truth, datasets, model_labels, build)
214        predictions[label] = held
215        rows.append(
216            {
217                "model": label,
218                "features": len(ALL_FEATURES),
219                **{
220                    key: value
221                    for name, values in held.items()
222                    for key, value in ((f"r2_{name}", r2_score(truth, values)), (f"mae_{name}", mae(truth, values)))
223                },
224            }
225        )
226    return OpaqueRun(table=pl.DataFrame(rows), predictions=predictions)
227
228
229def _doubly_held_out(
230    design: np.ndarray,
231    truth: np.ndarray,
232    first: np.ndarray,
233    second: np.ndarray,
234    build: Builder,
235    jobs: int = -1,
236) -> np.ndarray:
237    """Predictions with **both** the dataset and the model of each cell removed from training.
238
239    The protocol the study is actually for -- *what will this model reach on this dataset*,
240    when neither has been run -- and the only one under which the equation and an opaque
241    regressor are denied the same things. Under LODO a forest still has the
242    held-out learner on nineteen other problems; under LOMO it still has the
243    held-out dataset. Here it has neither, and neither does the equation.
244
245    That makes this the like-for-like comparison in the study, and it is also the protocol on
246    which the trivial per-model baselines cannot be computed at all: a model held out of every
247    fold has no rows to average. A feature-based predictor still predicts.
248
249    One refit per observed cell rather than per group -- 484 of them against 20 or 25 -- which
250    for a forest is the difference between seconds and minutes, so the cells are run in
251    parallel. Each estimator is given a single thread for the same reason: at 434 training rows
252    the thread pool costs more than the trees, and 484 pools of 16 would oversubscribe the
253    machine by an order of magnitude.
254    """
255    cells = [
256        (row, column)
257        for row in np.unique(first)
258        for column in np.unique(second)
259        if ((first == row) & (second == column)).any()
260    ]
261
262    def one(row: object, column: object) -> tuple[np.ndarray, np.ndarray]:
263        test = (first == row) & (second == column)
264        train = (first != row) & (second != column)
265        if not train.any():
266            return test, np.zeros(int(test.sum()))
267        estimator = build(1)
268        estimator.fit(design[train], truth[train])
269        held = np.asarray(estimator.predict(design[test]), dtype=np.float64)
270        return test, np.clip(held, float(truth[train].min()), float(truth[train].max()))
271
272    # A threading backend rather than processes: the estimators release the GIL inside their
273    # own fits, and 484 process spawns would each have to pickle the design across.
274    completed: list[tuple[np.ndarray, np.ndarray]] = Parallel(n_jobs=jobs, backend="threading")(
275        delayed(one)(row, column) for row, column in cells
276    )  # pyright: ignore[reportAssignmentType]
277    predictions = np.zeros_like(truth)
278    for test, values in completed:
279        predictions[test] = values
280    return predictions