Calibrate against several LC setups at once (v4.3.0)

CHANGELOG.md

@@ -6,6 +6,36 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to
 [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## [4.3.0] - 2026-08-28
+
+### Added
+
+- `MultiHeadRidgeCalibration`, a calibration that uses several LC setups of the multitask model
+  instead of only the best-correlating one. It keeps the existing ranking, calibrates the 80 best
+  heads individually with `SplineTransformerCalibration`, and fits a ridge regression from those
+  calibrated estimates onto the observed retention times.
+
+  On the eight PRIDE setups that no DeepLC model was trained on it lowered the held-out error on
+  all eight, by a median of 13 % relative to the observed gradient (0.01248 to 0.01090 MAE/span):
+  0.7 % on a setup that pools several fractions into one run, 7.5 to 17 % on four others, 28 to
+  38 % on the three where the single best head fitted worst, the largest gain being the smallest
+  reference of 230 peptidoforms. Fitting is not slower than the current path despite the extra
+  spline fits, because the head ranking is vectorised (a median of 1.0 s against 2.3 s in this
+  test), and prediction is unchanged since the full head matrix is computed either way.
+
+  Opt in per call, the default is untouched:
+
+  ```python
+  from deeplc import MultiHeadRidgeCalibration, predict_and_calibrate
+
+  rt = predict_and_calibrate(psms, psm_list_reference=reference,
+                             calibration=MultiHeadRidgeCalibration())
+  ```
+
+- `Calibration.uses_all_heads`, telling `calibrate` and `predict_and_calibrate` to hand a
+  calibration the whole `(n, n_heads)` prediction matrix rather than a single column. False for
+  every calibration that works on one series, which is all of them except the above.
+
 ## [4.2.0] - 2026-08-28
 
 ### Changed

deeplc/__init__.py

@@ -2,6 +2,7 @@
 
 from importlib.metadata import version
 
+from deeplc.calibration import MultiHeadRidgeCalibration
 from deeplc.core import (
     calibrate,
     finetune,
@@ -14,6 +15,7 @@
 
 __version__: str = version("deeplc")
 __all__: list[str] = [
+    "MultiHeadRidgeCalibration",
     "calibrate",
     "predict",
     "predict_and_calibrate",

deeplc/calibration.py

@@ -21,6 +21,11 @@ class Calibration(ABC):
 
     selected_model_head: int | None = None
 
+    #: Whether ``fit`` and ``transform`` take the whole ``(n, n_heads)`` prediction matrix of a
+    #: multitask model instead of a single series. False for every calibration here except
+    #: :class:`MultiHeadRidgeCalibration`.
+    uses_all_heads: bool = False
+
     @abstractmethod
     def __init__(self, *args, **kwargs):
         """Initialize the calibration model."""
@@ -325,6 +330,160 @@ def transform(self, source: np.ndarray) -> np.ndarray:
         return np.array(cal_preds)
 
 
+class MultiHeadRidgeCalibration(Calibration):
+    """
+    Calibrate against several LC setups of a multitask model at once.
+
+    The default path keeps one setup head: it ranks all heads by Pearson correlation to the
+    reference and fits a spline on the winner. A gradient that no trained setup matches exactly
+    is then described by the closest single setup, and the rest of the matrix is discarded.
+
+    This calibration keeps that ranking, calibrates the ``n_heads`` best heads individually with
+    :class:`SplineTransformerCalibration`, and fits a ridge regression from those calibrated
+    estimates onto the observed retention times. Every head therefore contributes an estimate
+    already in the unit of the reference, and the ridge decides how much to trust each one.
+
+    On the eight PRIDE setups that no DeepLC model was trained on, this lowered the held-out
+    error on all eight: a median of 13 % relative to the observed gradient (0.01248 to 0.01090
+    MAE/span), from 1 % on a setup whose retention times are not a single gradient to 38 % on the
+    smallest reference of 230 peptidoforms. The cost is ``n_heads`` spline fits and one ridge on
+    the reference; prediction is unchanged, because the full matrix is computed either way.
+
+    Parameters
+    ----------
+    n_heads
+        How many of the best-correlating heads to combine. 80 sits on the flat part of the
+        optimum for references from 230 to 2,000 peptidoforms; below about ten the gain shrinks,
+        and beyond a few hundred it slowly reverses.
+    alphas
+        Ridge strengths offered to the internal cross-validation. The default spans 1e-3 to 1e6,
+        wide enough for the fit to collapse towards an average when the reference is small.
+
+    """
+
+    uses_all_heads = True
+
+    def __init__(self, n_heads: int = 80, alphas: np.ndarray | None = None) -> None:
+        """Initialize MultiHeadRidgeCalibration."""
+        super().__init__()
+        if n_heads < 1:
+            raise ValueError(f"n_heads must be at least 1, got {n_heads}")
+        self.n_heads = n_heads
+        self.alphas = np.logspace(-3, 6, 19) if alphas is None else np.asarray(alphas)
+        self._head_idx: np.ndarray | None = None
+        self._head_calibrations: list[SplineTransformerCalibration] = []
+        self._ridge = None
+
+    @property
+    def is_fitted(self) -> bool:
+        """True once the heads are selected, calibrated and weighted."""
+        return self._head_idx is not None and self._ridge is not None
+
+    def fit(self, target: np.ndarray, source: np.ndarray) -> None:
+        """
+        Select, calibrate and weight the heads.
+
+        Parameters
+        ----------
+        target
+            Observed retention times of the reference, shape ``(n,)``.
+        source
+            Reference predictions for every head, shape ``(n, n_heads_total)``. A 1-D array is
+            accepted and treated as a single head, so a single-task model still works.
+
+        """
+        from sklearn.linear_model import RidgeCV
+
+        source = np.asarray(source, dtype=np.float64)
+        if source.ndim == 1:
+            source = source[:, None]
+        target = np.asarray(target, dtype=np.float64).ravel()
+        if source.shape[0] != target.shape[0]:
+            raise CalibrationError(
+                f"source has {source.shape[0]} rows and target {target.shape[0]}"
+            )
+        finite = np.isfinite(target) & np.isfinite(source).all(axis=1)
+        if int(finite.sum()) < 3:
+            raise CalibrationError("Fewer than three reference points with finite values.")
+        source, target = source[finite], target[finite]
+
+        order = _rank_heads_by_correlation(source, target)
+        # never fit more weights than half the reference: a 230-peptide reference cannot support
+        # eighty of them, and the ridge would be extrapolating its own regularisation
+        n_heads = int(min(self.n_heads, source.shape[1], max(1, len(target) // 2)))
+        self._head_idx = order[:n_heads]
+        self.selected_model_head = int(order[0])
+
+        calibrated = np.empty((len(target), n_heads), dtype=np.float64)
+        self._head_calibrations = []
+        for position, head in enumerate(self._head_idx):
+            head_calibration = SplineTransformerCalibration()
+            column = source[:, head].astype(np.float32)
+            head_calibration.fit(target=target.astype(np.float32), source=column)
+            calibrated[:, position] = np.asarray(
+                head_calibration.transform(column), dtype=np.float64
+            )
+            self._head_calibrations.append(head_calibration)
+
+        n_splits = int(min(5, max(2, len(target) // 20)))
+        self._ridge = RidgeCV(alphas=self.alphas, cv=n_splits).fit(calibrated, target)
+        LOGGER.info(
+            "Calibrated on %d of %d heads with ridge strength %.4g; head %d correlates best.",
+            n_heads,
+            source.shape[1],
+            float(getattr(self._ridge, "alpha_", float("nan"))),
+            self.selected_model_head,
+        )
+
+    def transform(self, source: np.ndarray) -> np.ndarray:
+        """
+        Calibrate predictions of the model this calibration was fitted with.
+
+        Parameters
+        ----------
+        source
+            Predictions for every head, shape ``(n, n_heads_total)``, as returned by
+            ``predict(..., return_matrix=True)``.
+
+        """
+        if not self.is_fitted:
+            raise CalibrationError("The model has not been fitted yet. Call fit() first.")
+        source = np.asarray(source, dtype=np.float64)
+        if source.ndim == 1:
+            source = source[:, None]
+        head_idx = cast(np.ndarray, self._head_idx)
+        if source.shape[1] <= int(head_idx.max()):
+            raise CalibrationError(
+                f"source has {source.shape[1]} heads, but the calibration was fitted on a model "
+                f"with at least {int(head_idx.max()) + 1}."
+            )
+        if source.shape[0] == 0:
+            return np.array([])
+        calibrated = np.column_stack(
+            [
+                np.asarray(cal.transform(source[:, head].astype(np.float32)), dtype=np.float64)
+                for cal, head in zip(self._head_calibrations, head_idx, strict=True)
+            ]
+        )
+        return np.asarray(self._ridge.predict(calibrated), dtype=np.float64)
+
+
+def _rank_heads_by_correlation(source: np.ndarray, target: np.ndarray) -> np.ndarray:
+    """
+    Head indices by decreasing Pearson correlation with the target, in one pass.
+
+    The same criterion as :func:`deeplc.core._best_correlating_head`, which takes the first
+    element of this order, but vectorised because thousands of heads are ranked at once.
+    """
+    centred = source - source.mean(axis=0)
+    target_centred = target - target.mean()
+    with np.errstate(invalid="ignore", divide="ignore"):
+        denominator = np.sqrt((centred**2).sum(axis=0) * (target_centred**2).sum())
+        correlation = (centred * target_centred[:, None]).sum(axis=0) / denominator
+    correlation = np.where(np.isfinite(correlation), correlation, -np.inf)
+    return np.argsort(-correlation)
+
+
 def _prepare_series(
     target: np.ndarray,
     source: np.ndarray,

deeplc/core.py

@@ -203,6 +203,12 @@ def calibrate(
     LOGGER.debug("Fitting calibration...")
     target_rt_cal = np.array(psm_list_reference["retention_time"], dtype=np.float32)
 
+    # A calibration that combines heads is given the whole matrix and picks its own; it sets
+    # selected_model_head itself, for callers that want to know which setup came out on top.
+    if getattr(calibration, "uses_all_heads", False):
+        calibration.fit(target=target_rt_cal, source=source_rt_cal)
+        return calibration
+
     # Select the best head for calibration if the model predicts for multiple LC setups
     if source_rt_cal.shape[1] > 1:
         calibration.selected_model_head = _best_correlating_head(source_rt_cal, target_rt_cal)
@@ -278,6 +284,10 @@ def predict_and_calibrate(
     else:
         LOGGER.info("Calibration is already fitted, skipping fitting step.")
 
+    if getattr(calibration, "uses_all_heads", False):
+        # the calibration combines several heads, so it takes the matrix as it is
+        return calibration.transform(predicted_rt)
+
     if predicted_rt.shape[1] > 1:
         if calibration.selected_model_head is None:
             raise ValueError(

docs/source/models.rst

@@ -22,6 +22,33 @@ The 4.0 default, ``multitask_model.pt`` (shared trunk, one head per setup), stay
 bundled as :data:`deeplc.core.LEGACY_MULTITASK_MODEL` and can be passed as
 ``model=`` to any core function to reproduce 4.0 and 4.1.0 predictions.
 
+Calibrating against several setups at once
+==========================================
+
+By default calibration keeps one output head: every head is ranked by Pearson correlation to the
+reference and a spline is fitted on the winner. A gradient that no trained setup matches exactly
+is then described by the closest single setup.
+
+:class:`~deeplc.calibration.MultiHeadRidgeCalibration` keeps that ranking but calibrates the 80
+best heads individually and fits a ridge regression from those calibrated estimates onto the
+observed retention times, so several setups contribute:
+
+.. code-block:: python
+
+   from deeplc import MultiHeadRidgeCalibration, predict_and_calibrate
+
+   calibrated_rt = predict_and_calibrate(
+       psm_list,
+       psm_list_reference=reference,
+       calibration=MultiHeadRidgeCalibration(),
+   )
+
+On eight PRIDE setups that no DeepLC model was trained on, this lowered the held-out error on all
+eight, by a median of 13 % relative to the gradient and by up to 38 % on the smallest reference
+(230 peptidoforms). The number of heads is the one parameter worth changing: 80 sits on a flat
+optimum between roughly 40 and 320, and the class never fits more weights than half the reference
+allows. Prediction costs nothing extra, because the full head matrix is computed either way.
+
 Training a model from scratch
 ==============================
 

pyproject.toml

@@ -1,6 +1,6 @@
 [project]
 name = "deeplc"
-version = "4.2.0"
+version = "4.3.0"
 description = "DeepLC: Retention time prediction for (modified) peptides using Deep Learning."
 readme = "README.md"
 license = { file = "LICENSE" }