Unexpected dt when using Universe.transfer_to_memory(step=x) (0.16.0-dev0)

[kaplajon]

Hi,

I'm using mda 0.16.0-dev0 to do some clustering of xtc trajectories with encore.cluster. I got some unexpected results when writing the clusters as xtc files, and I have pinpointed the problem to the way I loaded the trajectory data.
The problem is that the times in the written trajectory does not correspond to the original trajectory (starting from 0 with the wrong offset).

I use u.transfer_to_memory(step=5) before the actual clustering, and I noticed that the time offset (dt) is not updated to reflect the step variable.

Expected behaviour

The gap between frames (dt) should be updated to reflect the step chosen in u.transfer_to_memory() so that the calculated times correspond to the same frames in the original trajectory.

Actual behaviour

The offset (dt) is not updated, and time values do not correspond to frames in the original trajectory.

Code to reproduce the behaviour

import numpy as np
import MDAnalysis as mda
from MDAnalysisTests.datafiles import TPR, XTC
	
u = mda.Universe(TPR, XTC)

times = [ts.time for ts in u.trajectory]

u.transfer_to_memory(step=2)
times2 = [ts.time for ts in u.trajectory]

print times[::2]
# [0.0, 200.00001525878906, 400.0000305175781, 600.0, 800.0000610351562]
print times2
# [0.0, 100.00000762939453, 200.00001525878906, 300.0000228881836, 400.0000305175781]

assert np.allclose(times[::2], times2)
        

Current version of MDAnalysis:

0.16.0-dev0

[richardjgowers]

Hi Jon, thanks for the report.

I think fixing this might just be a case of fixing this line so it's dt = ts.dt * step.

But it also ties in with #1041 where we're not really copying the entire Timestep into the memory reader, but instead just coordinates. A more robust fix for this (and future corner cases) might be to copy ts.time for each frame slurped into the memory reader.

[orbeckst]

In principle we should have something like a time auxiliary.

Maybe to a quick fix first and then work on a more comprehensive solution along the lines of #1041

[kaplajon]

Hi,

Thanks for the feedback. I'll try the quickfix to see of I can get something useful out. :)

/Jon

[kaplajon]

Quickfix tested, works as expected. Thanks!

[orbeckst]

@kaplajon do you want to submit a small PR with the fix?

We'll then also ask for a test for the broken functionality so it's going to be a tiny bit more work but you get to be an author/contributor on the project.

We can then later tackle the bigger problems along the lines of #1041 but the new test would already be extremely valuable for that endeavor.

[kaplajon]

@orbeckst , I could at least try (but probably not until the end of the week). I'm not very familiar with team development and all the whereabouts of github, but I guess it's never too late to learn.

Should I build the test case around XTC and TPR? It seems that the standard case in MDA is DCD trajectories.

[orbeckst]

To get started, have a look at the docs for Contributing Code – it talks about cloning the repo and making a pull request (PR).

If you can recreate the test with PSF and DCD then that's even better (because it is faster). The test could look like the following. This code would be in testsuite/MDAnalysisTests/coordinates/test_memory.py — the standard memory rader test is a bit involved and rather than figuring out immediately how to use the existing test, just do something simple:

from numpy.testing import assert_array_almost_equal
def test_memory_step():
   u = mda.Universe(PSF, DCD)
   times = [ts.time for ts in u.trajectory]
   u.transfer_to_memory(step=2)
   times2 = [ts.time for ts in u.trajectory]

   assert array_almost_equal(times[::2], times2)

As part of the discussion on the PR we can then figure out how to streamline your test. But this can be done as part of the process of getting the code merged and other people can also chime in; all of this is much easier if we have code in front of us.

This test test_memory_step should fail without your patch.

I would start a PR with the test. Check that it fails.

Then implement the patch #1310 (comment)

I think fixing this might just be a case of fixing this line so it's dt = ts.dt * step.

and check that your test passes.

You can submit a PR even if you're still working on the code. New commits will just add to the PR.

Ask if you have questions!

[kaplajon]

Thanks,
there is a test_slicing_step_without_start_stop() in core/test_universe.py, that calls transfer_to_memory. Could that be extended to include another assert call, or could my test go into the same class?

Btw, what is the reason some of the tests use PDB_small and others PSF? E.g. all test_slicing...() functions use PDB_small, while test_frame_interval_convention() uses PSF. Just curious.

[richardjgowers]

Add another function / method to the class rather than making a test do two things. Ideally a test failing points to a single cause. The use of test files is a little haphazard, try and follow the existing patterns within the class if you can. Thanks for helping out!

…

On Tue, 25 Apr 2017, 2:16 p.m. Jon Kapla, ***@***.***> wrote: Thanks, there is a test_slicing_step_without_start_stop() in core/test_universe.py <https://github.com/MDAnalysis/mdanalysis/blob/develop/testsuite/MDAnalysisTests/core/test_universe.py#L397>, that calls transfer_to_memory. Could that be extended to include another assert call, or could my test go into the same class? Btw, what is the reason some of the tests use PDB_small and others PSF? E.g. all test_slicing...() functions use PDB_small, while test_frame_interval_convention() uses PSF. Just curious. — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#1310 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AI0jB7tV8iTcvXQTbQrAY5xz_kaYftNmks5rzfJCgaJpZM4M9GGP> .

[kaplajon]

Thanks, I'm on to it. One thing though: If I do dt=self.trajectory.ts.dt * step in core/universe.py: transfer_to_memory(), step has to default to 1, and right now it defaults to None. Should I do step=1 in the def? Would that break anything else?

[kain88-de]

None is equivalent to a 1 here. Please keep it like this as this follows standard python API's as it is.

[kaplajon]

Well, thats what I thought, but I get failed TestInMemoryUniverse tests with TypeError: unsupported operand type(s) for *: 'float' and 'NoneType'. Would it be better to just test if step is None and set it to one if it is?

[kain88-de]

Could you please open a PR with your code as it is right now. That makes it easier for us to see why the error occurs. It doesn't matter if the code doesn't work completely yet but giving advice is just much easier for us that way.

[kaplajon]

Sure, I'll push my changes tonight and open the PR. The status right now is that it should work. Den 26 apr. 2017 3:30 em skrev "Max Linke" <[email protected]>: Could you please open a PR with your as it is right now. That makes it easier for us to see why the error occurs. It doesn't matter if the code doesn't work completely yet but giving advice is just much easier for us that way. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#1310 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AIaZhC07wBfiySiO9R_njDRupOCUI-qKks5rz0bkgaJpZM4M9GGP> .