Survival probability: Intermittency, Step, Residues, Overlapping SP Selection Example

package/MDAnalysis/analysis/waterdynamics.py

@@ -1273,15 +1273,42 @@ def run(self, tau_max=20, start=0, stop=None, step=1, verbose=False):
         if tau_max > (stop - start):
             raise ValueError("Too few frames selected for given tau_max.")
 
-        # load all frames to an array of sets
+        # preload the frames (selected ids) to a list of sets
         selected_ids = []
-        for ts in self.universe.trajectory[start:stop]:
-            self.print(verbose, "Loading frame:", ts)
-            selected_ids.append(set(self.universe.select_atoms(self.selection).ids))
+        if step <= tau_max + 1:
+            # load all frames
+            for ts in self.universe.trajectory[start:stop]:
+                self.print(verbose, "Loading frame:", ts)
+                selected_ids.append(set(self.universe.select_atoms(self.selection).ids))
+        else:
+            # skip frames that will not be used
+            # Example: step 5 and tau 2: L, L, L, S, S, L, L, L, S, S, ... with L = Load, and S = Skip
+            loaded_counter = 0
+            no_of_frames_to_load = tau_max + 1
+            skipped = 0
+            no_of_frames_to_skip = step - (tau_max + 1)
+            for ts in self.universe.trajectory[start:stop]:
+                if skipped == no_of_frames_to_skip:
+                    skipped = 0
+                    loaded_counter = 0
+
+                if loaded_counter == no_of_frames_to_load:
+                    self.print(verbose, "Skipping loading frame:", ts)
+                    skipped += 1
+                    continue
+
+                self.print(verbose, "Loading frame:", ts)
+                selected_ids.append(set(self.universe.select_atoms(self.selection).ids))
+
+                loaded_counter += 1
+
 
         tau_timeseries = np.arange(1, tau_max + 1)
         sp_timeseries_data = [[] for _ in range(tau_max)]
 
+        # frames not analysed are skipped because they were not loaded
+        step = tau_max + 1 if step >= (tau_max + 1) else step
+
         for t in range(0, len(selected_ids), step):
             Nt = len(selected_ids[t])
 

testsuite/MDAnalysisTests/analysis/test_waterdynamics.py

@@ -31,7 +31,7 @@
 import numpy as np
 from mock import patch
 from mock import Mock
-from numpy.testing import assert_almost_equal
+from numpy.testing import assert_almost_equal, assert_equal
 
 SELECTION1 = "byres name OH2"
 SELECTION2 = "byres name P1"
@@ -120,3 +120,15 @@ def test_SurvivalProbability_alwaysPresent(universe):
         sp = waterdynamics.SurvivalProbability(universe, "")
         sp.run(tau_max=3, start=0, stop=6)
         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)
+        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'
+        assert universe.trajectory.n_frames > 6
+        assert_equal(select_atoms_mock.call_count, 6)