ml_meta_perf.analysis

  1"""Correlation screening: which terms carry signal, and which merely repeat each other.
  2
  3Two numbers are reported per term. The plain correlation against MCC is the obvious
  4one. The within-group correlation is the one that actually decides whether a term is
  5worth anything: both the term and the target are centred inside each group before
  6correlating, so a term is credited only for variance that group identity does not
  7already explain. A dataset feature scored within-model, or a model feature scored
  8within-dataset, has to earn its correlation the hard way.
  9
 10Study chapter: [3. Term generation and selection][study-chapter] -- the rationale, in
 11prose, with the figures.
 12
 13[study-chapter]: https://github.com/mariolpantunes/ml-meta-perf/blob/main/assets/docs/03-term-selection.md
 14"""
 15
 16from __future__ import annotations
 17
 18import numpy as np
 19import polars as pl
 20
 21from ml_meta_perf.fit import Standardizer
 22from ml_meta_perf.stats import pearson, r2_score, spearman
 23from ml_meta_perf.terms import Library
 24from ml_meta_perf.validate import leave_one_group_out
 25
 26
 27def _center_within_groups(values: np.ndarray, groups: np.ndarray) -> np.ndarray:
 28    """Subtract each group's own mean, leaving only within-group variation."""
 29    centered = values.astype(np.float64).copy()
 30    for label in np.unique(groups):
 31        mask = groups == label
 32        centered[mask] -= centered[mask].mean()
 33    return centered
 34
 35
 36def screen(
 37    library: Library,
 38    target: np.ndarray,
 39    groups: np.ndarray | None = None,
 40) -> pl.DataFrame:
 41    """Rank every term in the library by its association with the target.
 42
 43    ``groups`` enables the within-group columns. Sorting is by absolute within-group
 44    Pearson when it is available, and by absolute Pearson otherwise.
 45    """
 46    matrix = library.matrix
 47    rows: list[dict[str, object]] = []
 48    within_target = _center_within_groups(target, groups) if groups is not None else None
 49
 50    for index, name in enumerate(library.names):
 51        column = matrix[:, index]
 52        record: dict[str, object] = {
 53            "term": name,
 54            "pearson": pearson(column, target),
 55            "spearman": spearman(column, target),
 56        }
 57        if within_target is not None:
 58            within_column = _center_within_groups(column, groups)  # pyright: ignore[reportArgumentType]
 59            record["within_pearson"] = pearson(within_column, within_target)
 60        rows.append(record)
 61
 62    frame = pl.DataFrame(rows)
 63    key = "within_pearson" if within_target is not None else "pearson"
 64    return frame.with_columns(pl.col(key).abs().alias("strength")).sort("strength", descending=True)
 65
 66
 67def redundancy_groups(library: Library, threshold: float = 0.99) -> list[list[str]]:
 68    """Cluster terms that are near-duplicates of one another.
 69
 70    Purely diagnostic: it explains why a library of a thousand terms behaves like a much
 71    smaller one, and why the collinearity guard during selection is not optional.
 72    """
 73    matrix = library.matrix
 74    normalized = matrix - matrix.mean(axis=0)
 75    scale = np.sqrt((normalized**2).sum(axis=0))
 76    scale = np.where(scale < 1e-12, 1.0, scale)
 77    normalized = normalized / scale
 78    correlation = np.abs(normalized.T @ normalized)
 79
 80    names = library.names
 81    unassigned = set(range(len(names)))
 82    clusters: list[list[str]] = []
 83    while unassigned:
 84        seed = min(unassigned)
 85        members = sorted(index for index in unassigned if correlation[seed, index] >= threshold)
 86        unassigned -= set(members)
 87        if len(members) > 1:
 88            clusters.append([names[index] for index in members])
 89    return clusters
 90
 91
 92def feature_reach(
 93    library: Library,
 94    target: np.ndarray,
 95    features: tuple[str, ...],
 96) -> pl.DataFrame:
 97    """What each raw feature is worth on its own, raw and after the grammar has worked on it.
 98
 99    Two questions the equation itself cannot answer, because by the time an equation exists
100    every feature is entangled with the others:
101
102    * **How much of MCC does this feature carry at all?** ``r2_raw`` is the squared
103      correlation of the untransformed column with the target -- what a straight line through
104      that one feature would explain. It is NaN when the grammar's outlier guard rejects the
105      raw column while retaining a compressed transform.
106    * **How much of that does the grammar unlock?** ``r2_best`` is the same quantity for the
107      best *single-feature* term the vocabulary can build over it -- ``log(f)``, ``1/f``,
108      ``f^2`` and the rest. ``gain`` is the difference, and it is the only place the value of
109      the transforms is visible in isolation.
110
111    Terms touching more than one feature are excluded on purpose: a product of two features
112    would be credited to both, and the point here is the marginal worth of each one.
113
114    This is a **heuristic**, not a bound. It ignores every interaction the equation is built
115    out of, so an equation can and does exceed the sum of these parts -- which is the useful
116    reading, since the excess is exactly what the multi-feature terms contribute.
117    """
118    rows: list[dict[str, object]] = []
119    columns = {name: library.matrix[:, index] for index, name in enumerate(library.names)}
120
121    for feature in features:
122        candidates = [
123            term.name for term in library.terms if term.features == (feature,) or set(term.features) == {feature}
124        ]
125        raw = columns.get(feature)
126        raw_r = pearson(raw, target) if raw is not None else float("nan")
127        best_name, best_r = "", 0.0
128        for name in candidates:
129            value = pearson(columns[name], target)
130            if abs(value) > abs(best_r):
131                best_name, best_r = name, value
132        rows.append(
133            {
134                "feature": feature,
135                "n_terms": len(candidates),
136                "raw_pearson": raw_r,
137                "best_term": best_name,
138                "best_pearson": best_r,
139                "r2_raw": raw_r**2 if raw_r == raw_r else float("nan"),
140                "r2_best": best_r**2,
141                "gain": best_r**2 - (raw_r**2 if raw_r == raw_r else 0.0),
142            }
143        )
144    return pl.DataFrame(rows).sort("r2_best", descending=True)
145
146
147def grammar_ceiling(
148    library: Library,
149    target: np.ndarray,
150    features: tuple[str, ...],
151    *,
152    penalty: float = 0.0,
153) -> dict[str, float]:
154    """A heuristic estimate of how far the grammar reaches, before any search runs.
155
156    Three levels, each a least-squares fit, each computable from the library alone:
157
158    * ``r2_raw_additive`` -- every untransformed feature term admitted by the grammar's
159      outlier guard. This is the raw-term rung of the vocabulary, not an external model.
160    * ``r2_best_per_feature`` -- the best admissible single-feature term per feature, one
161      each. What the vocabulary adds over its raw-term rung, with no interactions yet.
162    * ``r2_all_single_feature`` -- *every* single-feature term at once. The most an equation
163      with no cross-feature terms could explain, and therefore the level above which an
164      equation's accuracy has to be coming from interactions.
165
166    The gap between the last of these and the fitted equation is the part of the study that
167    the additive mean-based reference cannot see: `validate.additive_mean_reference` combines a per-dataset value plus
168    a per-model value, whereas this bounds a sum of per-*feature* functions. An equation above
169    ``r2_all_single_feature`` is one whose mixed terms are doing real work.
170
171    Fitted in-sample (IS) and unregularised by default: it is a description of the vocabulary's
172    reach on this data, not a predictor, and shrinkage would understate it.
173    """
174    index = {name: position for position, name in enumerate(library.names)}
175
176    def fit_r2(names: list[str]) -> float:
177        present = [index[name] for name in names if name in index]
178        if not present:
179            return float("nan")
180        design = library.matrix[:, present]
181        design = design - design.mean(axis=0)
182        scale = design.std(axis=0)
183        design = design / np.where(scale < 1e-12, 1.0, scale)
184        gram = design.T @ design + penalty * np.eye(design.shape[1])
185        offset = float(target.mean())
186        weights = np.linalg.solve(gram, design.T @ (target - offset))
187        return r2_score(target, design @ weights + offset)
188
189    single = [term for term in library.terms if len(set(term.features)) == 1]
190    reach = feature_reach(library, target, features)
191    return {
192        "r2_raw_additive": fit_r2([name for name in features if name in index]),
193        "r2_best_per_feature": fit_r2([str(name) for name in reach["best_term"].to_list() if name]),
194        "r2_all_single_feature": fit_r2([term.name for term in single]),
195        "n_single_feature_terms": float(len(single)),
196    }
197
198
199def _min_norm_lstsq(design: np.ndarray, target: np.ndarray) -> np.ndarray:
200    """Minimum-norm least squares that does not depend on one LAPACK build converging.
201
202    `np.linalg.lstsq` raises "SVD did not converge" on some OpenBLAS builds for the
203    saturated design, which is numerically rank-deficient (about 100 independent columns of
204    229 in the leave-X-IIoTID-out fold). The SVD of the transpose is the same problem and
205    converges there, so it is the fallback; the truncation is `lstsq`'s own default.
206    """
207    try:
208        return np.linalg.lstsq(design, target, rcond=None)[0]
209    except np.linalg.LinAlgError:
210        right, singular, left = np.linalg.svd(design.T, full_matrices=False)
211        keep = singular > singular[0] * np.finfo(float).eps * max(design.shape)
212        return (right[:, keep] / singular[keep]) @ (left[keep] @ target)
213
214
215def saturated_fit(
216    library: Library,
217    target: np.ndarray,
218    groups: np.ndarray,
219) -> dict[str, float]:
220    """What happens if every candidate term is handed to least squares at once.
221
222    The control for the whole selection stage. If an unpenalised fit over the entire library
223    transferred well, the beam search and the length rule would be machinery in search of a
224    problem, and the honest thing would be to delete them. It does not: the design is far
225    wider than 20 held-out groups can support, so the fit is excellent and the transfer is
226    catastrophic.
227
228    Reported at both bounds because the gap between them is itself the finding. Unclipped, a
229    held-out dataset outside the convex hull of the others is extrapolated without limit and
230    the pooled R2 goes to three or four figures negative; clipped to the training fold's own
231    observed range -- exactly what `validate._clip_to_training` does for every reported
232    number -- it lands somewhere merely bad. Quoting only the clipped figure would hide how
233    far the unconstrained fit actually goes.
234
235    No ridge penalty and no screening: this is the null procedure the study is measured
236    against, so giving it either would be measuring something else.
237    """
238    matrix = library.matrix
239    standardizer = Standardizer.fit(matrix)
240    offset = float(target.mean())
241    weights = _min_norm_lstsq(standardizer.apply(matrix), target - offset)
242    in_sample = r2_score(target, standardizer.apply(matrix) @ weights + offset)
243
244    held = np.zeros_like(target)
245    bounded = np.zeros_like(target)
246    for _, train, test in leave_one_group_out(groups):
247        fold = Standardizer.fit(matrix[train])
248        design = fold.apply(matrix[train])
249        centre = float(target[train].mean())
250        fitted = _min_norm_lstsq(design, target[train] - centre)
251        held[test] = fold.apply(matrix[test]) @ fitted + centre
252        bounded[test] = np.clip(held[test], float(target[train].min()), float(target[train].max()))
253    return {
254        "terms": float(len(library)),
255        "r2_in_sample": in_sample,
256        "r2_loo_dataset_clipped": r2_score(target, bounded),
257        "r2_loo_dataset_unclipped": r2_score(target, held),
258    }