Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Waterdynamics.SurvivalProbability start/stop/step/tau_max #2525

Merged
5 changes: 5 additions & 0 deletions package/CHANGELOG
Original file line number Diff line number Diff line change
Expand Up @@ -82,6 +82,11 @@ Enhancements
* Improve the distance search in water bridge analysis with capped_distance (PR #2480)

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).
* Removes `save()` function from contacts, diffusionmap, hole, LinearDensity,
and rms (Issue #1745).
* Removes; `save_table()` from :class:`HydrogenBondAnalysis`,
Expand Down
85 changes: 41 additions & 44 deletions package/MDAnalysis/analysis/waterdynamics.py
Original file line number Diff line number Diff line change
Expand Up @@ -1237,28 +1237,17 @@ class SurvivalProbability(object):


.. versionadded:: 0.11.0

.. versionchanged:: 1.0.0
Removes 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.
"""

def __init__(self, universe, selection, t0=None, tf=None, dtmax=None, verbose=False):
def __init__(self, universe, selection, verbose=False):
self.universe = universe
self.selection = selection
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 print(self, verbose, *args):
if self.verbose:
print(args)
Expand Down Expand Up @@ -1316,55 +1305,63 @@ def _correct_intermittency(self, intermittency, selected_ids):
seen_frames_ago[atomid] = 0


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,
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There is no longer backward compatibility warnings? why?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@alejob by this do you mean the removal of t0, tf, and dtmax? These have been slated for deprecation for a long time. As part of the upcoming 1.0.0 release, we are aiming to remove all historical deprecated code.

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`.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If stop is exclusive does this docstring need updating?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not sure here. tau_max should be independent of stop, i.e. previously the case was that stop was set in the code to be stop+1. The changes done here only change it so that stop=stop, tau and tau_max's behaviour shouldn't be affected (if that makes sense).

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Exactly, tau_max is independent of stop.

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.

# 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.")
Expand Down
22 changes: 11 additions & 11 deletions testsuite/MDAnalysisTests/analysis/test_waterdynamics.py
Original file line number Diff line number Diff line change
Expand Up @@ -91,7 +91,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, [2 / 3.0, 1 / 3.0, 0])


Expand All @@ -100,31 +100,31 @@ 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, [2 / 3.0, 1 / 3.0, 0])


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))


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))


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])
# 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'
Expand All @@ -135,7 +135,7 @@ 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)

Expand All @@ -149,7 +149,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 == set([9, 8]) for x in sp.selected_ids))
assert_almost_equal(sp.sp_timeseries, [1, 1, 1])

Expand All @@ -162,7 +162,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, [0, 0, 0])


Expand All @@ -175,7 +175,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 == set([2, 3]) for x in sp.selected_ids))
assert_almost_equal(sp.sp_timeseries, [1, 1])

Expand All @@ -189,7 +189,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 == set([1]) for x in sp.selected_ids))
assert len(sp.selected_ids) == beforepopsing
assert_almost_equal(sp.sp_timeseries, [1])
assert_almost_equal(sp.sp_timeseries, [1])