ml_meta_perf.stats

  1"""The handful of statistics the project needs, written out rather than imported.
  2
  3scipy would supply all of this, but it is a large dependency for four short functions,
  4and the rank machinery here has to be exercised by the test suite regardless.
  5
  6Study chapter: [5. Evaluation][study-chapter] -- the rationale, in
  7prose, with the figures.
  8
  9[study-chapter]: https://github.com/mariolpantunes/ml-meta-perf/blob/main/assets/docs/05-evaluation.md
 10"""
 11
 12from __future__ import annotations
 13
 14import numpy as np
 15
 16
 17def rankdata(values: np.ndarray) -> np.ndarray:
 18    """Ranks from 1, with ties receiving their average rank.
 19
 20    Average ranks matter here: MCC can pile up at exactly 1.0,
 21    and ordinal ranks would invent an ordering among tied perfect scores.
 22    """
 23    order = np.argsort(values, kind="stable")
 24    ranks = np.empty(values.shape[0], dtype=np.float64)
 25    ranks[order] = np.arange(1, values.shape[0] + 1, dtype=np.float64)
 26    unique, inverse, counts = np.unique(values, return_inverse=True, return_counts=True)
 27    if unique.shape[0] != values.shape[0]:
 28        sums = np.zeros(unique.shape[0], dtype=np.float64)
 29        np.add.at(sums, inverse, ranks)
 30        ranks = (sums / counts)[inverse]
 31    return ranks
 32
 33
 34def rank_columns(matrix: np.ndarray) -> np.ndarray:
 35    """`rankdata` applied down every column at once, with the same tie handling.
 36
 37    `ml_meta_perf.search.guided_screen` needs the ranks of every candidate term in the library, in
 38    every fold. Doing that a column at a time cost 32k calls and a tenth of the study's
 39    runtime; the work is identical but the Python loop is not. ``test_stats`` asserts the
 40    two agree column by column, so `rankdata` stays the definition and this stays a
 41    restatement of it.
 42
 43    Ties are averaged by locating each run of equal values in the sorted order: a run
 44    spanning sorted positions ``i..j`` takes rank ``(i + j) / 2 + 1``, which is the mean of
 45    the consecutive ranks it would otherwise receive.
 46    """
 47    rows = matrix.shape[0]
 48    order = np.argsort(matrix, axis=0, kind="stable")
 49    ordered = np.take_along_axis(matrix, order, axis=0)
 50
 51    starts = np.empty(matrix.shape, dtype=bool)
 52    starts[0] = True
 53    np.not_equal(ordered[1:], ordered[:-1], out=starts[1:])
 54    ends = np.empty(matrix.shape, dtype=bool)
 55    ends[-1] = True
 56    ends[:-1] = starts[1:]
 57
 58    position = np.arange(rows, dtype=np.float64)[:, None]
 59    first = np.maximum.accumulate(np.where(starts, position, -1.0), axis=0)
 60    last = np.minimum.accumulate(np.where(ends, position, float(rows))[::-1], axis=0)[::-1]
 61
 62    ranks = np.empty(matrix.shape, dtype=np.float64)
 63    np.put_along_axis(ranks, order, (first + last) / 2.0 + 1.0, axis=0)
 64    return ranks
 65
 66
 67def pearson(first: np.ndarray, second: np.ndarray) -> float:
 68    """Linear correlation, returning 0.0 when either side is constant."""
 69    a = first - first.mean()
 70    b = second - second.mean()
 71    scale = float(np.sqrt(float(a @ a) * float(b @ b)))
 72    return float(a @ b / scale) if scale > 1e-15 else 0.0
 73
 74
 75def pearson_columns(matrix: np.ndarray, target: np.ndarray) -> np.ndarray:
 76    """`pearson` between every column of ``matrix`` and ``target``, in one pass.
 77
 78    `ml_meta_perf.search.guided_screen` needs this correlation for every candidate term, twice
 79    over -- once raw and once on ranks -- in every fold, which is tens of thousands of calls
 80    to `pearson` for arithmetic that is a single matrix-vector product. Constant columns
 81    return 0.0 here exactly as they do there. ``test_stats`` asserts the two agree column by
 82    column, so `pearson` stays the definition and this stays a restatement of it.
 83    """
 84    centred = matrix - matrix.mean(axis=0)
 85    other = target - target.mean()
 86    norms = np.sqrt(np.einsum("ij,ij->j", centred, centred) * float(other @ other))
 87    with np.errstate(divide="ignore", invalid="ignore"):
 88        return np.where(norms > 1e-15, (centred.T @ other) / np.where(norms > 1e-15, norms, 1.0), 0.0)
 89
 90
 91def spearman(first: np.ndarray, second: np.ndarray) -> float:
 92    """Rank correlation: monotone association, insensitive to the MCC ceiling."""
 93    if first.shape[0] < 2:
 94        return 0.0
 95    return pearson(rankdata(first), rankdata(second))
 96
 97
 98def r2_score(truth: np.ndarray, prediction: np.ndarray) -> float:
 99    """Coefficient of determination against the mean of ``truth``.
100
101    Note the denominator is always the variance of the *evaluated* rows, so an R2 from
102    the 20 aggregated dataset means is not comparable with one from the 484 raw rows.
103    That is exactly why E1 is reported on both scales.
104    """
105    residual = float(((truth - prediction) ** 2).sum())
106    total = float(((truth - truth.mean()) ** 2).sum())
107    return 1.0 - residual / total if total > 1e-15 else 0.0
108
109
110def mae(truth: np.ndarray, prediction: np.ndarray) -> float:
111    """Mean absolute error, in MCC units."""
112    return float(np.abs(truth - prediction).mean())
113
114
115def rmse(truth: np.ndarray, prediction: np.ndarray) -> float:
116    """Root mean squared error, in MCC units."""
117    return float(np.sqrt(float(((truth - prediction) ** 2).mean())))
118
119
120def smape(truth: np.ndarray, prediction: np.ndarray) -> float:
121    """Symmetric mean absolute percentage error, as a percentage in [0, 200].
122
123    **Read this one with care on MCC.** SMAPE divides by ``|truth| + |prediction|``, and
124    rows with MCC exactly 0 each contribute the full
125    200% unless the prediction is also exactly 0, so the metric is dominated by the rows
126    the equation is already known to handle worst rather than by its typical error. It is
127    reported because it was asked for and because it is scale-free, but MAE is the
128    honest headline for a target that legitimately passes through zero.
129
130    Rows where both truth and prediction are zero contribute 0, not a division by zero.
131    """
132    denominator = np.abs(truth) + np.abs(prediction)
133    safe = np.where(denominator > 0.0, denominator, 1.0)
134    ratio = np.where(denominator > 0.0, 2.0 * np.abs(truth - prediction) / safe, 0.0)
135    return float(100.0 * ratio.mean())