SurvivalProbability - Readability, performance, and algorithm changes

[bieniekmateusz]

Small changes:

• Performance (16 minutes to 60 seconds):
  -- issue due to the unnecessary function overhead,
  -- iterating over the dataset only once (RAM memory optimisation)
• More thorough documentation of parameters and additional sanity checks in init
• Removed temporarily the progress meter which only worked for loading the simulation dataset, confusing the reader
• Allows the user to access the distribution for each tau through sp.sp_timeseries
• Returns the taus along with the timeseries (opening access for dt). Not returning tau = 0

Validity Test Cases Changes:

• Replaced the arbitrary test cases with defined datasets and predictable taus
• Added a test case where no atoms IDs are found. In this case, we return tau = NaN. This is because 0 means that the initial molecular leaves after a given tau. However, if there is no initial molecules in the first place, this would be the wrong conclusion.

Algorithmic Changes made in this Pull Request - Request for Contribution from the Original Author:

• We modified the algorithm to be applicable to other atom groups (e.g. ions) besides water. Particularly, we changed the way the algorithm behaves in the case where there is no molecules found for the reference frame. Previously, one would expect the reference frame to always find some water molecules. If this was not the case, each tau was diluted - divided by (Nt + 1). Because it was unlikely to happen with any reasonable water selection, this case was not important. However, let us consider using SP the survival of ions in a given area, with frames where no ions are found in the first place (reference frame). For such frames, diluting the taus means that the ions leave earlier - whereas in fact, these ions might rarely diffuse into the selected area. The rare diffusion into the selected area should not affect how long they stay (ie their SurvivalProbability). Please let us know if we missed anything.

Questions:

• Since this algorithm is applicable to any mobile atoms in a system, we suggest this package is moved out of waterdynamics into a more general package.

Future considerations:

• Adding progress bar
• Adding dt to allow sample the dataset with a faster shift (and possibly discuss data correlations for larger taus)

Edit: Acknowledgement - This work was done together with @p-j-smith, but we submitted most of it on my account.

PR Checklist

• [X ] Tests?
• Docs?
• CHANGELOG updated?
• Issue raised/referenced?

[richardjgowers]

@bieniekmateusz thanks, this looks like it will make some good changes.

WRT "faster shift for dt", are you talking about doing every 2nd frame rather than every frame?

[bieniekmateusz]

@richardjgowers Thanks. Also, I am still to add your test case for the t0case which I hope to add in the next version with the progress meter.

Your question: Yes, every nth frame. We are considering letting the user specify the jump over the trajectory. However, we want to first check how, for example, tau = 20 is affected by being calculated for overlapping datasets (t=0-20, t=1-21, ...). We are not sure how (and if) that affects the SP yet. Any thoughts on it are welcome.

[richardjgowers]

For the case where no molecules are found, it makes sense to disregard this run as you can't measure anything. It's not that the molecule had 0 survival probability.

So for hydrogen bond autocorrelation I let people define a total duration (tau_max here) and a number of frames within that duration to consider. I'm not completely happy with that solution either, so it would be good to define a clear & comprehensive way/language to define this sort of stuff.

I need to read through this module fully but..

Reducing the resolution of frames should be fine in some cases. We can actually check after the calculation whether the calculated tau makes sense with respect to how we sampled it. Then issue some sort of warning if the data is probably junk (and suggest better settings).

Overlapping frames should be fine. The concentration in an area will remain constant, so by starting a new sample in the same window will identify different "starting molecules". It will be slightly correlated with the old sample though (so isn't 100% new data)

[codecov]

Codecov Report

Merging #1995 into develop will decrease coverage by 0.08%.
The diff coverage is 57.77%.

@@             Coverage Diff             @@
##           develop    #1995      +/-   ##
===========================================
- Coverage    88.93%   88.84%   -0.09%     
===========================================
  Files          144      144              
  Lines        17490    17497       +7     
  Branches      2693     2702       +9     
===========================================
- Hits         15554    15545       -9     
- Misses        1323     1332       +9     
- Partials       613      620       +7

Impacted Files	Coverage Δ
package/MDAnalysis/analysis/waterdynamics.py	80.94% <57.77%> (-3.51%)	⬇️

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update d2e22ff...8106a6c. Read the comment docs.

[bieniekmateusz]

I've updated the last test case (t0) to avoid arbitrary datasets and to test tf in the same test case.

[bieniekmateusz]

I did the rebase which was not too much trouble at the end. Cheers

[bieniekmateusz]

Error in testsuite/MDAnalysisTests/topology/test_pqr.py:93:

Could it be that the rebasing introduced this issue since I have not touched the topology files. Do you have any suggestions on what I should do to correct the tests? Thanks

[bieniekmateusz]

I've updated the example to match the data being stored in the object's fields.

[richardjgowers]

LGTM, thanks @bieniekmateusz

[bieniekmateusz]

and @p-j-smith!

We're happy to hear that. And thanks to you too!