Allow analysis classes to run on arbitrary selections of frames

[jbarnoud]

It would be convenient to be able to run some analyses on an arbitrary set of frames. I recently encounter a case where I wanted to build an RSMD matrix with mda.analysis.diffusionmap.DistanceMatrix on a selection of frame I obtained from clustering. For now, the only ways I could find to do that was to either write a new trajectory (on file or in memory), or to hack a SelectionReader. Adding a frames argument to run alonside start, stop, and step (see #1979) would allow to pass an array of frame indices, or a mask, or an already indexed trajectory using #1934. The interface would be the same as for ag.write.

# API as envisioned

u = mda.Universe(TPR, XTC)
ag = u.select_atoms(...)

labels = my_clustering_method(ag)  # returns np.array( [1, 1, 0, 0, 2, ...])
frames_group0 = u.trajectory[labels == 0]

mat0 = MySuperAnalysis(ag).run(frames=frames_group0)
# or
mat0 = MySuperAnalysis(ag).run(frames=(labels == 0))

# Should prevent non-sensical calls such as mixing `frames` with `start`, `stop`, `step`.
mat0 = MySuperAnalysis(ag).run(frames=frames_group0, step=3)  # ValueError

Some analyses rely on the time step to be constant, so the use of the frames must be limited for these analyses. An slices trajectory can be used, but not an indexed one. In these case, a NotImplementedError should be raised. This could be achieved by using a class attribute to switch off arbitrary frame selections:

class MyConstentDTAnalysis(AnalysisBase):
    needs_constant_dt = True

   ...

[orbeckst]

👍 I like the idea of using the new trajectory iterables as a concise way to work with trajectories. I suppose, this is now being called "working with frames".

[orbeckst]

Sidenote: I think this approach

frames_group0 = u.trajectory[labels == 0]

is beautiful.

Would it make sense to add a copy() method to the trajectory (EDIT 👎 )

t2 = u.trajectory.copy()
frames2 = t2[labels == 2]

so that we can have multiple frames iterables that do not interfere with each other, e.g., for cross correlation time analysis... zip(frames1, frames2) and whatever else one finds in itertools.

EDIT: Yes, this is flawed... it is not clear where one AtomGroup should get its coordinates from unless we bring back AtomGroup.coordinates(ts).

[kain88-de]

For this you can use the universe copy method and individual selections for each.

[orbeckst]

We have a Universe.copy() ... apparently so. And we also have trajectory.copy(). I didn't know... (and the docs don't tell me when this was added, so must have been a while.)

[jbarnoud]

Copying only the reader is going to badly mess with the Universe and Atomgroup. Except if it spawn a new Universe, but it is seriously non-intuitive and non trivial to do as a reader does not have knowledge of the universe.

…

On 07/13/2018 07:15 PM, Oliver Beckstein wrote: We have a |Universe.copy()| <https://www.mdanalysis.org/docs/documentation_pages/core/universe.html?highlight=universe%20copy#MDAnalysis.core.universe.Universe.copy> ... apparently so. And we also have |trajectory.copy()|. I didn't know... (and the docs don't tell me when this was added, so must have been a while.) — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#1985 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/ABUWuirRToWgVHvgpZceYmS3FhdvakOmks5uGNXEgaJpZM4VOeEi>.

[orbeckst]

Fair enough ;-)

[kain88-de]

and the docs don't tell me when this was added, so must have been a while

#1538

Nope fairly recent. It's mentioned in the changelog though.

[orbeckst]

Thanks for the insight.

On Jul 13, 2018, at 12:28 PM, Max Linke ***@***.***> wrote #1538 <#1538> Nope fairly recent. It's mentioned in the changelog though.

I tend to look for versionchanged in the docs.

[orbeckst]

@rsexton2 , as discussed, this might be an outline for changing _setup_frames():

from MDAnalysis.tests import datafiles as data
import MDAnalysis as mda

u = mda.Universe(data.TPR, data.XTC)
frames = [0, 2, 4, 5]

def setup_frames(start, stop, step, frames):
    if frames is not None:
         if any([start, stop, step]):
            raise ValueError("start/stop/step cannot be combined with frames")
         slicer = frames
     else:
         slicer = slice(start, stop, step)
     sliced_traj = u.trajectory[slicer]
     n_frames = len(sliced_traj)
     return n_frames, sliced_traj

[IAlibay]

@orbeckst that solution might need care if combined with check_slice_indices?

mdanalysis/package/MDAnalysis/coordinates/base.py

Lines 1694 to 1696 in 943da48

	The returned values `start`, `stop` and `step` give the expected result
	when passed in :func:`range` but gives unexpected behavior when passed
	in a :class:`slice` when ``stop=None`` and ``step=-1``

[orbeckst]

The way that @rsexton2 and I discussed it (which should be reflected in the snippet) is that using frames is exclusive with using any of start, stop, step. That's a lot easier to handle than thinking about what start/stop/step mean for a sliced trajectory. Fundamentally, I see this as two alternative ways of slicing and not something that should be combined, especially because a user who really only wants "every third frame of my custom frame selection" can just do so themselves with frames=myframes[::3].

[orbeckst]

For right now we also did not include the original suggestion to be able to pass a sliced trajectory in frames. Perhaps that should be allowed because we allow writing trajectories with AtomGroup.write(..., frames=...)

ag.write("traj.xtc", frames=ag.universe.trajectory[::2])

where frames can be a sliced trajectory (a FrameIterator) or a array-like of indices or a boolean array of length n_frames, or a slice(). There's something to be said for being consistent but I don't know if there are some odd interactions possible if we pass in a trajectory that belongs to a different universe than the atomgroup or universe that the analysis is supposed to run with.

[orbeckst]

For the purpose of PR #3415 @rsexton2 and I decided to keep things simple and only implement fancy indexing via the frames kwarg and not using a sliced trajectory. If anyone thinks the latter should be included then please say something. Reviews of the PR are also very welcome!