Fraction of best ligands

cinnabar/__init__.py

@@ -8,6 +8,7 @@
 __version__ = version("cinnabar")
 
 from cinnabar import stats
+from cinnabar.classification_metrics import compute_fraction_best_ligands
 from cinnabar.femap import FEMap, unit
 from cinnabar.measurements import Measurement, ReferenceState
 # from cinnabar. import plotting

cinnabar/_due.py

@@ -78,4 +78,4 @@ def _donothing_func(*args: Any, **kwargs: Any) -> None:
 # py-indent-offset: 4
 # tab-width: 4
 # indent-tabs-mode: nil
-# End:
+# End:

cinnabar/classification_metrics.py

@@ -0,0 +1,129 @@
+# This code is part of cinnabar and is licensed under the MIT license.
+# For details, see https://github.com/OpenFreeEnergy/cinnabar
+
+import math
+from typing import Iterable
+
+import numpy as np
+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: float = 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
+        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_best_ligands)]
+
+    fraction_best_ligands = sum(overlap_coefficients) / num_best_ligands
+
+    return fraction_best_ligands

cinnabar/stats.py

@@ -5,15 +5,16 @@
 import scipy
 import sklearn.metrics
 
-from cinnabar._due import due, Doi
+from cinnabar._due import Doi, due
 
 due.cite(
     Doi("10.1021/acs.jcim.9b00528"),
     description="Compute maximum likelihood estimate of free energies and covariance in their estimates",
     path="cinnabar.stats.mle",
-    cite_module=True
+    cite_module=True,
 )
 
+
 def bootstrap_statistic(
     y_true: np.ndarray,
     y_pred: np.ndarray,

cinnabar/tests/test_classification_metrics.py

@@ -0,0 +1,101 @@
+import numpy as np
+import pytest
+
+from cinnabar.classification_metrics import (
+    _compute_overlap_coefficient,
+    _create_2d_histogram,
+    compute_fraction_best_ligands,
+)
+
+
+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)

cinnabar/tests/test_due.py

@@ -1,11 +1,10 @@
 import pytest
 
-
 from cinnabar._due import due
 
-
 pytest.importorskip("duecredit")
 
+
 def test_duecredit_mle():
     """Make sure duecredit is captured when the stats module is used"""
     mle_key = ("cinnabar.stats.mle", "10.1021/acs.jcim.9b00528")

docs/api.rst

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

docs/concepts/estimators.rst

@@ -69,4 +69,3 @@ References
 ~~~~~~~~~~~
 
 .. [1] Xu, H., 2019. Optimal measurement network of pairwise differences. Journal of Chemical Information and Modeling, 59(11), pp.4720-4728.
-

docs/concepts/femap.rst

@@ -39,5 +39,3 @@ References
 ~~~~~~~~~~~
 
 .. [1] Xu, H., 2019. Optimal measurement network of pairwise differences. Journal of Chemical Information and Modeling, 59(11), pp.4720-4728.
-
-

docs/concepts/plotting.rst

@@ -65,4 +65,3 @@ References
 ~~~~~~~~~~~
 
 .. [1] Hahn, D.F., Bayly, C.I., Boby, M.L., Macdonald, H.E.B., Chodera, J.D., Gapsys, V., Mey, A.S., Mobley, D.L., Benito, L.P., Schindler, C.E. and Tresadern, G., 2022. Best practices for constructing, preparing, and evaluating protein-ligand binding affinity benchmarks [article v1. 0]. Living journal of computational molecular science, 4(1), p.1497.
-

docs/conf.py

@@ -97,9 +97,9 @@
 # documentation.
 
 html_theme = "ofe_sphinx_theme"
-html_theme_options = {
-    "accent_color": "FeelingSpicy",
-}
+# html_theme_options = {
+#     "accent_color": "FeelingSpicy",
+# }
 
 # Add any paths that contain custom static files (such as style sheets) here,
 # relative to this directory. They are copied after the builtin static files,

docs/installation.rst

@@ -37,4 +37,4 @@ If you're a developer, you will likely want to create a local editable installat
 
 .. code-block:: bash
 
-    $ python -m pip install --no-deps -e .
+    $ python -m pip install --no-deps -e .

news/TEMPLATE.rst

@@ -20,4 +20,4 @@
 
 **Security:**
 
-* <news item>
+* <news item>

news/arsenic_update.rst

@@ -20,4 +20,4 @@
 
 **Security:**
 
-* <news item>
+* <news item>

news/fraction_of_best.rst

@@ -0,0 +1,23 @@
+**Added:**
+
+* add the ``compute_fraction_best_ligands`` function to compute the fraction of best ligands metric `PR#164 <https://github.com/OpenFreeEnergy/cinnabar/pull/164>`_.
+
+**Changed:**
+
+* <news item>
+
+**Deprecated:**
+
+* <news item>
+
+**Removed:**
+
+* <news item>
+
+**Fixed:**
+
+* <news item>
+
+**Security:**
+
+* <news item>