Fraction of best ligands

cinnabar/__init__.py

@@ -10,4 +10,5 @@
 from cinnabar import stats
 from cinnabar.femap import FEMap, unit
 from cinnabar.measurements import Measurement, ReferenceState
+from cinnabar.classification_metrics import compute_fraction_best_ligands
 # from cinnabar. import plotting

cinnabar/classification_metrics.py

@@ -0,0 +1,130 @@
+# This code is part of cinnabar and is licensed under the MIT license.
+# For details, see https://github.com/OpenFreeEnergy/cinnabar
+
+from typing import Iterable
+import numpy as np
+import math
+from numpy.typing import NDArray
+
+
+def _create_2d_histogram(y_true: Iterable[float], y_pred: Iterable[float]) -> tuple[NDArray, NDArray, NDArray]:
+    """
+    Create a 2D histogram from two arrays of data.
+
+    Parameters
+    ----------
+    y_true : array-like
+        The true values.
+    y_pred : array-like
+        The predicted values.
+
+    Returns
+    -------
+    histogram : ndarray
+        The 2D histogram of the input data.
+    bins_true : ndarray
+        The bin edges along the y_true axis.
+    bins_pred : ndarray
+        The bin edges along the y_pred axis.
+
+    Raises
+    ------
+    ValueError
+        If `y_true` and `y_pred` have different lengths.
+    """
+
+    y_true = np.asarray(y_true)
+    y_pred = np.asarray(y_pred)
+
+    if y_true.shape != y_pred.shape:
+        raise ValueError("y_true and y_pred must have the same length.")
+
+    y_true_sorted = np.sort(y_true)
+    y_pred_sorted = np.sort(y_pred)
+    # Calculate bin edges using midpoints between sorted values
+    bins_true = np.concatenate(([y_true.min()], (y_true_sorted[:-1] + y_true_sorted[1:]) / 2, [y_true.max()]))
+    bins_pred = np.concatenate(([y_pred.min()], (y_pred_sorted[:-1] + y_pred_sorted[1:]) / 2, [y_pred.max()]))
+
+    # Note a perfect prediction will have all counts in the diagonal bins
+    histogram, bins_true, bins_pred = np.histogram2d(y_true, y_pred, bins=[bins_true, bins_pred])
+
+    return histogram, bins_true, bins_pred
+
+
+def _compute_overlap_coefficient(histogram: NDArray, ranking: int) -> float:
+    """
+    Compute the overlap coefficient from a 2D histogram.
+
+    The overlap coefficient is calculated based on the counts in the histogram
+    for the top N ranked ligands (most active).
+
+    Parameters
+    ----------
+    histogram : ndarray
+        A 2D histogram array where the counts are stored.
+    ranking : int
+        The number of rankings to consider when computing overlap.
+
+    Returns
+    -------
+    float
+        The overlap coefficient.
+
+    Raises
+    ------
+    ValueError
+        If `top_n_ligands` is greater than the number of ligands in the histogram.
+    """
+    if ranking < 1:
+        raise ValueError("Ranking must be greater than 0.")
+
+    if histogram.shape[0] < ranking:
+        raise ValueError("Ranking must be less than the number of ligands.")
+
+    overlap = np.sum(histogram[:ranking, :ranking])
+
+    return overlap / ranking
+
+
+def compute_fraction_best_ligands(y_true: Iterable[float], y_pred: Iterable[float], fraction=0.5) -> float:
+    """
+    Compute the fraction of the best ligands metric introduced by Chris Bayly.
+
+    This function calculates the fraction of the best ligands by computing overlap
+    coefficients for each ranking up to the number of ligands and then averaging up to the specified fraction.
+
+    Parameters
+    ----------
+    y_true : array-like
+        The true values.
+    y_pred : array-like
+        The predicted values.
+    fraction : float, optional
+        The fraction of ligands to consider as the best (default is 0.5).
+
+    Returns
+    -------
+    float
+        The computed fraction of the best ligands.
+
+    Raises
+    ------
+    ValueError
+        If `fraction` is not between 0 and 1.
+    """
+
+    if not (0 <= fraction <= 1):
+        raise ValueError("Fraction must be between 0 and 1.")
+
+    histogram = _create_2d_histogram(y_true, y_pred)[0]
+
+    num_ligands = histogram.shape[0]
+    num_best_ligands = math.floor(num_ligands * fraction)
+
+    overlap_coefficients = [_compute_overlap_coefficient(histogram, i + 1) for i in range(num_ligands)]
+    best_coefficients = overlap_coefficients[:num_best_ligands]
+
+
+    fraction_best_ligands = sum(best_coefficients) / num_best_ligands
+
+    return fraction_best_ligands

cinnabar/tests/test_classification_metrics.py

@@ -0,0 +1,94 @@
+from cinnabar.classification_metrics import _compute_overlap_coefficient, _create_2d_histogram,compute_fraction_best_ligands
+import numpy as np
+import pytest
+
+def test_2d_histogram_wrong_shape():
+    with pytest.raises(ValueError, match="same length"):
+        _create_2d_histogram([1, 2, 3], [1, 2])
+
+
+def test_2d_histogram_simple():
+    # test creating a histogram with simple data
+    # the perfect case with a diagonal histogram
+    hist, xedges, yedges = _create_2d_histogram([1, 2, 3], [1, 2, 3])
+    assert hist.shape == (3, 3)
+    assert np.all(hist.diagonal() == 1)
+    # there should only be 3 non-zero entries
+    assert np.sum(hist) == 3
+    # check bin edges
+    ref_edges = [1, 1.5, 2.5, 3]
+    assert np.allclose(xedges, ref_edges)
+    assert np.allclose(yedges, ref_edges)
+
+    # test a case with non-diagonal histogram
+    hist, xedges, yedges = _create_2d_histogram([1, 3.2, 3], [1, 2, 3])
+    assert hist.shape == (3, 3)
+    # there should only be 3 non-zero entries
+    assert np.sum(hist) == 3
+    assert hist[0, 0] == 1
+    assert hist[1, 2] == 1
+    assert hist[2, 1] == 1
+
+
+def test_overlap_coefficient_wrong_ranking():
+    hist = np.array([[1, 0], [0, 1]])
+    with pytest.raises(ValueError, match="greater than 0"):
+        _compute_overlap_coefficient(hist, 0)
+    with pytest.raises(ValueError, match="less than the number of ligands"):
+        _compute_overlap_coefficient(hist, 3)
+
+
+def test_overlap_coefficient_simple():
+    # perfect overlap
+    hist = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]])
+    overlap = _compute_overlap_coefficient(hist, 3)
+    assert overlap == 1.0
+
+    # half overlap
+    hist = np.array([[1, 0, 0], [0, 0, 1], [0, 1, 0]])
+    overlap = _compute_overlap_coefficient(hist, 2)
+    assert overlap == 0.5
+
+    # no overlap
+    hist = np.array([[0, 0, 1], [1, 0, 0], [0, 1, 0]])
+    overlap = _compute_overlap_coefficient(hist, 1)
+    assert overlap == 0.0
+
+
+def test_fraction_best_ligands_bad_fraction():
+    with pytest.raises(ValueError, match="between 0 and 1"):
+        compute_fraction_best_ligands([1, 2, 3], [1, 2, 3], fraction=1.5)
+
+
+def test_fraction_best_ligands_simple():
+    # perfect prediction
+    fraction = compute_fraction_best_ligands([1, 2, 3, 4], [1, 2, 3, 4], fraction=0.5)
+    assert fraction == 1.0
+
+    # 50% correct prediction
+    fraction = compute_fraction_best_ligands([1, 2, 3, 4], [2.5, 2, 4, 3], fraction=0.5)
+    assert fraction == 0.5
+
+    # 75% correct prediction
+    fraction = compute_fraction_best_ligands([1, 2, 3, 4], [1, 4, 2, 3], fraction=0.5)
+    assert fraction == 0.75
+
+    # no correct prediction
+    fraction = compute_fraction_best_ligands([1, 2, 3, 4], [4, 3, 2, 1], fraction=0.5)
+    assert fraction == 0.0
+
+def test_fraction_best_ligands_regression(fe_map):
+    # regression test for a real dataset
+    fe_map.generate_absolute_values()
+
+    # we need to compare the absolute experimental and calculated values so generate them
+    abs_dataframe = fe_map.get_absolute_dataframe()
+    exp_df = abs_dataframe[~abs_dataframe['computational']]
+    calc_df = abs_dataframe[abs_dataframe['computational']]
+
+    # get the calculated and experimental values in the same order
+    merged = exp_df.merge(calc_df, on='label', suffixes=('_exp', '_calc'))
+    y_true = merged['DG (kcal/mol)_exp'].values
+    y_pred = merged['DG (kcal/mol)_calc'].values
+    fraction = compute_fraction_best_ligands(y_true, y_pred, fraction=0.5)
+    assert fraction == pytest.approx(0.7216416707838275)

docs/api.rst

@@ -4,6 +4,7 @@ API Documentation
 .. autosummary::
    :toctree: generated
 
+   cinnabar.classification_metrics
    cinnabar.cli
    cinnabar.femap
    cinnabar.measurements