diff --git a/cinnabar/__init__.py b/cinnabar/__init__.py index 742d4b71..8b11262e 100644 --- a/cinnabar/__init__.py +++ b/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 diff --git a/cinnabar/_due.py b/cinnabar/_due.py index d10005cc..f26b06d7 100644 --- a/cinnabar/_due.py +++ b/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: \ No newline at end of file +# End: diff --git a/cinnabar/classification_metrics.py b/cinnabar/classification_metrics.py new file mode 100644 index 00000000..d5265339 --- /dev/null +++ b/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 diff --git a/cinnabar/stats.py b/cinnabar/stats.py index 292f7a08..7a864c3c 100644 --- a/cinnabar/stats.py +++ b/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, diff --git a/cinnabar/tests/test_classification_metrics.py b/cinnabar/tests/test_classification_metrics.py new file mode 100644 index 00000000..178ce96a --- /dev/null +++ b/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) diff --git a/cinnabar/tests/test_due.py b/cinnabar/tests/test_due.py index dfd8a844..cc0c3f4e 100644 --- a/cinnabar/tests/test_due.py +++ b/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") diff --git a/docs/api.rst b/docs/api.rst index 7bc1a046..94c1d4f2 100644 --- a/docs/api.rst +++ b/docs/api.rst @@ -4,6 +4,7 @@ API Documentation .. autosummary:: :toctree: generated + cinnabar.classification_metrics cinnabar.cli cinnabar.femap cinnabar.measurements diff --git a/docs/concepts/estimators.rst b/docs/concepts/estimators.rst index a42efaac..be1aa1bb 100644 --- a/docs/concepts/estimators.rst +++ b/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. - diff --git a/docs/concepts/femap.rst b/docs/concepts/femap.rst index ac2f5b78..e160de82 100644 --- a/docs/concepts/femap.rst +++ b/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. - - diff --git a/docs/concepts/plotting.rst b/docs/concepts/plotting.rst index f9141df2..571c44cd 100644 --- a/docs/concepts/plotting.rst +++ b/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. - diff --git a/docs/conf.py b/docs/conf.py index 50f06ec6..56aabfba 100644 --- a/docs/conf.py +++ b/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, diff --git a/docs/installation.rst b/docs/installation.rst index 96943caf..ce2a9b43 100644 --- a/docs/installation.rst +++ b/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 . \ No newline at end of file + $ python -m pip install --no-deps -e . diff --git a/news/TEMPLATE.rst b/news/TEMPLATE.rst index 9a2b78ad..790d30b1 100644 --- a/news/TEMPLATE.rst +++ b/news/TEMPLATE.rst @@ -20,4 +20,4 @@ **Security:** -* \ No newline at end of file +* diff --git a/news/arsenic_update.rst b/news/arsenic_update.rst index cfe9f434..f1ea2172 100644 --- a/news/arsenic_update.rst +++ b/news/arsenic_update.rst @@ -20,4 +20,4 @@ **Security:** -* \ No newline at end of file +* diff --git a/news/fraction_of_best.rst b/news/fraction_of_best.rst new file mode 100644 index 00000000..7c4c8655 --- /dev/null +++ b/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 `_. + +**Changed:** + +* + +**Deprecated:** + +* + +**Removed:** + +* + +**Fixed:** + +* + +**Security:** + +*