diff --git a/package/CHANGELOG b/package/CHANGELOG index e3c5e01ff44..d33c2b4c92d 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -135,6 +135,11 @@ Enhancements the capability to allow intermittent behaviour (PR #2256) Changes + * Removes the deprecated `t0`, `tf`, and `dtmax` from + :class:Waterdynamics.SurvivalProbability. Instead the `start`, `stop` and + `tau_max` keywords should be passed to + :meth:`Waterdynamics.SurvivalProbability.run`. Furthermore, the `stop` + keyword is now exclusive instead of inclusive (Issue #2524). * :meth:`align.fasta2select` now writes `alnfilename` and `treefilename` to the current working directory (Issue #2701) * Removes duplicate `NUMBER_TO_ELEMENTS` table from topology._elements, diff --git a/package/MDAnalysis/analysis/waterdynamics.py b/package/MDAnalysis/analysis/waterdynamics.py index 70d42a26a93..b7f0e206677 100644 --- a/package/MDAnalysis/analysis/waterdynamics.py +++ b/package/MDAnalysis/analysis/waterdynamics.py @@ -289,7 +289,7 @@ universe = MDAnalysis.Universe(pdb, trajectory) select = "byres name OH2 and sphzone 12.3 (resid 42 or resid 26) " sp = SP(universe, select, verbose=True) - sp.run(start=0, stop=100, tau_max=20) + sp.run(start=0, stop=101, tau_max=20) tau_timeseries = sp.tau_timeseries sp_timeseries = sp.sp_timeseries @@ -304,6 +304,12 @@ plt.plot(tau_timeseries, sp_timeseries) plt.show() +One should note that the `stop` keyword as used in the above example has an +`exclusive` behaviour, i.e. here the final frame used will be 100 not 101. +This behaviour is aligned with :class:`AnalysisBase` but currently differs from +other :mod:`MDAnalysis.analysis.waterdynamics` classes, which all exhibit +`inclusive` behaviour for their final frame selections. + Another example applies to the situation where you work with many different "residues". Here we calculate the SP of a potassium ion around each lipid in a membrane and average the results. In this example, if the SP analysis were run without treating each lipid @@ -1171,85 +1177,94 @@ class SurvivalProbability(object): When True, prints progress and comments to the console. + Notes + ----- + Currently :class:`SurvivalProbability` is the only on in + :mod:`MDAnalysis.analysis.waterdynamics` to support an `exclusive` + behaviour (i.e. similar to the current behaviour of :class:`AnalysisBase` + to the `stop` keyword passed to :meth:`SurvivalProbability.run`. Unlike + other :mod:`MDAnalysis.analysis.waterdynamics` final frame definitions + which are `inclusive`. + + .. versionadded:: 0.11.0 - .. versionchanged:: 0.21.0 - Using the MDAnalysis.lib.correlations.py to carry out the intermittency and autocorrelation calculations .. versionchanged:: 1.0.0 - Changed `selection` keyword to `select` + Using the MDAnalysis.lib.correlations.py to carry out the intermittency + and autocorrelation calculations. + Changed `selection` keyword to `select`. + Removed support for the deprecated `t0`, `tf`, and `dtmax` keywords. + These should instead be passed to :meth:`SurvivalProbability.run` as + the `start`, `stop`, and `tau_max` keywords respectively. + The `stop` keyword as passed to :meth:`SurvivalProbability.run` has now + changed behaviour and will act in an `exclusive` manner (instead of it's + previous `inclusive` behaviour), """ - def __init__(self, universe, select, t0=None, tf=None, dtmax=None, verbose=False): + def __init__(self, universe, select, verbose=False): self.universe = universe self.selection = select self.verbose = verbose - # backward compatibility - self.start = self.stop = self.tau_max = None - if t0 is not None: - self.start = t0 - warnings.warn("t0 is deprecated, use run(start=t0) instead", category=DeprecationWarning) - - if tf is not None: - self.stop = tf - warnings.warn("tf is deprecated, use run(stop=tf) instead", category=DeprecationWarning) - - if dtmax is not None: - self.tau_max = dtmax - warnings.warn("dtmax is deprecated, use run(tau_max=dtmax) instead", category=DeprecationWarning) - - - def run(self, tau_max=20, start=0, stop=None, step=1, residues=False, intermittency=0, verbose=False): + def run(self, tau_max=20, start=None, stop=None, step=None, residues=False, + intermittency=0, verbose=False): """ Computes and returns the Survival Probability (SP) timeseries Parameters ---------- start : int, optional - Zero-based index of the first frame to be analysed + Zero-based index of the first frame to be analysed, Default: None + (first frame). stop : int, optional - Zero-based index of the last frame to be analysed (inclusive) + Zero-based index of the last frame to be analysed (exclusive), + Default: None (last frame). step : int, optional - Jump every `step`-th frame. This is compatible but independant of the taus used, and it is good to consider - using the `step` equal to `tau_max` to remove the overlap. - Note that `step` and `tau_max` work consistently with intermittency. + Jump every `step`-th frame. This is compatible but independant of + the taus used, and it is good to consider using the `step` equal + to `tau_max` to remove the overlap. Note that `step` and `tau_max` + work consistently with intermittency. Default: None + (use every frame). tau_max : int, optional - Survival probability is calculated for the range 1 <= `tau` <= `tau_max` + Survival probability is calculated for the range + 1 <= `tau` <= `tau_max`. residues : Boolean, optional - If true, the analysis will be carried out on the residues (.resids) rather than on atom (.ids). - A single atom is sufficient to classify the residue as within the distance. + If true, the analysis will be carried out on the residues + (.resids) rather than on atom (.ids). A single atom is sufficient + to classify the residue as within the distance. intermittency : int, optional - The maximum number of consecutive frames for which an atom can leave but be counted as present if it returns - at the next frame. An intermittency of `0` is equivalent to a continuous survival probability, which does - not allow for the leaving and returning of atoms. For example, for `intermittency=2`, any given atom may - leave a region of interest for up to two consecutive frames yet be treated as being present at all frames. - The default is continuous (0). + The maximum number of consecutive frames for which an atom can + leave but be counted as present if it returns at the next frame. + An intermittency of `0` is equivalent to a continuous survival + probability, which does not allow for the leaving and returning of + atoms. For example, for `intermittency=2`, any given atom may leave + a region of interest for up to two consecutive frames yet be + treated as being present at all frames. The default is continuous + (0). verbose : Boolean, optional - Print the progress to the console + Print the progress to the console. Returns ------- tau_timeseries : list tau from 1 to `tau_max`. Saved in the field tau_timeseries. sp_timeseries : list - survival probability for each value of `tau`. Saved in the field sp_timeseries. + survival probability for each value of `tau`. Saved in the field + sp_timeseries. sp_timeseries_data: list raw datapoints from which the average is taken (sp_timeseries). Time dependancy and distribution can be extracted. - """ - # backward compatibility (and priority) - start = self.start if self.start is not None else start - stop = self.stop if self.stop is not None else stop - tau_max = self.tau_max if self.tau_max is not None else tau_max - # sanity checks - if stop is not None and stop >= len(self.universe.trajectory): - raise ValueError("\"stop\" must be smaller than the number of frames in the trajectory.") + .. versionchanged:: 1.0.0 + To math other analysis methods, the `stop` keyword is now exclusive + rather than inclusive. + """ - if stop is None: - stop = len(self.universe.trajectory) - else: - stop = stop + 1 + start, stop, step = self.universe.trajectory.check_slice_indices( + start, + stop, + step + ) if tau_max > (stop - start): raise ValueError("Too few frames selected for given tau_max.") diff --git a/testsuite/MDAnalysisTests/analysis/test_waterdynamics.py b/testsuite/MDAnalysisTests/analysis/test_waterdynamics.py index 3d2325e1a21..230ac63f1fd 100644 --- a/testsuite/MDAnalysisTests/analysis/test_waterdynamics.py +++ b/testsuite/MDAnalysisTests/analysis/test_waterdynamics.py @@ -95,7 +95,7 @@ def test_SurvivalProbability_intermittency1and2(universe): ids = [(9, 8), (), (8,), (9,), (8,), (), (9, 8), (), (8,), (9, 8)] select_atoms_mock.side_effect = lambda selection: Mock(ids=ids.pop()) # atom IDs fed set by set sp = waterdynamics.SurvivalProbability(universe, "") - sp.run(tau_max=3, stop=9, verbose=True, intermittency=2) + sp.run(tau_max=3, stop=10, verbose=True, intermittency=2) assert all(x == {9, 8} for x in sp._intermittent_selected_ids) assert_almost_equal(sp.sp_timeseries, [1, 1, 1, 1]) @@ -109,7 +109,7 @@ def test_SurvivalProbability_intermittency2lacking(universe): ids = [(9,), (), (), (), (9,), (), (), (), (9,)] select_atoms_mock.side_effect = lambda selection: Mock(ids=ids.pop()) # atom IDs fed set by set sp = waterdynamics.SurvivalProbability(universe, "") - sp.run(tau_max=3, stop=8, verbose=True, intermittency=2) + sp.run(tau_max=3, stop=9, verbose=True, intermittency=2) assert_almost_equal(sp.sp_timeseries, [1, 0, 0, 0]) @@ -122,7 +122,7 @@ def test_SurvivalProbability_intermittency1_step5_noSkipping(universe): ids = [(2, 3), (3,), (2, 3), (3,), (2,), (3,), (2, 3), (3,), (2, 3), (2, 3)] select_atoms_mock.side_effect = lambda selection: Mock(ids=ids.pop()) # atom IDs fed set by set sp = waterdynamics.SurvivalProbability(universe, "") - sp.run(tau_max=2, stop=9, verbose=True, intermittency=1, step=5) + sp.run(tau_max=2, stop=10, verbose=True, intermittency=1, step=5) assert all((x == {2, 3} for x in sp._intermittent_selected_ids)) assert_almost_equal(sp.sp_timeseries, [1, 1, 1]) @@ -137,7 +137,7 @@ def test_SurvivalProbability_intermittency1_step5_Skipping(universe): beforepopsing = len(ids) - 2 select_atoms_mock.side_effect = lambda selection: Mock(ids=ids.pop()) # atom IDs fed set by set sp = waterdynamics.SurvivalProbability(universe, "") - sp.run(tau_max=1, stop=9, verbose=True, intermittency=1, step=5) + sp.run(tau_max=1, stop=10, verbose=True, intermittency=1, step=5) assert all((x == {1} for x in sp._intermittent_selected_ids)) assert len(sp._selected_ids) == beforepopsing assert_almost_equal(sp.sp_timeseries, [1, 1]) @@ -235,7 +235,7 @@ def test_SurvivalProbability_t0tf(universe): ids = [(0, ), (0, ), (7, 6, 5), (6, 5, 4), (5, 4, 3), (4, 3, 2), (3, 2, 1), (0, )] select_atoms_mock.side_effect = lambda selection: Mock(ids=ids.pop(2)) # atom IDs fed set by set sp = waterdynamics.SurvivalProbability(universe, "") - sp.run(tau_max=3, start=2, stop=6) + sp.run(tau_max=3, start=2, stop=7) assert_almost_equal(sp.sp_timeseries, [1, 2 / 3.0, 1 / 3.0, 0]) @@ -244,7 +244,7 @@ def test_SurvivalProbability_definedTaus(universe): ids = [(9, 8, 7), (8, 7, 6), (7, 6, 5), (6, 5, 4), (5, 4, 3), (4, 3, 2), (3, 2, 1)] select_atoms_mock.side_effect = lambda selection: Mock(ids=ids.pop()) # atom IDs fed set by set sp = waterdynamics.SurvivalProbability(universe, "") - sp.run(tau_max=3, start=0, stop=6, verbose=True) + sp.run(tau_max=3, start=0, stop=7, verbose=True) assert_almost_equal(sp.sp_timeseries, [1, 2 / 3.0, 1 / 3.0, 0]) @@ -252,7 +252,7 @@ def test_SurvivalProbability_zeroMolecules(universe): # no atom IDs found with patch.object(universe, 'select_atoms', return_value=Mock(ids=[])) as select_atoms_mock: sp = waterdynamics.SurvivalProbability(universe, "") - sp.run(tau_max=3, start=3, stop=6, verbose=True) + sp.run(tau_max=3, start=3, stop=7, verbose=True) assert all(np.isnan(sp.sp_timeseries[1:])) @@ -260,7 +260,7 @@ def test_SurvivalProbability_alwaysPresent(universe): # always the same atom IDs found, 7 and 8 with patch.object(universe, 'select_atoms', return_value=Mock(ids=[7, 8])) as select_atoms_mock: sp = waterdynamics.SurvivalProbability(universe, "") - sp.run(tau_max=3, start=0, stop=6, verbose=True) + sp.run(tau_max=3, start=0, stop=7, verbose=True) assert all(np.equal(sp.sp_timeseries, 1)) @@ -268,7 +268,7 @@ def test_SurvivalProbability_stepLargerThanDtmax(universe): # Testing if the frames are skipped correctly with patch.object(universe, 'select_atoms', return_value=Mock(ids=(1,))) as select_atoms_mock: sp = waterdynamics.SurvivalProbability(universe, "") - sp.run(tau_max=2, step=5, stop=9, verbose=True) + sp.run(tau_max=2, step=5, stop=10, verbose=True) assert_equal(sp.sp_timeseries, [1, 1, 1]) # with tau_max=2 for all the frames we only read 6 of them # this is because the frames which are not used are skipped, and therefore 'select_atoms' @@ -279,6 +279,6 @@ def test_SurvivalProbability_stepLargerThanDtmax(universe): def test_SurvivalProbability_stepEqualDtMax(universe): with patch.object(universe, 'select_atoms', return_value=Mock(ids=(1,))) as select_atoms_mock: sp = waterdynamics.SurvivalProbability(universe, "") - sp.run(tau_max=4, step=5, stop=9, verbose=True) + sp.run(tau_max=4, step=5, stop=10, verbose=True) # all frames from 0, with 9 inclusive assert_equal(select_atoms_mock.call_count, 10)