From b79c77da31e50975c537d72a5c95cc77047f0a3e Mon Sep 17 00:00:00 2001 From: Hugo MacDermott-Opeskin Date: Mon, 10 Feb 2025 15:40:41 +1100 Subject: [PATCH] distopia 0.4.0 compatibility changes (#4734) * start of docs * improve distopia stup * add backend conditionals * fix wrong backend * shims * working except for dihedrals * move to 0.3.0 * bump ci * bump ci * Update package/MDAnalysis/lib/distances.py Co-authored-by: Oliver Beckstein * fix versioning * add version checking * remove erroneous triclinic note * remove unnesecary import * add explanatory notes * change tests in line with discussion * changelog * Apply suggestions from code review Co-authored-by: Oliver Beckstein * add more tests * improve docs * improve docs again * hmmm * add type for calc_bond_distance_ortho() * just stuff answers back in if we used a buffer * add cram for distance array also * fix changelogs * black * imporve mock test * fix typo * add debug printing * add explanatory comment on use of serial backend comparison * add testing * black * fix typo * fix test * fix dih * black * fix dihedrals * fix reference values * change function names for 0.4.0 * black * fix typo * fix wrong versions --------- Co-authored-by: Oliver Beckstein --- .github/actions/setup-deps/action.yaml | 2 +- package/CHANGELOG | 6 +- package/MDAnalysis/lib/_distopia.py | 149 ++++++++++-- package/MDAnalysis/lib/distances.py | 121 ++++++++-- .../MDAnalysisTests/lib/test_distances.py | 218 ++++++++++++------ 5 files changed, 377 insertions(+), 119 deletions(-) diff --git a/.github/actions/setup-deps/action.yaml b/.github/actions/setup-deps/action.yaml index 91dd56b6d7..9b4c804265 100644 --- a/.github/actions/setup-deps/action.yaml +++ b/.github/actions/setup-deps/action.yaml @@ -59,7 +59,7 @@ inputs: dask: default: 'dask' distopia: - default: 'distopia>=0.2.0,<0.3.0' + default: 'distopia>=0.4.0' h5py: default: 'h5py>=2.10' hole2: diff --git a/package/CHANGELOG b/package/CHANGELOG index 7b5ce3b4c8..e0b75cda0a 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -15,7 +15,7 @@ The rules for this file: ------------------------------------------------------------------------------- ??/??/?? IAlibay, ChiahsinChu, RMeli, tanishy7777, talagayev, tylerjereddy, - marinegor + marinegor, hmacdope * 2.9.0 @@ -27,12 +27,14 @@ Fixes the function to prevent shared state. (Issue #4655) Enhancements + * Improve distopia backend support in line with new functionality available + in distopia >= 0.3.1 (PR #4734) * Addition of 'water' token for water selection (Issue #4839) * Enables parallelization for analysis.density.DensityAnalysis (Issue #4677, PR #4729) * Enables parallelization for analysis.contacts.Contacts (Issue #4660) * Enable parallelization for analysis.nucleicacids.NucPairDist (Issue #4670) * Add check and warning for empty (all zero) coordinates in RDKit converter (PR #4824) - * Added `precision` for XYZWriter (Issue #4775, PR #4771) + * Added `precision` for XYZWriter (Issue #4775, PR #4771) Changes diff --git a/package/MDAnalysis/lib/_distopia.py b/package/MDAnalysis/lib/_distopia.py index 297ce4a3b5..a6ec059168 100644 --- a/package/MDAnalysis/lib/_distopia.py +++ b/package/MDAnalysis/lib/_distopia.py @@ -28,48 +28,51 @@ as a selectable backend. """ import warnings +from packaging.version import Version + +MIN_DISTOPIA_VERSION = Version("0.4.0") # check for distopia try: import distopia except ImportError: HAS_DISTOPIA = False + distopia_version = Version("0.0.0") else: HAS_DISTOPIA = True - # check for compatibility: currently needs to be >=0.2.0,<0.3.0 (issue - # #4740) No distopia.__version__ available so we have to do some probing. - needed_funcs = ["calc_bonds_no_box_float", "calc_bonds_ortho_float"] - has_distopia_020 = all([hasattr(distopia, func) for func in needed_funcs]) - if not has_distopia_020: + # check for compatibility: currently needs to be >=0.4.0, + # some versions of `distopia` don't have a version attribute + try: + distopia_version = Version(distopia.__version__) + except AttributeError: + distopia_version = Version("0.0.0") + if distopia_version < MIN_DISTOPIA_VERSION: warnings.warn( - "Install 'distopia>=0.2.0,<0.3.0' to be used with this " - "release of MDAnalysis. Your installed version of " - "distopia >=0.3.0 will NOT be used.", + f"distopia version {distopia_version} is too old; " + f"need at least {MIN_DISTOPIA_VERSION}, Your installed version of " + "distopia will NOT be used.", category=RuntimeWarning, ) - del distopia HAS_DISTOPIA = False -from .c_distances import ( - calc_bond_distance_triclinic as _calc_bond_distance_triclinic_serial, -) import numpy as np def calc_bond_distance_ortho( - coords1, coords2: np.ndarray, box: np.ndarray, results: np.ndarray + coords1: np.ndarray, + coords2: np.ndarray, + box: np.ndarray, + results: np.ndarray, ) -> None: - distopia.calc_bonds_ortho_float(coords1, coords2, box[:3], results=results) - # upcast is currently required, change for 3.0, see #3927 + distopia.distances_ortho(coords1, coords2, box[:3], results=results) def calc_bond_distance( coords1: np.ndarray, coords2: np.ndarray, results: np.ndarray ) -> None: - distopia.calc_bonds_no_box_float(coords1, coords2, results=results) - # upcast is currently required, change for 3.0, see #3927 + distopia.distances_no_box(coords1, coords2, results=results) def calc_bond_distance_triclinic( @@ -78,8 +81,112 @@ def calc_bond_distance_triclinic( box: np.ndarray, results: np.ndarray, ) -> None: - # redirect to serial backend - warnings.warn( - "distopia does not support triclinic boxes, using serial backend instead." + distopia.distances_triclinic(coords1, coords2, box, results=results) + + +def calc_angle( + coords1: np.ndarray, + coords2: np.ndarray, + coords3: np.ndarray, + results: np.ndarray, +) -> None: + distopia.angles_no_box(coords1, coords2, coords3, results=results) + + +def calc_angle_ortho( + coords1: np.ndarray, + coords2: np.ndarray, + coords3: np.ndarray, + box: np.ndarray, + results: np.ndarray, +) -> None: + distopia.angles_ortho(coords1, coords2, coords3, box[:3], results=results) + + +def calc_angle_triclinic( + coords1: np.ndarray, + coords2: np.ndarray, + coords3: np.ndarray, + box: np.ndarray, + results: np.ndarray, +) -> None: + + distopia.angles_triclinic(coords1, coords2, coords3, box, results=results) + + +def calc_dihedral( + coords1: np.ndarray, + coords2: np.ndarray, + coords3: np.ndarray, + coords4: np.ndarray, + results: np.ndarray, +) -> None: + distopia.dihedrals_no_box( + coords1, coords2, coords3, coords4, results=results + ) + + +def calc_dihedral_ortho( + coords1: np.ndarray, + coords2: np.ndarray, + coords3: np.ndarray, + coords4: np.ndarray, + box: np.ndarray, + results: np.ndarray, +) -> None: + distopia.dihedrals_ortho( + coords1, coords2, coords3, coords4, box[:3], results=results ) - _calc_bond_distance_triclinic_serial(coords1, coords2, box, results) + + +def calc_dihedral_triclinic( + coords1: np.ndarray, + coords2: np.ndarray, + coords3: np.ndarray, + coords4: np.ndarray, + box: np.ndarray, + results: np.ndarray, +) -> None: + distopia.dihedrals_triclinic( + coords1, coords2, coords3, coords4, box, results=results + ) + + +def calc_distance_array( + coords1: np.ndarray, coords2: np.ndarray, results: np.ndarray +) -> None: + distopia.distance_array_no_box(coords1, coords2, results=results) + + +def calc_distance_array_ortho( + coords1: np.ndarray, + coords2: np.ndarray, + box: np.ndarray, + results: np.ndarray, +) -> None: + distopia.distance_array_ortho(coords1, coords2, box[:3], results=results) + + +def calc_distance_array_triclinic( + coords1: np.ndarray, + coords2: np.ndarray, + box: np.ndarray, + results: np.ndarray, +) -> None: + distopia.distance_array_triclinic(coords1, coords2, box, results=results) + + +def calc_self_distance_array(coords: np.ndarray, results: np.ndarray) -> None: + distopia.self_distance_array_no_box(coords, results=results) + + +def calc_self_distance_array_ortho( + coords: np.ndarray, box: np.ndarray, results: np.ndarray +) -> None: + distopia.self_distance_array_ortho(coords, box[:3], results=results) + + +def calc_self_distance_array_triclinic( + coords: np.ndarray, box: np.ndarray, results: np.ndarray +) -> None: + distopia.self_distance_array_triclinic(coords, box, results=results) diff --git a/package/MDAnalysis/lib/distances.py b/package/MDAnalysis/lib/distances.py index bbd1e56dd5..5fa35ecc90 100644 --- a/package/MDAnalysis/lib/distances.py +++ b/package/MDAnalysis/lib/distances.py @@ -48,6 +48,7 @@ "OpenMP" :mod:`c_distances_openmp` parallel implementation in C/Cython with OpenMP + "distopia" `_distopia` SIMD-accelerated implementation ========== ========================= ====================================== Use of the distopia library @@ -61,29 +62,33 @@ you wish to use is covered by distopia. For more information see the `distopia documentation`_. -.. Table:: Functions available using the `distopia`_ backend. - :align: center - - +-------------------------------------+-----------------------------------+ - | Functions | Notes | - +=====================================+===================================+ - | MDAnalysis.lib.distances.calc_bonds | Doesn't support triclinic boxes | - +-------------------------------------+-----------------------------------+ +.. table:: Functions available using the `distopia`_ backend. + :align: center + + +-----------------------------------------------+ + | Functions | + +===============================================+ + | MDAnalysis.lib.distances.calc_bonds | + +-----------------------------------------------+ + | MDAnalysis.lib.distances.calc_angles | + +-----------------------------------------------+ + | MDAnalysis.lib.distances.calc_dihedrals | + +-----------------------------------------------+ + | MDAnalysis.lib.distances.distance_array | + +-----------------------------------------------+ + | MDAnalysis.lib.distances.self_distance_array | + +-----------------------------------------------+ If `distopia`_ is installed, the functions in this table will accept the key 'distopia' for the `backend` keyword argument. If the distopia backend is selected the `distopia` library will be used to calculate the distances. Note that for functions listed in this table **distopia is not the default backend -if and must be selected.** +and must be selected.** -.. Note:: - Distopia does not currently support triclinic simulation boxes. If you - specify `distopia` as the backend and your simulation box is triclinic, - the function will fall back to the default `serial` backend. .. Note:: Due to the use of Instruction Set Architecture (`ISA`_) specific SIMD - intrinsics in distopia via `VCL2`_, the precision of your results may + intrinsics in distopia via `HWY`_, the precision of your results may depend on the ISA available on your machine. However, in all tested cases distopia satisfied the accuracy thresholds used to the functions in this module. Please document any issues you encounter with distopia's accuracy @@ -92,7 +97,7 @@ .. _distopia: https://github.com/MDAnalysis/distopia .. _distopia documentation: https://www.mdanalysis.org/distopia .. _ISA: https://en.wikipedia.org/wiki/Instruction_set_architecture -.. _VCL2: https://github.com/vectorclass/version2 +.. _HWY: https://github.com/google/highway .. _relevant distopia issue: https://github.com/MDAnalysis/mdanalysis/issues/3915 .. versionadded:: 0.13.0 @@ -101,6 +106,8 @@ :class:`~MDAnalysis.core.groups.AtomGroup` or an :class:`np.ndarray` .. versionchanged:: 2.5.0 Interface to the `distopia`_ package added. +.. versionchanged:: 2.9.0 + Distopia support greatly expanded (with distopia ≥ 0.4.0). Functions --------- @@ -300,7 +307,7 @@ def distance_array( ``numpy.float64``. Avoids creating the array which saves time when the function is called repeatedly. - backend : {'serial', 'OpenMP'}, optional + backend : {'serial', 'OpenMP', 'distopia'}, optional Keyword selecting the type of acceleration. Returns @@ -318,6 +325,8 @@ def distance_array( .. versionchanged:: 2.3.0 Can now accept an :class:`~MDAnalysis.core.groups.AtomGroup` as an argument in any position and checks inputs using type hinting. + .. versionchanged:: 2.9.0 + Added support for the `distopia` backend. """ confnum = configuration.shape[0] refnum = reference.shape[0] @@ -332,6 +341,14 @@ def distance_array( distances = _check_result_array(result, (refnum, confnum)) if len(distances) == 0: return distances + + if backend == "distopia": + # distopia requires that all the input arrays are the same type, + # while MDAnalysis allows for mixed types, this should be changed + # pre 3.0.0 release see issue #3707 + distances = distances.astype(np.float32) + box = np.asarray(box).astype(np.float32) if box is not None else None + if box is not None: boxtype, box = check_box(box) if boxtype == "ortho": @@ -353,6 +370,13 @@ def distance_array( backend=backend, ) + if backend == "distopia": + # mda expects the result to be in float64, so we need to convert it back + # to float64, change for 3.0, see #3707 + distances = distances.astype(np.float64) + if result is not None: + result[:] = distances + return distances @@ -388,7 +412,7 @@ def self_distance_array( Preallocated result array which must have the shape ``(n*(n-1)/2,)`` and dtype ``numpy.float64``. Avoids creating the array which saves time when the function is called repeatedly. - backend : {'serial', 'OpenMP'}, optional + backend : {'serial', 'OpenMP', 'distopia'}, optional Keyword selecting the type of acceleration. Returns @@ -411,6 +435,8 @@ def self_distance_array( .. versionchanged:: 2.3.0 Can now accept an :class:`~MDAnalysis.core.groups.AtomGroup` as an argument in any position and checks inputs using type hinting. + .. versionchanged:: 2.9.0 + Added support for the `distopia` backend. """ refnum = reference.shape[0] distnum = refnum * (refnum - 1) // 2 @@ -424,6 +450,14 @@ def self_distance_array( distances = _check_result_array(result, (distnum,)) if len(distances) == 0: return distances + + if backend == "distopia": + # distopia requires that all the input arrays are the same type, + # while MDAnalysis allows for mixed types, this should be changed + # pre 3.0.0 release see issue #3707 + distances = distances.astype(np.float32) + box = np.asarray(box).astype(np.float32) if box is not None else None + if box is not None: boxtype, box = check_box(box) if boxtype == "ortho": @@ -445,6 +479,13 @@ def self_distance_array( backend=backend, ) + if backend == "distopia": + # mda expects the result to be in float64, so we need to convert it back + # to float64, change for 3.0, see #3707 + distances = distances.astype(np.float64) + if result is not None: + result[:] = distances + return distances @@ -1630,13 +1671,17 @@ def calc_bonds( """ numatom = coords1.shape[0] bondlengths = _check_result_array(result, (numatom,)) + if backend == "distopia": + # distopia requires that all the input arrays are the same type, + # while MDAnalysis allows for mixed types, this should be changed + # pre 3.0.0 release see issue #3707 + bondlengths = bondlengths.astype(np.float32) + box = np.asarray(box).astype(np.float32) if box is not None else None if numatom > 0: if box is not None: boxtype, box = check_box(box) if boxtype == "ortho": - if backend == "distopia": - bondlengths = bondlengths.astype(np.float32) _run( "calc_bond_distance_ortho", args=(coords1, coords2, box, bondlengths), @@ -1649,15 +1694,18 @@ def calc_bonds( backend=backend, ) else: - if backend == "distopia": - bondlengths = bondlengths.astype(np.float32) _run( "calc_bond_distance", args=(coords1, coords2, bondlengths), backend=backend, ) if backend == "distopia": + # mda expects the result to be in float64, so we need to convert it back + # to float64, change for 3.0, see #3707 bondlengths = bondlengths.astype(np.float64) + if result is not None: + result[:] = bondlengths + return bondlengths @@ -1719,7 +1767,7 @@ def calc_angles( Preallocated result array of dtype ``numpy.float64`` and shape ``(n,)`` (for ``n`` coordinate triplets). Avoids recreating the array in repeated function calls. - backend : {'serial', 'OpenMP'}, optional + backend : {'serial', 'OpenMP', 'distopia'}, optional Keyword selecting the type of acceleration. Returns @@ -1746,6 +1794,13 @@ def calc_angles( numatom = coords1.shape[0] angles = _check_result_array(result, (numatom,)) + if backend == "distopia": + # distopia requires that all the input arrays are the same type, + # while MDAnalysis allows for mixed types, this should be changed + # pre 3.0.0 release see issue #3707 + angles = angles.astype(np.float32) + box = np.asarray(box).astype(np.float32) if box is not None else None + if numatom > 0: if box is not None: boxtype, box = check_box(box) @@ -1768,6 +1823,12 @@ def calc_angles( backend=backend, ) + if backend == "distopia": + # mda expects the result to be in float64, so we need to convert it back + # to float64, change for 3.0, see #3707 + angles = angles.astype(np.float64) + if result is not None: + result[:] = angles return angles @@ -1841,7 +1902,7 @@ def calc_dihedrals( Preallocated result array of dtype ``numpy.float64`` and shape ``(n,)`` (for ``n`` coordinate quadruplets). Avoids recreating the array in repeated function calls. - backend : {'serial', 'OpenMP'}, optional + backend : {'serial', 'OpenMP', 'distopia'}, optional Keyword selecting the type of acceleration. Returns @@ -1872,6 +1933,13 @@ def calc_dihedrals( numatom = coords1.shape[0] dihedrals = _check_result_array(result, (numatom,)) + if backend == "distopia": + # distopia requires that all the input arrays are the same type, + # while MDAnalysis allows for mixed types, this should be changed + # pre 3.0.0 release see issue #3707 + dihedrals = dihedrals.astype(np.float32) + box = np.asarray(box).astype(np.float32) if box is not None else None + if numatom > 0: if box is not None: boxtype, box = check_box(box) @@ -1893,7 +1961,12 @@ def calc_dihedrals( args=(coords1, coords2, coords3, coords4, dihedrals), backend=backend, ) - + if backend == "distopia": + # mda expects the result to be in float64, so we need to convert it back + # to float64, change for 3.0, see #3707 + dihedrals = dihedrals.astype(np.float64) + if result is not None: + result[:] = dihedrals return dihedrals diff --git a/testsuite/MDAnalysisTests/lib/test_distances.py b/testsuite/MDAnalysisTests/lib/test_distances.py index 8844ef9b84..67d6728a13 100644 --- a/testsuite/MDAnalysisTests/lib/test_distances.py +++ b/testsuite/MDAnalysisTests/lib/test_distances.py @@ -20,11 +20,14 @@ # MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations. # J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787 # -import itertools import sys +from unittest.mock import patch +import pytest +import numpy as np +from numpy.testing import assert_equal, assert_almost_equal, assert_allclose +import itertools from itertools import combinations_with_replacement as comb -from unittest.mock import Mock, patch - +from types import ModuleType import MDAnalysis import numpy as np import pytest @@ -34,6 +37,15 @@ from numpy.testing import assert_allclose, assert_almost_equal, assert_equal +def distopia_conditional_backend(): + # functions that allow distopia acceleration need to be tested with + # distopia backend argument but distopia is an optional dep. + if HAS_DISTOPIA: + return ["serial", "openmp", "distopia"] + else: + return ["serial", "openmp"] + + class TestCheckResultArray(object): ref = np.zeros(1, dtype=np.float64) @@ -367,7 +379,7 @@ def ref_system_universe(ref_system): ) -@pytest.mark.parametrize("backend", ["serial", "openmp"]) +@pytest.mark.parametrize("backend", distopia_conditional_backend()) class TestDistanceArray(object): @staticmethod def _dist(x, ref): @@ -515,7 +527,7 @@ def Triclinic_Universe(): return universe -@pytest.mark.parametrize("backend", ["serial", "openmp"]) +@pytest.mark.parametrize("backend", distopia_conditional_backend()) class TestDistanceArrayDCD_TRIC(object): # reasonable precision so that tests succeed on 32 and 64 bit machines # (the reference values were obtained on 64 bit) @@ -688,7 +700,7 @@ def test_atomgroup_matches_numpy_tric( ) -@pytest.mark.parametrize("backend", ["serial", "openmp"]) +@pytest.mark.parametrize("backend", distopia_conditional_backend()) class TestSelfDistanceArrayDCD_TRIC(object): prec = 5 @@ -835,7 +847,6 @@ def test_atomgroup_matches_numpy_tric( ) -@pytest.mark.parametrize("backend", ["serial", "openmp"]) class TestTriclinicDistances(object): """Unit tests for the Triclinic PBC functions. Tests: @@ -879,6 +890,7 @@ def S_mol_single(TRIC): S_mol2 = TRIC.atoms[390].position return S_mol1, S_mol2 + @pytest.mark.parametrize("backend", ["serial", "openmp"]) @pytest.mark.parametrize("S_mol", [S_mol, S_mol_single], indirect=True) def test_transforms(self, S_mol, tri_vec_box, box, backend): # To check the cython coordinate transform, the same operation is done in numpy @@ -920,9 +932,12 @@ def test_transforms(self, S_mol, tri_vec_box, box, backend): err_msg="Round trip 2 failed in transform", ) + @pytest.mark.parametrize("backend", distopia_conditional_backend()) def test_selfdist(self, S_mol, box, tri_vec_box, backend): S_mol1, S_mol2 = S_mol - R_coords = distances.transform_StoR(S_mol1, box, backend=backend) + # need to compare to serial here as transform_StoR does not have + # distopia backend support + R_coords = distances.transform_StoR(S_mol1, box, backend="serial") # Transform functions are tested elsewhere so taken as working here dists = distances.self_distance_array( R_coords, box=box, backend=backend @@ -948,7 +963,9 @@ def test_selfdist(self, S_mol, box, tri_vec_box, backend): ) # Do it again for input 2 (has wider separation in points) - R_coords = distances.transform_StoR(S_mol2, box, backend=backend) + # need to compare to serial here as transform_StoR does not have + # distopia backend support + R_coords = distances.transform_StoR(S_mol2, box, backend="serial") # Transform functions are tested elsewhere so taken as working here dists = distances.self_distance_array( R_coords, box=box, backend=backend @@ -973,11 +990,13 @@ def test_selfdist(self, S_mol, box, tri_vec_box, backend): err_msg="self_distance_array failed with input 2", ) + @pytest.mark.parametrize("backend", distopia_conditional_backend()) def test_distarray(self, S_mol, tri_vec_box, box, backend): S_mol1, S_mol2 = S_mol - - R_mol1 = distances.transform_StoR(S_mol1, box, backend=backend) - R_mol2 = distances.transform_StoR(S_mol2, box, backend=backend) + # need to compare to serial here as transform_StoR does not have + # distopia backend support + R_mol1 = distances.transform_StoR(S_mol1, box, backend="serial") + R_mol2 = distances.transform_StoR(S_mol2, box, backend="serial") # Try with box dists = distances.distance_array( @@ -998,6 +1017,7 @@ def test_distarray(self, S_mol, tri_vec_box, box, backend): dists, manual, self.prec, err_msg="distance_array failed with box" ) + @pytest.mark.parametrize("backend", distopia_conditional_backend()) def test_pbc_dist(self, S_mol, box, backend): S_mol1, S_mol2 = S_mol results = np.array([[37.629944]]) @@ -1012,6 +1032,7 @@ def test_pbc_dist(self, S_mol, box, backend): err_msg="distance_array failed to retrieve PBC distance", ) + @pytest.mark.parametrize("backend", distopia_conditional_backend()) def test_pbc_wrong_wassenaar_distance(self, backend): box = [2, 2, 2, 60, 60, 60] tri_vec_box = mdamath.triclinic_vectors(box) @@ -1068,50 +1089,31 @@ def convert_position_dtype_if_ndarray(a, b, c, d, dtype): ) -def distopia_conditional_backend(): - # functions that allow distopia acceleration need to be tested with - # distopia backend argument but distopia is an optional dep. - if HAS_DISTOPIA: - return ["serial", "openmp", "distopia"] - else: - return ["serial", "openmp"] +def test_HAS_DISTOPIA_distopia_too_old(): + # mock a version of distopia that is too old + sys.modules.pop("distopia", None) + sys.modules.pop("MDAnalysis.lib._distopia", None) + + module_name = "distopia" + mocked_module = ModuleType(module_name) + # too old version + mocked_module.__version__ = "0.1.0" + sys.modules[module_name] = mocked_module + + import MDAnalysis.lib._distopia + assert not MDAnalysis.lib._distopia.HAS_DISTOPIA -def test_HAS_DISTOPIA_incompatible_distopia(): - # warn if distopia is the wrong version and set HAS_DISTOPIA to False sys.modules.pop("distopia", None) sys.modules.pop("MDAnalysis.lib._distopia", None) - # fail any Attribute access for calc_bonds_ortho_float, - # calc_bonds_no_box_float but pretend to have the distopia - # 0.3.0 functions (from - # https://github.com/MDAnalysis/distopia/blob/main/distopia/__init__.py - # __all__): - mock_distopia_030 = Mock( - spec=[ - "calc_bonds_ortho", - "calc_bonds_no_box", - "calc_bonds_triclinic", - "calc_angles_no_box", - "calc_angles_ortho", - "calc_angles_triclinic", - "calc_dihedrals_no_box", - "calc_dihedrals_ortho", - "calc_dihedrals_triclinic", - "calc_distance_array_no_box", - "calc_distance_array_ortho", - "calc_distance_array_triclinic", - "calc_self_distance_array_no_box", - "calc_self_distance_array_ortho", - "calc_self_distance_array_triclinic", - ] - ) - with patch.dict("sys.modules", {"distopia": mock_distopia_030}): - with pytest.warns( - RuntimeWarning, match="Install 'distopia>=0.2.0,<0.3.0' to" - ): - import MDAnalysis.lib._distopia - assert not MDAnalysis.lib._distopia.HAS_DISTOPIA + # new enough version + mocked_module.__version__ = "0.4.0" + sys.modules[module_name] = mocked_module + + import MDAnalysis.lib._distopia + + assert MDAnalysis.lib._distopia.HAS_DISTOPIA class TestCythonFunctions(object): @@ -1248,6 +1250,25 @@ def test_bonds(self, box, backend, dtype, pos, request): err_msg="PBC check #w with box", ) + @pytest.mark.parametrize("dtype", (np.float32, np.float64)) + @pytest.mark.parametrize("backend", distopia_conditional_backend()) + def test_calc_bonds_results_inplace_all_backends( + self, + backend, + dtype, + ): + N = 10 + c0 = np.ones(3 * N, dtype=dtype).reshape(N, 3) * 2 + c1 = np.ones(3 * N, dtype=dtype).reshape(N, 3) * 3 + + result = np.zeros(N, dtype=np.float64) + distances.calc_bonds(c0, c1, result=result, backend=backend) + expected = np.ones(N, dtype=dtype) * 3 ** (1 / 2) + # test the result array is updated in place + assert_almost_equal( + result, expected, self.prec, err_msg="calc_bonds inplace failed" + ) + @pytest.mark.parametrize("backend", distopia_conditional_backend()) def test_bonds_badbox(self, positions, backend): a, b, c, d = positions @@ -1314,7 +1335,7 @@ def test_bonds_single_coords(self, shift, periodic, backend): @pytest.mark.parametrize("dtype", (np.float32, np.float64)) @pytest.mark.parametrize("pos", ["positions", "positions_atomgroups"]) - @pytest.mark.parametrize("backend", ["serial", "openmp"]) + @pytest.mark.parametrize("backend", distopia_conditional_backend()) def test_angles(self, backend, dtype, pos, request): a, b, c, d = request.getfixturevalue(pos) a, b, c, d = convert_position_dtype_if_ndarray(a, b, c, d, dtype) @@ -1344,7 +1365,27 @@ def test_angles(self, backend, dtype, pos, request): err_msg="Small angle failed in calc_angles", ) - @pytest.mark.parametrize("backend", ["serial", "openmp"]) + @pytest.mark.parametrize("dtype", (np.float32, np.float64)) + @pytest.mark.parametrize("backend", distopia_conditional_backend()) + def test_calc_angles_results_inplace_all_backends( + self, + backend, + dtype, + ): + N = 10 + c0 = np.ones(3 * N, dtype=dtype).reshape(N, 3) * 2 + c1 = np.ones(3 * N, dtype=dtype).reshape(N, 3) * 3 + c2 = np.ones(3 * N, dtype=dtype).reshape(N, 3) * 4 + + result = np.zeros(N, dtype=np.float64) + distances.calc_angles(c0, c1, c2, result=result, backend=backend) + expected = np.ones(N, dtype=dtype) * np.pi + # test the result array is updated in place + assert_almost_equal( + result, expected, self.prec, err_msg="calc_angles inplace failed" + ) + + @pytest.mark.parametrize("backend", distopia_conditional_backend()) def test_angles_bad_result(self, positions, backend): a, b, c, d = positions badresult = np.zeros(len(a) - 1) # Bad result array @@ -1370,7 +1411,7 @@ def test_angles_bad_result(self, positions, backend): ) @pytest.mark.parametrize("shift", shifts) @pytest.mark.parametrize("periodic", [True, False]) - @pytest.mark.parametrize("backend", ["serial", "openmp"]) + @pytest.mark.parametrize("backend", distopia_conditional_backend()) def test_angles_single_coords(self, case, shift, periodic, backend): def manual_angle(x, y, z): return mdamath.angle(y - x, y - z) @@ -1391,7 +1432,7 @@ def manual_angle(x, y, z): @pytest.mark.parametrize("dtype", (np.float32, np.float64)) @pytest.mark.parametrize("pos", ["positions", "positions_atomgroups"]) - @pytest.mark.parametrize("backend", ["serial", "openmp"]) + @pytest.mark.parametrize("backend", distopia_conditional_backend()) def test_dihedrals(self, backend, dtype, pos, request): a, b, c, d = request.getfixturevalue(pos) a, b, c, d = convert_position_dtype_if_ndarray(a, b, c, d, dtype) @@ -1404,8 +1445,9 @@ def test_dihedrals(self, backend, dtype, pos, request): ) assert np.isnan(dihedrals[0]), "Zero length dihedral failed" assert np.isnan(dihedrals[1]), "Straight line dihedral failed" + # 180 degree dihedral can be either pi or (-pi for distopia) assert_almost_equal( - dihedrals[2], + np.abs(dihedrals[2]), np.pi, self.prec, err_msg="180 degree dihedral failed", @@ -1417,7 +1459,39 @@ def test_dihedrals(self, backend, dtype, pos, request): err_msg="arbitrary dihedral angle failed", ) - @pytest.mark.parametrize("backend", ["serial", "openmp"]) + @pytest.mark.parametrize("dtype", (np.float32, np.float64)) + @pytest.mark.parametrize("backend", distopia_conditional_backend()) + def test_calc_dihedrals_results_inplace_all_backends( + self, + backend, + dtype, + ): + N = 10 + c0 = np.ones(3 * N, dtype=dtype).reshape(N, 3) * 0 # 0,0,0 + c1 = np.ones(3 * N, dtype=dtype).reshape(N, 3) * 1 # 1,1,1 + # now make a 90 degree angle, with 2,2,1 + c2 = np.ones(3 * N, dtype=dtype).reshape(N, 3) + c2[:, 0] = 2 + c2[:, 1] = 2 + # now back to 0 on z axis + c3 = np.ones(3 * N, dtype=dtype).reshape(N, 3) * 0 + c3[:, 0] = 3 + c3[:, 1] = 3 + + result = np.zeros(N, dtype=np.float64) + distances.calc_dihedrals( + c0, c1, c2, c3, result=result, backend=backend + ) + expected = np.ones(N, dtype=dtype) * 0 + # test the result array is updated in place + assert_almost_equal( + result, + expected, + self.prec, + err_msg="calc_dihedrals inplace failed", + ) + + @pytest.mark.parametrize("backend", distopia_conditional_backend()) def test_dihedrals_wronglength(self, positions, wronglength, backend): a, b, c, d = positions with pytest.raises(ValueError): @@ -1432,7 +1506,7 @@ def test_dihedrals_wronglength(self, positions, wronglength, backend): with pytest.raises(ValueError): distances.calc_dihedrals(a, b, c, wronglength, backend=backend) - @pytest.mark.parametrize("backend", ["serial", "openmp"]) + @pytest.mark.parametrize("backend", distopia_conditional_backend()) def test_dihedrals_bad_result(self, positions, backend): a, b, c, d = positions badresult = np.zeros(len(a) - 1) # Bad result array @@ -1491,7 +1565,7 @@ def test_dihedrals_bad_result(self, positions, backend): ) @pytest.mark.parametrize("shift", shifts) @pytest.mark.parametrize("periodic", [True, False]) - @pytest.mark.parametrize("backend", ["serial", "openmp"]) + @pytest.mark.parametrize("backend", distopia_conditional_backend()) def test_dihedrals_single_coords(self, case, shift, periodic, backend): def manual_dihedral(a, b, c, d): return mdamath.dihedral(b - a, c - b, d - c) @@ -1526,7 +1600,7 @@ def test_numpy_compliance_bonds(self, positions, backend): err_msg="Cython bonds didn't match numpy calculations", ) - @pytest.mark.parametrize("backend", ["serial", "openmp"]) + @pytest.mark.parametrize("backend", distopia_conditional_backend()) def test_numpy_compliance_angles(self, positions, backend): a, b, c, d = positions # Checks that the cython functions give identical results to the numpy versions @@ -1544,7 +1618,7 @@ def test_numpy_compliance_angles(self, positions, backend): err_msg="Cython angles didn't match numpy calcuations", ) - @pytest.mark.parametrize("backend", ["serial", "openmp"]) + @pytest.mark.parametrize("backend", distopia_conditional_backend()) def test_numpy_compliance_dihedrals(self, positions, backend): a, b, c, d = positions # Checks that the cython functions give identical results to the numpy versions @@ -1555,6 +1629,8 @@ def test_numpy_compliance_dihedrals(self, positions, backend): dihedrals_numpy = np.array( [mdamath.dihedral(x, y, z) for x, y, z in zip(ab, bc, cd)] ) + # 180 (and 360) degree dihedral can be either pi or (-pi for distopia) + dihedrals[2] = np.abs(dihedrals[2]) assert_almost_equal( dihedrals, dihedrals_numpy, @@ -1920,7 +1996,7 @@ def test_input_unchanged_calc_bonds_atomgroup( assert_equal([crd.positions for crd in crds], refs) @pytest.mark.parametrize("box", boxes) - @pytest.mark.parametrize("backend", ["serial", "openmp"]) + @pytest.mark.parametrize("backend", distopia_conditional_backend()) def test_input_unchanged_calc_angles(self, coords, box, backend): crds = coords[:3] refs = [crd.copy() for crd in crds] @@ -1930,7 +2006,7 @@ def test_input_unchanged_calc_angles(self, coords, box, backend): assert_equal(crds, refs) @pytest.mark.parametrize("box", boxes) - @pytest.mark.parametrize("backend", ["serial", "openmp"]) + @pytest.mark.parametrize("backend", distopia_conditional_backend()) def test_input_unchanged_calc_angles_atomgroup( self, coords_atomgroups, box, backend ): @@ -1942,7 +2018,7 @@ def test_input_unchanged_calc_angles_atomgroup( assert_equal([crd.positions for crd in crds], refs) @pytest.mark.parametrize("box", boxes) - @pytest.mark.parametrize("backend", ["serial", "openmp"]) + @pytest.mark.parametrize("backend", distopia_conditional_backend()) def test_input_unchanged_calc_dihedrals(self, coords, box, backend): crds = coords refs = [crd.copy() for crd in crds] @@ -1952,7 +2028,7 @@ def test_input_unchanged_calc_dihedrals(self, coords, box, backend): assert_equal(crds, refs) @pytest.mark.parametrize("box", boxes) - @pytest.mark.parametrize("backend", ["serial", "openmp"]) + @pytest.mark.parametrize("backend", distopia_conditional_backend()) def test_input_unchanged_calc_dihedrals_atomgroup( self, coords_atomgroups, box, backend ): @@ -2013,7 +2089,7 @@ def empty_coord(): return np.empty((0, 3), dtype=np.float32) @pytest.mark.parametrize("box", boxes) - @pytest.mark.parametrize("backend", ["serial", "openmp"]) + @pytest.mark.parametrize("backend", distopia_conditional_backend()) def test_empty_input_distance_array(self, empty_coord, box, backend): res = distances.distance_array( empty_coord, empty_coord, box=box, backend=backend @@ -2021,7 +2097,7 @@ def test_empty_input_distance_array(self, empty_coord, box, backend): assert_equal(res, np.empty((0, 0), dtype=np.float64)) @pytest.mark.parametrize("box", boxes) - @pytest.mark.parametrize("backend", ["serial", "openmp"]) + @pytest.mark.parametrize("backend", distopia_conditional_backend()) def test_empty_input_self_distance_array(self, empty_coord, box, backend): res = distances.self_distance_array( empty_coord, box=box, backend=backend @@ -2092,7 +2168,7 @@ def test_empty_input_calc_bonds(self, empty_coord, box, backend): assert_equal(res, np.empty((0,), dtype=np.float64)) @pytest.mark.parametrize("box", boxes) - @pytest.mark.parametrize("backend", ["serial", "openmp"]) + @pytest.mark.parametrize("backend", distopia_conditional_backend()) def test_empty_input_calc_angles(self, empty_coord, box, backend): res = distances.calc_angles( empty_coord, empty_coord, empty_coord, box=box, backend=backend @@ -2100,7 +2176,7 @@ def test_empty_input_calc_angles(self, empty_coord, box, backend): assert_equal(res, np.empty((0,), dtype=np.float64)) @pytest.mark.parametrize("box", boxes) - @pytest.mark.parametrize("backend", ["serial", "openmp"]) + @pytest.mark.parametrize("backend", distopia_conditional_backend()) def test_empty_input_calc_dihedrals(self, empty_coord, box, backend): res = distances.calc_dihedrals( empty_coord, @@ -2283,7 +2359,7 @@ def test_output_type_calc_bonds(self, incoords, box, backend): @pytest.mark.parametrize( "incoords", [3 * [coords[0]]] + list(comb(coords[1:], 3)) ) - @pytest.mark.parametrize("backend", ["serial", "openmp"]) + @pytest.mark.parametrize("backend", distopia_conditional_backend()) def test_output_type_calc_angles(self, incoords, box, backend): res = distances.calc_angles(*incoords, box=box, backend=backend) maxdim = max([crd.ndim for crd in incoords]) @@ -2299,7 +2375,7 @@ def test_output_type_calc_angles(self, incoords, box, backend): @pytest.mark.parametrize( "incoords", [4 * [coords[0]]] + list(comb(coords[1:], 4)) ) - @pytest.mark.parametrize("backend", ["serial", "openmp"]) + @pytest.mark.parametrize("backend", distopia_conditional_backend()) def test_output_type_calc_dihedrals(self, incoords, box, backend): res = distances.calc_dihedrals(*incoords, box=box, backend=backend) maxdim = max([crd.ndim for crd in incoords])