ml_meta_perf.identity

  1"""Model identity: the part of a classifier the meta-features do not describe.
  2
  3Chapter 4 measures a gap. Dataset identity explains 0.354 of MCC variance and the twelve
  4dataset meta-features recover 99% of that; model identity explains 0.282 and the six model
  5meta-features need not encode every distinction between learner histories --
  6the difference between *what a model costs*, which the features state, and *what a model
  7is good at*, which they do not.
  8
  9The remaining headroom for a per-model table is small, but still useful as a ceiling:
 10it bounds what any model-side proposal could recover beyond the generated descriptors.
 11The module is kept for that measurement, not because the correction it computes is worth
 12publishing -- see the standing objections in chapter 4.
 13
 14This module writes that difference down. Under leave-one-dataset-out (LODO), every one of the 25
 15models appears in every training fold, so the residual of the fitted equation can be
 16averaged per model and the average used on the held-out dataset. Nothing about the
 17held-out dataset enters it, and the result is two numbers per model:
 18
 19    MCC = [the equation]  +  b_model  +  c_model * log(gravity)
 20
 21``b`` is a level -- this model runs above or below what its meta-features predict -- and
 22``c`` is a slope, how that offset changes with how separable the data is. The pair is a
 23*table*, not an equation, and that is the honest cost of the accuracy it buys.
 24
 25**Two properties of this construction decide how it may be used.**
 26
 27It is not a predictor for an unseen *model*. A model with no training rows has no
 28effect, `ModelEffects.apply` returns zero for it, and the correction degenerates to the
 29equation it corrects. The leave-one-model-out (LOMO) column measures exactly that and shows no
 30gain, by construction rather than by accident.
 31
 32It is fitted on the equation's residual, not jointly with it. Alternating the two would
 33turn the measurement into a backfitted model with a different training procedure
 34(Hastie & Tibshirani, 1990). The reported experiment does not fit or evaluate that
 35variant. The equation remains a statement about MCC, and the table explains only what
 36the finished equation leaves behind.
 37
 38The construction is a **factorial regression** in the sense used for genotype-by-
 39environment trials: a two-way table modelled with covariates on one side and free
 40coefficients on the other, one of which multiplies an environmental covariate (Denis,
 411988; van Eeuwijk, Denis & Kang, 1996). Where the AMMI oracle of chapter 4 estimates
 42both sides freely and therefore predicts nothing, this estimates the dataset side from
 43meta-features and keeps the model side free -- which is what makes it usable on a
 44dataset nobody has run. In machine-learning terms it is the collaborative half of a
 45hybrid recommender: the equation is content-based over meta-features, the table is a
 46per-item effect learned from the observed grid, and algorithm selection has been posed
 47as collaborative filtering before (Misir & Sebag, 2017; Fusi et al., 2018; Yang et al.,
 482019).
 49
 50Study chapter: [4. The equation][study-chapter]
 51-- the rationale, in prose, with the figures.
 52
 53[study-chapter]: https://github.com/mariolpantunes/ml-meta-perf/blob/main/assets/docs/04-equation.md#limitations-of-the-equation
 54"""
 55
 56from __future__ import annotations
 57
 58from dataclasses import dataclass
 59
 60import numpy as np
 61import polars as pl
 62
 63from ml_meta_perf.model import MCC_LOWER, MCC_UPPER, Equation
 64from ml_meta_perf.terms import Atom, composition_atom
 65from ml_meta_perf.validate import CrossValidation, leave_one_group_out
 66
 67#: Denominator added to a model's row count when averaging its residuals. A model with
 68#: 19 rows is shrunk by 19/24 toward zero, which is the empirical-Bayes estimator for a
 69#: group mean under a common prior (Efron & Morris, 1975) with the prior-to-noise ratio
 70#: fixed rather than estimated. Fixed because estimating it from 20 groups is noisier
 71#: than choosing it: LODO R² moves by 0.015 across 0 to 20, which is
 72#: less than the fold noise it would be tuned against.
 73INTERCEPT_SHRINKAGE = 5.0
 74
 75#: The same idea for the slopes, in units of the standardised carrying feature. R2 moves
 76#: by 0.005 across 0.5 to 10.
 77SLOPE_SHRINKAGE = 1.0
 78
 79
 80@dataclass(frozen=True)
 81class ModelEffects:
 82    """Per-model level and slope, in the units the equation is written in.
 83
 84    ``intercepts`` and ``slopes`` are keyed by model name and read directly: a prediction
 85    for model *m* on a dataset is the equation's value plus ``intercepts[m]`` plus
 86    ``slopes[m]`` times the carrying feature. Models absent from the table contribute
 87    nothing, so an equation corrected this way never becomes *worse* than uncorrected on
 88    a model it has not seen -- it simply stops being corrected.
 89    """
 90
 91    #: The dataset feature the slope multiplies, and the transform it enters under.
 92    #: ``None`` when the effects are levels only.
 93    atom: Atom | None
 94    intercepts: dict[str, float]
 95    slopes: dict[str, float]
 96
 97    @property
 98    def n_models(self) -> int:
 99        return len(self.intercepts)
100
101    def apply(self, columns: dict[str, np.ndarray], models: np.ndarray) -> np.ndarray:
102        """The correction to add to the equation's output, row by row."""
103        levels = np.array([self.intercepts.get(str(name), 0.0) for name in models])
104        if self.atom is None:
105            return levels
106        carrier = self.atom.evaluate(columns)
107        gradient = np.array([self.slopes.get(str(name), 0.0) for name in models])
108        return levels + gradient * carrier
109
110    def table(self) -> pl.DataFrame:
111        """One row per model, largest level first.
112
113        ``effect_at_median`` is the whole correction evaluated at the median of the
114        carrying feature, which is the number to compare against MCC: a level and a
115        slope in different units cannot be ranked against each other, and reading the
116        level alone overstates models whose slope cancels it.
117        """
118        names = sorted(self.intercepts, key=lambda name: -self.intercepts[name])
119        return pl.DataFrame(
120            {
121                "model": names,
122                "level": [self.intercepts[name] for name in names],
123                "slope": [self.slopes.get(name, 0.0) for name in names],
124            }
125        )
126
127    def __str__(self) -> str:
128        carrier = "" if self.atom is None else f" + c_model * {self.atom.name}"
129        lines = [f"correction = b_model{carrier}"]
130        for row in self.table().iter_rows(named=True):
131            lines.append(f"      {row['model']:<22s} b={row['level']:+.4f}  c={row['slope']:+.4f}")
132        return "\n".join(lines)
133
134
135def carrying_atoms(features: tuple[str, ...], columns: dict[str, np.ndarray]) -> list[Atom]:
136    """The candidate carriers for the slope, one per dataset feature.
137
138    Each feature enters under the library's own compression rule -- log when strictly
139    positive, raw otherwise -- so the slope is read in the same units as any composite
140    term in the equation, and no second convention is introduced here.
141    """
142    return [composition_atom(feature, columns) for feature in features]
143
144
145def fit_effects(
146    residual: np.ndarray,
147    columns: dict[str, np.ndarray],
148    models: np.ndarray,
149    features: tuple[str, ...] = (),
150    *,
151    intercept_shrinkage: float = INTERCEPT_SHRINKAGE,
152    slope_shrinkage: float = SLOPE_SHRINKAGE,
153) -> ModelEffects:
154    """Fit per-model levels, and one per-model slope, to an equation's residual.
155
156    Two stages, in this order. The level is the shrunk mean residual of the model's own
157    rows. The slope is then fitted to what the levels leave behind, and the feature that
158    carries it is chosen by training residual sum of squares over every candidate in
159    ``features`` -- which is a search, so it happens inside the fold like any other.
160
161    Pass no ``features`` for levels only.
162    """
163    names = [str(name) for name in np.unique(models)]
164    labels = np.array([str(name) for name in models])
165    intercepts = {
166        name: float(residual[labels == name].sum() / (int(np.sum(labels == name)) + intercept_shrinkage))
167        for name in names
168    }
169    inner = residual - np.array([intercepts[name] for name in labels])
170
171    candidates = carrying_atoms(features, columns) if features else []
172    best: tuple[float, Atom | None, dict[str, float]] = (float(inner @ inner), None, {})
173    for atom in candidates:
174        with np.errstate(all="ignore"):
175            carrier = atom.evaluate(columns)
176        if not np.all(np.isfinite(carrier)):
177            continue
178        # Centred and scaled so the shrinkage means the same thing whichever feature is
179        # carrying, then folded back into raw units below.
180        centre, spread = float(carrier.mean()), float(carrier.std())
181        if spread < 1e-12:
182            continue
183        scaled = (carrier - centre) / spread
184        gradient = {
185            name: float(
186                scaled[labels == name]
187                @ inner[labels == name]
188                / (float(scaled[labels == name] @ scaled[labels == name]) + slope_shrinkage)
189            )
190            for name in names
191        }
192        fitted = np.array([gradient[name] for name in labels]) * scaled
193        rss = float((inner - fitted) @ (inner - fitted))
194        if rss < best[0]:
195            raw = {name: value / spread for name, value in gradient.items()}
196            best = (rss, atom, raw)
197
198    _, atom, slopes = best
199    if atom is None:
200        return ModelEffects(atom=None, intercepts=intercepts, slopes={})
201    # The slope's centring belongs in the level, so that both numbers read in raw units.
202    centre = float(atom.evaluate(columns).mean())
203    shifted = {name: intercepts[name] - slopes[name] * centre for name in names}
204    return ModelEffects(atom=atom, intercepts=shifted, slopes=slopes)
205
206
207def correct_out_of_fold(
208    fold: CrossValidation,
209    columns: dict[str, np.ndarray],
210    target: np.ndarray,
211    groups: np.ndarray,
212    models: np.ndarray,
213    features: tuple[str, ...] = (),
214    *,
215    intercept_shrinkage: float = INTERCEPT_SHRINKAGE,
216    slope_shrinkage: float = SLOPE_SHRINKAGE,
217) -> np.ndarray:
218    """Re-score a finished cross-validation with per-model effects fitted in each fold.
219
220    The effects are cheap and the beam search is not, so this consumes the per-fold
221    equations a completed `ml_meta_perf.validate.cross_validate_fixed_form` already recorded rather
222    than searching again. Each fold's effects are fitted on that fold's training rows
223    only, from that fold's own equation, and applied to the rows it held out -- so the
224    protocol is the one the path was run under, unchanged.
225    """
226    predictions = np.zeros_like(target)
227    for label, train, test in leave_one_group_out(groups):
228        equation = fold.equations.get(label)
229        if equation is None:
230            predictions[test] = fold.predictions[test]
231            continue
232        trained = {name: values[train] for name, values in columns.items()}
233        held = {name: values[test] for name, values in columns.items()}
234        effects = fit_effects(
235            target[train] - equation.evaluate(trained),
236            trained,
237            models[train],
238            features,
239            intercept_shrinkage=intercept_shrinkage,
240            slope_shrinkage=slope_shrinkage,
241        )
242        corrected = equation.evaluate(held) + effects.apply(held, models[test])
243        predictions[test] = np.clip(corrected, MCC_LOWER, MCC_UPPER)
244    return predictions
245
246
247def predict(
248    equation: Equation,
249    effects: ModelEffects,
250    columns: dict[str, np.ndarray],
251    models: np.ndarray,
252) -> np.ndarray:
253    """Equation plus table, clipped into the range MCC can take."""
254    return np.clip(equation.evaluate(columns) + effects.apply(columns, models), MCC_LOWER, MCC_UPPER)
255
256
257def carrier_stability(
258    columns: dict[str, np.ndarray],
259    target: np.ndarray,
260    groups: np.ndarray,
261    models: np.ndarray,
262    fold: CrossValidation,
263    features: tuple[str, ...],
264) -> pl.DataFrame:
265    """Which feature carried the slope, counted across folds.
266
267    The carrier is chosen inside each fold, so it is a fitted quantity like any selected
268    term and reporting it without its selection frequency would present a choice made 20
269    times out of 20 identically to one made 6 times.
270    """
271    counts: dict[str, int] = {}
272    for label, train, _ in leave_one_group_out(groups):
273        equation = fold.equations.get(label)
274        if equation is None:
275            continue
276        trained = {name: values[train] for name, values in columns.items()}
277        effects = fit_effects(target[train] - equation.evaluate(trained), trained, models[train], features)
278        name = "none" if effects.atom is None else effects.atom.name
279        counts[name] = counts.get(name, 0) + 1
280    total = max(sum(counts.values()), 1)
281    return (
282        pl.DataFrame({"carrier": list(counts), "folds": list(counts.values())})
283        .with_columns((pl.col("folds") / total).alias("frequency"))
284        .sort("folds", descending=True)
285    )