diff --git a/evaluation/config.py b/evaluation/config.py index 1040fd2..f48ccec 100644 --- a/evaluation/config.py +++ b/evaluation/config.py @@ -1,77 +1,74 @@ -"""Configuration file for the OpenADMET CYP blind challenge.""" +"""Configuration file for the OpenADMET CYP blind challenge. + +Ported from the challenge backend, so it should match exactly what you see on the leaderboard! +""" + +from functools import partial -import numpy as np from scipy.stats import kendalltau, spearmanr from sklearn.metrics import ( accuracy_score, - balanced_accuracy_score, f1_score, - mean_absolute_error, matthews_corrcoef, + mean_absolute_error, + precision_score, r2_score, - roc_auc_score, + recall_score, ) +from .custom_scoring_functions import rae_soft_threshold_absolute_error -def rae(y_true, y_pred): - """Relative absolute error (RAE) metric for regression tasks.""" - return np.sum(np.abs(y_true - y_pred)) / np.sum(np.abs(y_true - np.mean(y_true))) - - -def _binarize(y_pred_proba, threshold: float = 0.5): - """Threshold predicted TDI probabilities into hard class labels.""" - return (np.asarray(y_pred_proba) >= threshold).astype(int) - - -def accuracy_from_proba(y_true, y_pred_proba): - """Accuracy of the TDI classifier after thresholding probabilities at 0.5.""" - return accuracy_score(y_true, _binarize(y_pred_proba)) - - -def balanced_accuracy_from_proba(y_true, y_pred_proba): - """Balanced accuracy of the TDI classifier after thresholding probabilities at 0.5.""" - return balanced_accuracy_score(y_true, _binarize(y_pred_proba)) - - -def f1_from_proba(y_true, y_pred_proba): - """F1 score of the TDI classifier after thresholding probabilities at 0.5.""" - return f1_score(y_true, _binarize(y_pred_proba)) - - -def mcc_from_proba(y_true, y_pred_proba): - """Matthews correlation coefficient after thresholding probabilities at 0.5.""" - return matthews_corrcoef(y_true, _binarize(y_pred_proba)) - - -def roc_auc_safe(y_true, y_pred_proba): - """ROC-AUC of the predicted TDI probabilities. - - Returns NaN if the bootstrap sample only contains one class, since - ROC-AUC is undefined in that case. - """ - if len(np.unique(y_true)) < 2: - return np.nan - return roc_auc_score(y_true, y_pred_proba) - +# Multi-endpoint macro-averaging . +# A pseudo-endpoint, scored alongside the real endpoints in every bootstrap sample, whose +# per-metric values are macro-averages across endpoints rather than raw per-endpoint scores. +MACRO_ENDPOINT_LABEL = "MA" # Activity dataset -ENDPOINTS = ["pEC50"] +IDENTIFIER_COLUMNS = ["SMILES", "Molecule_Name"] +REGRESSION_ENDPOINTS = [ + "CYP1A2_pIC50_direct_inhibition", + "CYP2C9_pIC50_direct_inhibition", + "CYP2D6_pIC50_direct_inhibition", + "CYP3A4_pIC50_direct_inhibition", +] +REGRESSION_CREDIBLE_INTERVALS_UPPER_SUFFIX = "_conf_high" +REGRESSION_CREDIBLE_INTERVALS_LOWER_SUFFIX = "_conf_low" +# Leave list empty if no classification tasks +CLASSIFICATION_ENDPOINTS = [ + "CYP2D6_is_TDI", + "CYP3A4_is_TDI", +] +ACTIVITY_ENDPOINTS = REGRESSION_ENDPOINTS + CLASSIFICATION_ENDPOINTS ENDPOINTS_TO_LOG_TRANSFORM: list[str] = [] +ACTIVITY_DATASET_SIZE = 750 ACTIVITY_METRICS = [ + ("ST-RAE", rae_soft_threshold_absolute_error), ("MAE", mean_absolute_error), - ("RAE", rae), ("R2", r2_score), - ("Spearman R", spearmanr), - ("Kendall's Tau", kendalltau), + ("Spearman_R", spearmanr), + ("Kendall_Tau", kendalltau), ] -BOOTSTRAP_SAMPLES = 1000 - -# TDI (time-dependent inhibition) classification dataset -TDI_ENDPOINT = "is_TDI" -TDI_METRICS = [ - ("Accuracy", accuracy_from_proba), - ("Balanced Accuracy", balanced_accuracy_from_proba), - ("F1", f1_from_proba), - ("MCC", mcc_from_proba), - ("ROC-AUC", roc_auc_safe), +# Rank correlations (Spearman_R, Kendall_Tau) are mathematically undefined +# whenever a bootstrap sample has zero variance in +# y_pred (e.g. a submission that predicts the same value for every compound) or y_true. +# 0.0 is the "no correlation" value on both metrics' [-1, 1] scale, matching how an +# unconditionally-constant predictor should be scored: no better than chance, not a +# hard failure. +METRIC_NAN_FALLBACK: dict[str, float] = { + "Spearman_R": 0.0, + "Kendall_Tau": 0.0, +} +# zero_division=0 matches sklearn's documented degenerate-case default, avoiding +# warnings/errors on bootstrap resamples with no positive predictions (TDI labels are +# imbalanced). matthews_corrcoef already returns 0.0 (not NaN) in its own degenerate +# case, so it needs no wrapping. +CLASSIFICATION_METRICS = [ + ("MCC", matthews_corrcoef), + ("Accuracy", accuracy_score), + ("Precision", partial(precision_score, zero_division=0)), + ("Recall", partial(recall_score, zero_division=0)), + ("F1", partial(f1_score, zero_division=0)), ] +SORT_REGRESSION_LEADERBOARD_BY = "ST-RAE" +SORT_CLASSIFICATION_LEADERBOARD_BY = "MCC" +BOOTSTRAP_SAMPLES = 1000 diff --git a/evaluation/custom_scoring_functions.py b/evaluation/custom_scoring_functions.py new file mode 100644 index 0000000..cfc1136 --- /dev/null +++ b/evaluation/custom_scoring_functions.py @@ -0,0 +1,300 @@ +"""Custom scoring functions for regression tasks.""" + +import numpy as np +import pandas as pd + +# Minimum confidence-interval width used in rae_inverse_confidence_weighting, to stop +# a near-zero interval from producing an arbitrarily large (or infinite) weight. +_MIN_CONFIDENCE_INTERVAL = 1e-6 + + +def rae(y_true: pd.Series | np.ndarray, y_pred: pd.Series | np.ndarray) -> float: + """Relative absolute error (RAE) metric for regression tasks. + + Args: + y_true (pd.Series | np.ndarray): True values. + y_pred (pd.Series | np.ndarray): Predicted values. + + Returns: + float: The relative absolute error (RAE) score. + + """ + return np.sum(np.abs(y_true - y_pred)) / np.sum(np.abs(y_true - np.mean(y_true))) + + +def _resolve_bounds( + y_true: pd.Series | np.ndarray, + y_true_upper: pd.Series | np.ndarray | None, + y_true_lower: pd.Series | np.ndarray | None, + confidence_interval: float | pd.Series | np.ndarray | None, +) -> tuple[pd.Series | np.ndarray, pd.Series | np.ndarray]: + """Resolve upper/lower bounds from either explicit bounds or a confidence interval. + + ``confidence_interval`` and explicit bounds (``y_true_upper``/``y_true_lower``) are + mutually exclusive. ``confidence_interval`` is treated as the full width of a band + centred on ``y_true`` (i.e. ``y_true +/- confidence_interval / 2``). + + A side of the band that is neither given explicitly nor derivable from + ``confidence_interval`` defaults to ``y_true`` itself — i.e. no tolerance on that + side, since there's no information to define one. If both sides default this way + (no bounds and no confidence_interval given at all), the band collapses to the + point estimate on both sides, which — for callers built on top of this, like + ``rae_soft_threshold_absolute_error`` — makes soft-thresholding a no-op and + recovers the plain (non-thresholded) behaviour exactly. + + Args: + y_true (pd.Series | np.ndarray): True values. + y_true_upper (pd.Series | np.ndarray | None): Optional upper bounds for true + values. + y_true_lower (pd.Series | np.ndarray | None): Optional lower bounds for true + values. + confidence_interval (float | pd.Series | np.ndarray | None): Optional + confidence interval (full width) for true values. + + Returns: + tuple[pd.Series | np.ndarray, pd.Series | np.ndarray]: The resolved + ``(y_true_lower, y_true_upper)`` bounds. + + Raises: + ValueError: If both explicit bounds and a confidence interval are provided. + + """ + if ( + y_true_upper is not None or y_true_lower is not None + ) and confidence_interval is not None: + raise ValueError( + "Cannot provide both upper/lower bounds and confidence interval for soft thresholding." + ) + + if confidence_interval is not None: + half_width = confidence_interval / 2 + y_true_upper = y_true + half_width + y_true_lower = y_true - half_width + else: + if y_true_upper is None: + y_true_upper = y_true + if y_true_lower is None: + y_true_lower = y_true + + return y_true_lower, y_true_upper + + +def rae_soft_threshold_absolute_error( + y_true: pd.Series | np.ndarray, + y_pred: pd.Series | np.ndarray, + y_true_upper: pd.Series | np.ndarray | None = None, + y_true_lower: pd.Series | np.ndarray | None = None, + confidence_interval: float | pd.Series | np.ndarray | None = None, +) -> float: + """RAE metric for regression tasks, with soft thresholding. + + If a confidence interval is provided, the absolute error is clipped to the distance + to the nearest bound defined by the confidence interval. The confidence interval can + be a single float or an array of the same shape as y_true/y_pred. If upper/lower + bounds are provided instead, the absolute error is clipped to the distance to the + nearest bound. Explicit bounds and a confidence interval are mutually exclusive. + + A prediction that falls inside the ``[y_true_lower, y_true_upper]`` tolerance band + contributes zero error — it is treated as indistinguishable from the true value + given measurement uncertainty. A prediction outside the band contributes only the + distance to the nearest edge, rather than the distance to the point estimate. + + Either ``y_true_upper`` or ``y_true_lower`` may be omitted on its own — that side + then defaults to ``y_true`` itself (no tolerance on that side; see + ``_resolve_bounds``). Omitting both, and no confidence_interval either, collapses + the band to a single point on both sides, which makes soft-thresholding a no-op: + the result is then identical to plain ``rae()``. + + The naive baseline in the denominator (a constant predictor at ``mean(y_true)``) is + put through the same soft-thresholding as the model's predictions, so both halves + of the ratio are computed under the same rule and RAE keeps its usual meaning: + ``1.0`` means the model is exactly as good as always predicting the mean, under + this tolerance-band error function. + + Args: + y_true (pd.Series | np.ndarray): True values. + y_pred (pd.Series | np.ndarray): Predicted values. + y_true_upper (pd.Series | np.ndarray | None): Optional upper bounds for true + values. Defaults to ``y_true`` (no tolerance above) if omitted. + y_true_lower (pd.Series | np.ndarray | None): Optional lower bounds for true + values. Defaults to ``y_true`` (no tolerance below) if omitted. + confidence_interval (float | pd.Series | np.ndarray | None): Optional + confidence interval (full width) for true values. + + Returns: + float: The relative absolute error (RAE) score, with soft-thresholded + absolute error in both the model error and the naive-baseline error. + + Raises: + ValueError: If both explicit bounds and a confidence interval are provided. + + """ + y_true_lower, y_true_upper = _resolve_bounds( + y_true, y_true_upper, y_true_lower, confidence_interval + ) + + above_upper = np.clip(y_pred - y_true_upper, a_min=0, a_max=None) + below_lower = np.clip(y_true_lower - y_pred, a_min=0, a_max=None) + soft_abs_error = above_upper + below_lower + + mean_true = np.mean(y_true) + baseline_above_upper = np.clip(mean_true - y_true_upper, a_min=0, a_max=None) + baseline_below_lower = np.clip(y_true_lower - mean_true, a_min=0, a_max=None) + soft_baseline_error = baseline_above_upper + baseline_below_lower + + return np.sum(soft_abs_error) / np.sum(soft_baseline_error) + + +def _weighted_absolute_error_below_threshold( + y_true: pd.Series | np.ndarray, + y_pred: pd.Series | np.ndarray | float, + threshold: float, + weighting: float, +) -> pd.Series | np.ndarray: + """Per-point absolute error, downweighted where both y_true and y_pred are below threshold. + + Args: + y_true (pd.Series | np.ndarray): True values. + y_pred (pd.Series | np.ndarray | float): Predicted values, or a single constant + prediction (e.g. ``mean(y_true)`` for a naive baseline) broadcast against + ``y_true``. + threshold (float): The threshold below which the absolute error is weighted. + weighting (float): The factor by which to weight the absolute error below the + threshold. + + Returns: + pd.Series | np.ndarray: Per-point weighted absolute error. + + """ + abs_error = np.abs(y_true - y_pred) + below_threshold = (y_true < threshold) & (y_pred < threshold) + weights = np.where(below_threshold, weighting, 1.0) + return weights * abs_error + + +def rae_weight_below_threshold( + y_true: pd.Series | np.ndarray, + y_pred: pd.Series | np.ndarray, + threshold: float = 4, + weighting: float = 0.25, +) -> float: + """RAE metric for regression tasks, weighted below a threshold. + + The absolute error is weighted by a factor if the ground truth and prediction are + both below a certain threshold. Requiring both — rather than just the true value — + to be below the threshold means a real, qualitative miss (e.g. true is below + threshold but predicted well above it, or vice versa) is still scored at full + weight; only errors confined to the low-confidence region are downweighted. + + The naive baseline in the denominator is treated as a constant predictor at + ``mean(y_true)`` and put through this *exact same* weighting rule (i.e. its + "prediction" for the below-threshold check is ``mean(y_true)`` itself, not + ``y_true`` alone) — so both halves of the ratio are computed identically and RAE + keeps its usual meaning: ``1.0`` means the model is exactly as good as always + predicting the mean, under this weighting. In practice, ``mean(y_true)`` for a + pIC50-like dataset is almost always above the threshold, so the baseline's + below-threshold condition is rarely satisfied and the denominator ends up close to + the unweighted RAE denominator — the weighting mostly changes the numerator. + + Args: + y_true (pd.Series | np.ndarray): True values. + y_pred (pd.Series | np.ndarray): Predicted values. + threshold (float): The threshold below which the absolute error is weighted. + Defaults to 4. + weighting (float): The factor by which to weight the absolute error below the + threshold. Defaults to 0.25. + + Returns: + float: The relative absolute error (RAE) score, with weighted absolute error + below the threshold in both the model error and the naive-baseline error. + + """ + mean_true = np.mean(y_true) + + weighted_error = _weighted_absolute_error_below_threshold( + y_true, y_pred, threshold, weighting + ) + weighted_baseline_error = _weighted_absolute_error_below_threshold( + y_true, mean_true, threshold, weighting + ) + + return np.sum(weighted_error) / np.sum(weighted_baseline_error) + + +def rae_inverse_confidence_weighting( + y_true: pd.Series | np.ndarray, + y_pred: pd.Series | np.ndarray, + y_true_upper: pd.Series | np.ndarray | None = None, + y_true_lower: pd.Series | np.ndarray | None = None, + confidence_interval: float | pd.Series | np.ndarray | None = None, +) -> float: + """RAE metric for regression tasks, with inverse confidence weighting. + + If a confidence interval is provided, the absolute error is weighted by the inverse + of the confidence interval. The confidence interval can be a single float or an + array of the same shape as y_true/y_pred. If upper/lower bounds are provided, the + absolute error is weighted by the inverse of the range of the upper and lower bounds. + If both or neither are provided, an error is raised. + + Unlike ``rae_soft_threshold_absolute_error``, this scales every point's + contribution continuously by how (un)certain its true value is, rather than + applying a hard tolerance band. + + The naive baseline in the denominator is the constant predictor that minimises the + *weighted squared* error — i.e. the weighted mean ``sum(weights * y_true) / + sum(weights)`` using these same per-point inverse-confidence weights — not the + plain, unweighted ``mean(y_true)``. This is the direct weighted analogue of how + plain RAE's baseline (the arithmetic mean) is used: it's the weighted + least-squares-optimal constant, *not* the weighted-L1-optimal constant (that would + be the weighted median) — plain RAE has this same quirk, using the mean rather + than the true minimiser of its own (L1) error function, so this preserves that + convention rather than fixing it. Since these weights depend only on ``y_true``'s + own uncertainty and not on any prediction, the weighted mean has a closed form, + unlike in ``rae_weight_below_threshold`` (where the weighting itself depends on the + candidate prediction, making a closed-form weighted mean ill-defined, so that + function keeps the plain mean instead). Both halves of the ratio are computed + under the same rule, so predicting the weighted mean exactly gives ``RAE = 1.0`` + by construction — though, as with plain RAE, some other constant (e.g. the + weighted median) can occasionally score lower than this baseline. + + A side of the bound left unspecified defaults to ``y_true`` itself (see + ``_resolve_bounds``) — e.g. omitting both bounds and confidence_interval entirely + collapses ``interval_width`` to ``0`` everywhere (clipped up to + ``_MIN_CONFIDENCE_INTERVAL``), giving every point the same weight, which reduces + this metric to plain ``rae()`` exactly (the constant weight cancels out of both + the numerator and denominator). + + Args: + y_true (pd.Series | np.ndarray): True values. + y_pred (pd.Series | np.ndarray): Predicted values. + y_true_upper (pd.Series | np.ndarray | None): Optional upper bounds for true + values. Defaults to ``y_true`` if omitted. + y_true_lower (pd.Series | np.ndarray | None): Optional lower bounds for true + values. Defaults to ``y_true`` if omitted. + confidence_interval (float | pd.Series | np.ndarray | None): Optional + confidence interval (full width) for true values. + + Returns: + float: The relative absolute error (RAE) score, with inverse confidence + weighting applied to both the model error and the naive-baseline error. + + Raises: + ValueError: If both explicit bounds and a confidence interval are provided. + + """ + y_true_lower, y_true_upper = _resolve_bounds( + y_true, y_true_upper, y_true_lower, confidence_interval + ) + interval_width = np.clip( + y_true_upper - y_true_lower, a_min=_MIN_CONFIDENCE_INTERVAL, a_max=None + ) + weights = 1.0 / interval_width + + abs_error = np.abs(y_true - y_pred) + + # Weighted-least-squares-optimal constant, i.e. the weighted analogue of + # mean(y_true) — the best a naive constant predictor can do under this weighting. + weighted_mean = np.sum(weights * y_true) / np.sum(weights) + baseline_error = np.abs(y_true - weighted_mean) + + return np.sum(weights * abs_error) / np.sum(weights * baseline_error) diff --git a/evaluation/evaluate_predictions.py b/evaluation/evaluate_predictions.py index 56998b7..d968280 100644 --- a/evaluation/evaluate_predictions.py +++ b/evaluation/evaluate_predictions.py @@ -1,4 +1,10 @@ -"""Functions for evaluating the predictions of the OpenADMET CYP blind challenge.""" +"""Functions for evaluating the predictions of the OpenADMET CYP blind challenge. + +Ported from the challenge backend, so it should match exactly what you see on the leaderboard! +""" + +import inspect +from typing import Callable import numpy as np import pandas as pd @@ -7,10 +13,13 @@ from .config import ( ACTIVITY_METRICS, BOOTSTRAP_SAMPLES, - ENDPOINTS, + CLASSIFICATION_ENDPOINTS, + CLASSIFICATION_METRICS, ENDPOINTS_TO_LOG_TRANSFORM, - TDI_ENDPOINT, - TDI_METRICS, + MACRO_ENDPOINT_LABEL, + METRIC_NAN_FALLBACK, + REGRESSION_CREDIBLE_INTERVALS_LOWER_SUFFIX, + REGRESSION_CREDIBLE_INTERVALS_UPPER_SUFFIX, ) from .utils import bootstrap_sampling, clip_and_log_transform @@ -20,31 +29,69 @@ # --------------------------------------------------------------------------- +def _metrics_for_endpoint(endpoint: str) -> list[tuple[str, Callable]]: + """Return the metric list to use for a given activity endpoint. + + Classification (TDI) endpoints are scored with ``CLASSIFICATION_METRICS`` + (MCC/Accuracy/Precision/Recall/F1); every other activity endpoint (regression, + direct-inhibition pIC50) is scored with ``ACTIVITY_METRICS``. + """ + return CLASSIFICATION_METRICS if endpoint in CLASSIFICATION_ENDPOINTS else ACTIVITY_METRICS + + def score_activity_predictions( - predictions: pd.DataFrame, ground_truth: pd.DataFrame + predictions: pd.DataFrame, ground_truth: pd.DataFrame, endpoints: list[str] ) -> pd.DataFrame: """Score the activity predictions against the ground truth. Metrics are calculated for bootstrapped samples of the dataset to allow for testing - the statistical significance of differences between submissions. + the statistical significance of differences between submissions. Each endpoint is + scored with the metric list appropriate to its type — regression endpoints get + ``ACTIVITY_METRICS``, classification (TDI) endpoints get ``CLASSIFICATION_METRICS`` + — see ``_metrics_for_endpoint``. + + Each endpoint is scored only on the compounds that have a ground-truth value for + that endpoint — a compound not tested for a given endpoint has ``y_true == NaN`` + there and is excluded from that endpoint's bootstrap sampling entirely. Every + compound with a ground-truth value is expected to have a prediction (participants + are asked to predict every compound, and submission_validation.py is the primary + check for missing predictions); a NaN prediction for such a compound is treated as + a validation failure here too. + + Regression endpoints carry credible-interval bound columns in ``ground_truth`` + (named ``f"{endpoint}{REGRESSION_CREDIBLE_INTERVALS_UPPER_SUFFIX}"`` / + ``f"{endpoint}{REGRESSION_CREDIBLE_INTERVALS_LOWER_SUFFIX}"``), used by the + soft-thresholded RAE metric (``ST-RAE``). When present, these are threaded through + to ``bootstrap_metrics`` alongside ``y_true``/``y_pred``; classification endpoints + have no such columns, so ``None`` is passed instead (harmless, since none of + ``CLASSIFICATION_METRICS`` consume them). + + This function does not compute the macro-averaged "MA" pseudo-endpoint. Callers + that want a track's "MA" row should call ``add_macro_endpoint`` on this function's + output with the same ``endpoints`` (and that track's own metrics). Args: predictions (pd.DataFrame): The predicted activity values. ground_truth (pd.DataFrame): The true activity values. + endpoints (list[str]): The endpoints to score, e.g. ``REGRESSION_ENDPOINTS`` + or ``CLASSIFICATION_ENDPOINTS`` — regression and classification are + independent submission tracks, so a given call only ever scores one + track's endpoints. Returns: pd.DataFrame: A DataFrame containing the scored bootstrapped activity - predictions. + predictions, one row per (endpoint, bootstrap sample) — no macro + pseudo-endpoint included. Raises: - ValueError: If the merged DataFrame contains NaN values after merging - predictions with ground truth. + ValueError: If a compound with a ground-truth value for an endpoint has no + prediction for that endpoint. """ logger.info("Scoring activity predictions against ground truth") merged_df = predictions.merge( - ground_truth, on="Molecule Name", suffixes=("_pred", "_true"), how="right" - ).sort_values("Molecule Name") + ground_truth, on="Molecule_Name", suffixes=("_pred", "_true"), how="right" + ).sort_values("Molecule_Name") logger.info( "Completed merging predictions with ground truth. Merged dataset contains {} " "rows and {} columns.", @@ -52,38 +99,218 @@ def score_activity_predictions( merged_df.shape[1], ) - if merged_df.isnull().any().any(): - logger.warning( - "Merged DataFrame contains NaN values after merging predictions with ground" - " truth. This may indicate missing predictions for some molecules." - ) - raise ValueError( - "Merged DataFrame contains NaN values after merging predictions with ground truth." - ) - all_endpoint_bootstrap_results_list = [] - for endpoint in ENDPOINTS: + for endpoint in endpoints: logger.info("Scoring endpoint: {}", endpoint) y_pred = merged_df[f"{endpoint}_pred"].to_numpy() y_true = merged_df[f"{endpoint}_true"].to_numpy() - if endpoint in ENDPOINTS_TO_LOG_TRANSFORM: + # Credible-interval bound columns aren't merge-suffixed: they only ever come + # from ground_truth (predictions never carry them), so they keep their plain + # names — see merge() above (suffixes only apply to overlapping columns). + upper_col = f"{endpoint}{REGRESSION_CREDIBLE_INTERVALS_UPPER_SUFFIX}" + lower_col = f"{endpoint}{REGRESSION_CREDIBLE_INTERVALS_LOWER_SUFFIX}" + y_true_upper = ( + merged_df[upper_col].to_numpy() if upper_col in merged_df.columns else None + ) + y_true_lower = ( + merged_df[lower_col].to_numpy() if lower_col in merged_df.columns else None + ) + + # pd.isna (not np.isnan) so this works for classification endpoints too — + # their ground-truth column can be object/bool dtype (e.g. CYP2D6_is_TDI has + # a couple of genuinely missing labels), which np.isnan can't handle. + has_ground_truth = ~pd.isna(y_true) + if not has_ground_truth.all(): + logger.debug( + "Excluding {} compound(s) with no ground truth for endpoint {}", + (~has_ground_truth).sum(), + endpoint, + ) + y_pred = y_pred[has_ground_truth] + y_true = y_true[has_ground_truth] + if y_true_upper is not None: + y_true_upper = y_true_upper[has_ground_truth] + if y_true_lower is not None: + y_true_lower = y_true_lower[has_ground_truth] + + # A submission itself must never contain NaN (submission_validation.py's + # nullable=False already rejects that) — this instead defends against a + # prediction going missing specifically for a compound that *does* have + # ground truth, which validation of the raw submission can't catch on its own. + if pd.isna(y_pred).any(): + raise ValueError( + f"Missing prediction(s) for endpoint '{endpoint}': every compound " + "with a ground-truth value must have a prediction." + ) + + if endpoint in CLASSIFICATION_ENDPOINTS: + # Safe only after NaN rows have already been dropped from both arrays. + y_true = y_true.astype(bool) + y_pred = y_pred.astype(bool) + elif endpoint in ENDPOINTS_TO_LOG_TRANSFORM: logger.debug("Applying log transformation to endpoint {}", endpoint) y_pred = clip_and_log_transform(y_pred) y_true = clip_and_log_transform(y_true) bootstrap_df = bootstrap_metrics( - y_pred, y_true, endpoint, ACTIVITY_METRICS, n_bootstrap_samples=BOOTSTRAP_SAMPLES + y_pred, + y_true, + endpoint, + n_bootstrap_samples=BOOTSTRAP_SAMPLES, + metrics=_metrics_for_endpoint(endpoint), + y_true_upper=y_true_upper, + y_true_lower=y_true_lower, ) all_endpoint_bootstrap_results_list.append(bootstrap_df) all_endpoint_bootstrap_results = pd.concat( all_endpoint_bootstrap_results_list, ignore_index=True ) - all_endpoint_bootstrap_results = all_endpoint_bootstrap_results.fillna(0) logger.info("Completed scoring activity predictions") return all_endpoint_bootstrap_results +def add_macro_endpoint( + all_endpoint_bootstrap_results: pd.DataFrame, + endpoints: list[str], + metrics: list[tuple[str, Callable]], +) -> pd.DataFrame: + """Narrow to one track's endpoints/metrics and append its macro "MA" row. + + ``score_activity_predictions`` scores every activity endpoint (regression and + classification) in one call, concatenating per-endpoint frames that have + *different* metric columns (regression rows have MAE/ST-RAE/..., classification rows + have MCC/Accuracy/...) — the concatenated result has both sets of columns, NaN + wherever a metric doesn't apply to that row's endpoint. This filters rows down to + just ``endpoints`` (one track's real endpoints) and columns down to just + ``metrics`` (that track's own metrics), so no cross-track NaN columns leak into + the result, then appends a macro-averaged "MA" row set (via + ``compute_macro_bootstrap_results``) when there's more than one endpoint to + average across. + + Args: + all_endpoint_bootstrap_results (pd.DataFrame): Output of + ``score_activity_predictions`` (or any frame with "Sample", "Endpoint", + and metric columns for multiple endpoints/tracks). + endpoints (list[str]): The track's real endpoints to keep, e.g. + ``REGRESSION_ENDPOINTS`` or ``CLASSIFICATION_ENDPOINTS``. May be empty, in + which case the result is empty (callers should generally avoid calling + this with an empty list rather than relying on that). + metrics (list[tuple[str, Callable]]): The track's own metric list, e.g. + ``ACTIVITY_METRICS`` or ``CLASSIFICATION_METRICS`` — only these columns + are kept. + + Returns: + pd.DataFrame: This track's real-endpoint rows, narrowed to its own metric + columns, plus a macro "MA" row set when ``len(endpoints) > 1``. + + """ + metric_names = [name for name, _ in metrics] + track_results = all_endpoint_bootstrap_results[ + all_endpoint_bootstrap_results["Endpoint"].isin(endpoints) + ][["Sample", "Endpoint", *metric_names]] + + if len(endpoints) > 1: + logger.info( + "Calculating macro-averaged metrics across endpoints for each bootstrap sample" + ) + macro_bootstrap_results = compute_macro_bootstrap_results( + track_results, metrics=metrics + ) + track_results = pd.concat( + [track_results, macro_bootstrap_results], ignore_index=True + ) + return track_results + + +def compute_macro_bootstrap_results( + all_endpoint_bootstrap_results: pd.DataFrame, + metrics: list[tuple[str, Callable]], +) -> pd.DataFrame: + """Compute per-bootstrap-sample macro-averaged metrics across all endpoints. + + For every bootstrap sample, every metric in ``metrics`` is macro-averaged across + endpoints with a plain arithmetic mean. + + Spearman_R was previously averaged via a Fisher z-transform (``arctanh`` / + ``tanh``), the standard variance-stabilising treatment for combining several + noisy *estimates of the same underlying correlation* (e.g. meta-analysis, or + averaging one endpoint's Spearman across repeated resamples of the same data). + That doesn't apply here: this average combines Spearman scores from *different* + endpoints (different isoforms), which are unrelated true correlations, not + repeated estimates of one. Fisher's z blows up near +/-1 (``arctanh(1) = inf``, + clipped in practice but still huge — e.g. ``arctanh(1 - 1e-7) ≈ 8.4`` vs. + ``arctanh(0) = 0``), so a submission with a near-perfect Spearman on a couple of + endpoints and ~0 on the rest could macro-average to ~0.99 instead of the + naively-expected ~0.5, letting a handful of easy/lucky endpoints dominate the + macro score. A plain mean — already used for ST-RAE/MAE/R2/Kendall_Tau — doesn't have + this failure mode, so Spearman_R now uses one too, for the same reason + Kendall_Tau always has: Fisher's z-transform has no standard extension to + Kendall's tau (different asymptotic sampling distribution), so there was never a + transform-based option for it here. + + Args: + all_endpoint_bootstrap_results (pd.DataFrame): Per-endpoint bootstrap metrics + for a single track, as returned by ``add_macro_endpoint``'s narrowing step + (or by concatenating per-endpoint ``bootstrap_metrics(...)`` results). + Must contain "Sample", "Endpoint", and one column per metric in + ``metrics``. + metrics (list[tuple[str, Callable]]): The metric list to macro-average — only + the names are used here (e.g. ``ACTIVITY_METRICS`` or + ``CLASSIFICATION_METRICS``). + + Returns: + pd.DataFrame: One row per bootstrap sample, with columns "Sample", + "Endpoint" (``MACRO_ENDPOINT_LABEL`` for every row), and the macro-averaged + value of each metric in ``metrics`` for that sample. + + """ + grouped = all_endpoint_bootstrap_results.groupby("Sample") + macro_results = pd.DataFrame(index=grouped.size().index) + for metric_name, _ in metrics: + logger.info( + "Computing macro-averaged metric {} across bootstrap iterations", + metric_name, + ) + macro_results[metric_name] = grouped[metric_name].mean() + macro_results = macro_results.reset_index() + macro_results["Endpoint"] = MACRO_ENDPOINT_LABEL + return macro_results + + +def pivot_endpoint_results_wide(by_endpoint_results: pd.DataFrame) -> pd.DataFrame: + """Pivot per-endpoint mean/std results into a single wide row. + + ``by_endpoint_results`` (as returned by ``average_bootstrap_results_by_endpoint``) + has one row per endpoint (indexed by endpoint name, including the synthetic + ``MACRO_ENDPOINT_LABEL`` ("MA") pseudo-endpoint computed by + ``compute_macro_bootstrap_results``) and one column per ``_mean`` / + ``_std``. This flattens it into a single-row DataFrame suitable for saving + as a submission's ``averaged-results.parquet``, with every endpoint's columns + consistently prefixed as ``f"{endpoint}_{metric}_{mean|std}"`` (e.g. + ``"CYP3A4_pIC50_active_MAE_mean"``, ``"MA_ST-RAE_mean"``). + + A leaderboard for any single endpoint (macro or real) is then built by narrowing + back down to that endpoint's columns and stripping the prefix, so ``primary_metric`` + is always a bare metric name (e.g. ``"ST-RAE"``) regardless of which endpoint a + given leaderboard targets. + + Args: + by_endpoint_results (pd.DataFrame): Per-endpoint mean/std results, indexed by + endpoint name (including ``MACRO_ENDPOINT_LABEL``). + + Returns: + pd.DataFrame: A single-row DataFrame with one column per + endpoint/metric/statistic combination. + + """ + wide_row: dict[str, float] = {} + for endpoint, row in by_endpoint_results.iterrows(): + for column, value in row.items(): + wide_row[f"{endpoint}_{column}"] = value + return pd.DataFrame([wide_row]) + + def average_bootstrap_results_by_endpoint( all_endpoint_bootstrap_results: pd.DataFrame, ) -> pd.DataFrame: @@ -108,12 +335,26 @@ def average_bootstrap_results_by_endpoint( return agg_df +def _metric_needs_credible_interval_bounds(metric_func: Callable) -> bool: + """True if ``metric_func`` accepts ``y_true_upper``/``y_true_lower`` keywords. + + Lets ``bootstrap_metrics`` dispatch the credible-interval bounds only to metrics + that use them (e.g. ``rae_soft_threshold_absolute_error``), while other metrics in + the same list (MAE, R2, ...) keep their plain two-argument call — introspecting + the signature avoids hardcoding metric names here. + """ + params = inspect.signature(metric_func).parameters + return "y_true_upper" in params and "y_true_lower" in params + + def bootstrap_metrics( y_pred: np.ndarray, y_true: np.ndarray, endpoint: str, - metrics: list[tuple[str, object]], n_bootstrap_samples: int, + metrics: list[tuple[str, Callable]] = ACTIVITY_METRICS, + y_true_upper: np.ndarray | None = None, + y_true_lower: np.ndarray | None = None, ) -> pd.DataFrame: """Calculate bootstrap metrics given predicted and true values. @@ -121,91 +362,80 @@ def bootstrap_metrics( y_pred (np.ndarray): The predicted values. y_true (np.ndarray): The true values. endpoint (str): The endpoint for which the metrics are being calculated. - metrics (list[tuple[str, object]]): List of (metric_name, metric_func) pairs. n_bootstrap_samples (int): The number of bootstrap samples to generate. + metrics (list[tuple[str, Callable]]): The ``(name, func)`` metric list to + compute for every bootstrap sample — ``ACTIVITY_METRICS`` for a + regression endpoint, ``CLASSIFICATION_METRICS`` for a classification + (TDI) endpoint. Defaults to ``ACTIVITY_METRICS``. + y_true_upper (np.ndarray | None): Per-compound upper credible-interval bound + for ``y_true``, aligned with ``y_true``/``y_pred``. Only consumed by + metrics whose signature accepts ``y_true_upper``/``y_true_lower`` (see + ``_metric_needs_credible_interval_bounds``), e.g. the soft-thresholded + RAE metric — ignored by every other metric. Required if ``metrics`` + includes such a metric, otherwise optional. + y_true_lower (np.ndarray | None): Per-compound lower credible-interval bound, + counterpart to ``y_true_upper``. Returns: pd.DataFrame: A DataFrame containing the bootstrap metrics for the given endpoint. + Raises: + RuntimeError: If a metric cannot be calculated, or returns a non-finite + value with no entry in ``METRIC_NAN_FALLBACK``, for any bootstrap sample + — rather than silently scoring that sample as 0 (which would misrepresent + a real failure as a perfect score for error metrics like MAE/ST-RAE). + Metrics listed in ``METRIC_NAN_FALLBACK`` (e.g. Spearman_R/Kendall_Tau, + which are mathematically undefined for a zero-variance bootstrap sample — + such as a submission predicting the same value for every compound) use + that fallback value instead of raising. This also covers a metric that + needs credible-interval bounds (e.g. ST-RAE) when none were supplied. + """ + metrics_with_bounds_flag = [ + (name, func, _metric_needs_credible_interval_bounds(func)) + for name, func in metrics + ] + bootstrap_metrics_list = [] for bootstrap_iteration, idx in enumerate( bootstrap_sampling(y_true.shape[0], n_bootstrap_samples) ): metric_values = {"Sample": bootstrap_iteration, "Endpoint": endpoint} - for metric_name, metric_func in metrics: + for metric_name, metric_func, needs_bounds in metrics_with_bounds_flag: try: - metric_value = metric_func(y_true[idx], y_pred[idx]) + if needs_bounds: + if y_true_upper is None or y_true_lower is None: + raise ValueError( + f"Metric '{metric_name}' requires credible-interval " + "bounds (y_true_upper/y_true_lower), but none were " + "provided to bootstrap_metrics." + ) + metric_value = metric_func( + y_true[idx], + y_pred[idx], + y_true_upper=y_true_upper[idx], + y_true_lower=y_true_lower[idx], + ) + else: + metric_value = metric_func(y_true[idx], y_pred[idx]) + if not isinstance(metric_value, (int, float)): + metric_value = metric_value.statistic except Exception as e: - logger.warning( - f"Error calculating metric {metric_name} for endpoint {endpoint}: {e}" - ) - metric_value = np.nan - if not isinstance(metric_value, (int, float)): - metric_value = metric_func(y_true[idx], y_pred[idx]).statistic + raise RuntimeError( + f"Error calculating metric '{metric_name}' for endpoint " + f"'{endpoint}' (bootstrap sample {bootstrap_iteration}): {e}" + ) from e + if not np.isfinite(metric_value): + if metric_name not in METRIC_NAN_FALLBACK: + raise RuntimeError( + f"Metric '{metric_name}' for endpoint '{endpoint}' " + f"(bootstrap sample {bootstrap_iteration}) returned a " + f"non-finite value: {metric_value}" + ) + metric_value = METRIC_NAN_FALLBACK[metric_name] metric_values[metric_name] = metric_value bootstrap_metrics_list.append(metric_values) bootstrap_df = pd.DataFrame(bootstrap_metrics_list) return bootstrap_df - - -# --------------------------------------------------------------------------- -# TDI (time-dependent inhibition) classification scoring -# --------------------------------------------------------------------------- - - -def score_tdi_predictions( - predictions: pd.DataFrame, ground_truth: pd.DataFrame -) -> pd.DataFrame: - """Score TDI classification predictions against the ground truth. - - ``predictions`` is expected to contain a ``TDI_probability`` column (the - predicted probability that a compound is TDI-positive); ``ground_truth`` - is expected to contain the binary ``is_TDI`` label. Metrics are computed - on bootstrapped samples of the dataset, mirroring the activity scoring - pipeline. - - Args: - predictions (pd.DataFrame): The predicted TDI probabilities. - ground_truth (pd.DataFrame): The true TDI labels. - - Returns: - pd.DataFrame: A DataFrame containing the scored bootstrapped TDI - predictions. - - Raises: - ValueError: If the merged DataFrame contains NaN values after merging - predictions with ground truth. - - """ - logger.info("Scoring TDI predictions against ground truth") - merged_df = predictions.merge( - ground_truth, on="Molecule Name", suffixes=("_pred", "_true"), how="right" - ).sort_values("Molecule Name") - logger.info( - "Completed merging predictions with ground truth. Merged dataset contains {} " - "rows and {} columns.", - merged_df.shape[0], - merged_df.shape[1], - ) - - if merged_df.isnull().any().any(): - logger.warning( - "Merged DataFrame contains NaN values after merging predictions with ground" - " truth. This may indicate missing predictions for some molecules." - ) - raise ValueError( - "Merged DataFrame contains NaN values after merging predictions with ground truth." - ) - - y_pred = merged_df["TDI_probability"].to_numpy() - y_true = merged_df[TDI_ENDPOINT].to_numpy() - - bootstrap_df = bootstrap_metrics( - y_pred, y_true, TDI_ENDPOINT, TDI_METRICS, n_bootstrap_samples=BOOTSTRAP_SAMPLES - ) - bootstrap_df = bootstrap_df.fillna(0) - logger.info("Completed scoring TDI predictions") - return bootstrap_df