-
Notifications
You must be signed in to change notification settings - Fork 658
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
Slicing fixes #916
Slicing fixes #916
Changes from 26 commits
93b52b0
9a58b1b
8145421
0649737
9e7d439
e9b6586
04dc85d
1983709
cb5852e
1b1d4f2
320b1dc
55910a0
83c7e4f
beffbea
8e24480
d687cc7
6893abd
6b27588
32fce06
68ae897
6f9e57d
e476ee5
d036859
0c268ed
d3d4dad
6833e3a
76ec586
0f62787
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -73,6 +73,7 @@ | |
import numpy as np | ||
import struct | ||
import types | ||
import warnings | ||
|
||
from ..core import flags | ||
from .. import units as mdaunits # use mdaunits instead of units to avoid a clash | ||
|
@@ -85,6 +86,7 @@ | |
from . import dcdtimeseries | ||
|
||
|
||
|
||
class Timestep(base.Timestep): | ||
#: Indices into :attr:`Timestep._unitcell` (``[A, gamma, B, beta, alpha, | ||
#: C]``, provided by the :class:`DCDReader` C code) to pull out | ||
|
@@ -498,22 +500,33 @@ def _read_frame(self, frame): | |
ts.frame = frame | ||
return ts | ||
|
||
def timeseries(self, asel, start=0, stop=-1, skip=1, format='afc'): | ||
def timeseries(self, asel, start=None, stop=None, step=None, skip=None, | ||
format='afc'): | ||
"""Return a subset of coordinate data for an AtomGroup | ||
|
||
:Arguments: | ||
*asel* | ||
:class:`~MDAnalysis.core.AtomGroup.AtomGroup` object | ||
*start, stop, skip* | ||
range of trajectory to access, start and stop are inclusive | ||
*start, stop, step* | ||
A range of the trajectory to access, with start being inclusive | ||
and stop being exclusive. | ||
*format* | ||
the order/shape of the return data array, corresponding | ||
to (a)tom, (f)rame, (c)oordinates all six combinations | ||
of 'a', 'f', 'c' are allowed ie "fac" - return array | ||
where the shape is (frame, number of atoms, | ||
coordinates) | ||
:Deprecated: | ||
*skip* | ||
Skip has been deprecated in favor of the standard keyword step. | ||
""" | ||
start, stop, skip = self.check_slice_indices(start, stop, skip) | ||
if skip is not None: | ||
step = skip | ||
warnings.warn("Skip is deprecated and will be removed in" | ||
"in 1.0. Use step instead.", | ||
category=DeprecationWarning) | ||
|
||
start, stop, step = self.check_slice_indices(start, stop, step) | ||
if len(asel) == 0: | ||
raise NoDataError("Timeseries requires at least one atom to analyze") | ||
if len(format) != 3 and format not in ['afc', 'acf', 'caf', 'cfa', 'fac', 'fca']: | ||
|
@@ -522,26 +535,36 @@ def timeseries(self, asel, start=0, stop=-1, skip=1, format='afc'): | |
# Check if the atom numbers can be grouped for efficiency, then we can read partial buffers | ||
# from trajectory file instead of an entire timestep | ||
# XXX needs to be implemented | ||
return self._read_timeseries(atom_numbers, start, stop, skip, format) | ||
return self._read_timeseries(atom_numbers, start, stop, step, format) | ||
|
||
def correl(self, timeseries, start=0, stop=-1, skip=1): | ||
def correl(self, timeseries, start=None, stop=None, step=None, skip=None): | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Did you test that step in There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. I have not addressed this, thanks. There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. @orbeckst When you are saying 'behaves the way it should', you mean that start stop and step reflect the same indices in a trajectory as numpy array slicing? If this is what you mean, it definitely does not, I can go in and change the pyrex code to do that in the same fashion that I did for The second interpretation of what you said, just handing the new defaults correctly, also hasn't been addressed, and is easier to test. There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Yes, all of these functions should behave consistently, eg in the same way as in timeseries and general trajectory slicing. If we are doing the cleanup everywhere else then we need to do it here, to.
Oliver Beckstein * [email protected] There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Alright. This PR was a bit of a pandora's box I guess 😃 There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more.
More like a hornets' nest. But yes, you're right. Much more than "yeah, just fix the docs". Oliver Beckstein * [email protected] |
||
"""Populate a TimeseriesCollection object with timeseries computed from the trajectory | ||
|
||
:Arguments: | ||
*timeseries* | ||
:class:`MDAnalysis.core.Timeseries.TimeseriesCollection` | ||
*start, stop, skip* | ||
subset of trajectory to use, with start and stop being inclusive | ||
*start, stop, step* | ||
A subset of the trajectory to use, with start being inclusive | ||
and stop being exclusive. | ||
:Deprecated: | ||
*skip* | ||
Skip has been deprecated in favor of the standard keyword step. | ||
""" | ||
start, stop, skip = self.check_slice_indices(start, stop, skip) | ||
if skip is not None: | ||
step = skip | ||
warnings.warn("Skip is deprecated and will be removed in" | ||
"in 1.0. Use step instead.", | ||
category=DeprecationWarning) | ||
|
||
start, stop, step = self.check_slice_indices(start, stop, step) | ||
atomlist = timeseries._getAtomList() | ||
format = timeseries._getFormat() | ||
lowerb, upperb = timeseries._getBounds() | ||
sizedata = timeseries._getDataSize() | ||
atomcounts = timeseries._getAtomCounts() | ||
auxdata = timeseries._getAuxData() | ||
return self._read_timecorrel(atomlist, atomcounts, format, auxdata, | ||
sizedata, lowerb, upperb, start, stop, skip) | ||
sizedata, lowerb, upperb, start, stop, step) | ||
|
||
def close(self): | ||
if self.dcdfile is not None: | ||
|
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -59,7 +59,7 @@ cdef extern from "correl.h": | |
|
||
import numpy as np | ||
|
||
def __read_timecorrel(object self, object atoms, object atomcounts, object format, object auxdata, int sizedata, int lowerb, int upperb, int start, int stop, int skip): | ||
def __read_timecorrel(object self, object atoms, object atomcounts, object format, object auxdata, int sizedata, int lowerb, int upperb, int start, int stop, int step): | ||
cdef dcdhandle* dcd | ||
cdef numpy.ndarray atomlist, atomcountslist, auxlist | ||
cdef numpy.ndarray data, temp | ||
|
@@ -69,8 +69,9 @@ def __read_timecorrel(object self, object atoms, object atomcounts, object forma | |
|
||
dcd = <dcdhandle*>PyCObject_AsVoidPtr(self._dcd_C_ptr) | ||
cdef int n_frames | ||
if (stop == -1): stop = dcd.nsets | ||
n_frames = (stop-start+1) / skip | ||
|
||
n_frames = ((stop-start) / step); | ||
if (stop-start) % step > 0 : n_frames += 1; | ||
cdef int numdata | ||
numdata = len(format) | ||
if numdata==0: | ||
|
@@ -109,18 +110,31 @@ def __read_timecorrel(object self, object atoms, object atomcounts, object forma | |
cdef int index, numskip | ||
cdef int i, j | ||
cdef float unitcell[6] | ||
cdef int remaining_frames = stop-start | ||
numskip = 0 | ||
for i from 0 <= i < n_frames: | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. This comment said "We should check for fixed atoms here but we have not implemented it". I don't think that you are addressing fixed atoms anywhere (it's an obscure feature in CHARMM DCD files where only the first frame contains the positions of the frozen atoms and all later frames do not have these positions, which are unchanging.) Therefore, please leave this comment in. |
||
if (skip > 1): | ||
if (step > 1 and i > 0): | ||
# Check if we have fixed atoms | ||
# XXX not done | ||
numskip = skip - (dcd.setsread % skip) - 1 | ||
# Figure out how many steps to step over, if step = n, np array | ||
# slicing treats this as skip over n-1, read the nth. | ||
numskip = step - 1 | ||
# If the number to skip is greater than the number of frames left | ||
# to be jumped over, just take one more step to reflect np slicing | ||
# if there is a remainder, guaranteed to have at least one more | ||
# frame. | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. This is good, thanks. |
||
if(remaining_frames < numskip): | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. What does this do? Can you add a comment to explain, please. There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. I will add a comment, thanks. There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Thanks – I am interested in the summary of thoughts that go into these few lines (not "subtract one from step and call it numskip"....). Clearly, this is one of the key points where you put thought into so for others to quickly catch up it is important to get them on the right track when reading the code. |
||
numskip = 1 | ||
|
||
rc = skip_dcdstep(dcd.fd, dcd.natoms, dcd.nfixed, dcd.charmm, numskip) | ||
if (rc < 0): | ||
raise IOError("Error skipping frame from DCD file") | ||
dcd.setsread = dcd.setsread + numskip | ||
# on first iteration, numskip = 0, first set is always read. | ||
dcd.setsread += numskip | ||
rc = read_dcdsubset(dcd.fd, dcd.natoms, lowerb, upperb, tempX, tempY, tempZ, unitcell, dcd.nfixed, dcd.first, dcd.freeind, dcd.fixedcoords, dcd.reverse, dcd.charmm) | ||
dcd.first=0 | ||
dcd.setsread = dcd.setsread + 1 | ||
dcd.first = 0 | ||
dcd.setsread += 1 | ||
remaining_frames = stop - dcd.setsread | ||
if (rc < 0): | ||
raise IOError("Error reading frame from DCD file") | ||
# Copy into data array based on format | ||
|
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Here add a deprecation warning, skip to be removed in 1.0.
Also add a
.. versionchanged:: 0.16.0
that explains howskip
is nowstep
and behaves differently; you can add a link to issue #818.