Skip to content
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

Merged
merged 28 commits into from
Aug 2, 2016
Merged
Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
28 commits
Select commit Hold shift + click to select a range
93b52b0
Reformatting commit
jdetle Jul 29, 2016
9a58b1b
RMSF fix
jdetle Jul 27, 2016
8145421
Fixed DCD defaults and documentation.
jdetle Jul 27, 2016
0649737
Fixed analysis docs
jdetle Jul 27, 2016
9e7d439
Check_slice index fix
jdetle Jul 27, 2016
e9b6586
DCD fix
jdetle Jul 27, 2016
04dc85d
Changed RMSFarray to reflect correct data
jdetle Jul 27, 2016
1983709
Updated changelog [skip ci]
jdetle Jul 27, 2016
cb5852e
Added test for slicing, slowly changing DCD source for timeseries sli…
jdetle Jul 27, 2016
1b1d4f2
DCD.py changes
jdetle Jul 27, 2016
320b1dc
Everything working.
jdetle Jul 27, 2016
55910a0
Change to ++
jdetle Jul 28, 2016
83c7e4f
Elaborated on docs
jdetle Jul 28, 2016
beffbea
Added generator but it appears to not be working, please help.
jdetle Jul 28, 2016
8e24480
Deprecated skip for step
jdetle Jul 28, 2016
d687cc7
Changed to warnings.warn() because decorator wasnt working
jdetle Jul 29, 2016
6893abd
Fixed silly correl bug
jdetle Jul 29, 2016
6b27588
Fixed issue caused by integer division, when stop-start / nframes has…
jdetle Jul 29, 2016
32fce06
Reduced number of tests by a signifcant margin, removed test of test
jdetle Jul 30, 2016
68ae897
Fixed incomplete logic
jdetle Jul 30, 2016
6f9e57d
Changes to DCD correl functions
jdetle Jul 31, 2016
e476ee5
Changed to logic to take advantage of new default, simplify code.
jdetle Jul 31, 2016
d036859
I made a silly mistake in Timeseries
jdetle Jul 31, 2016
0c268ed
More testing
jdetle Jul 31, 2016
d3d4dad
Added comments.
jdetle Jul 31, 2016
6833e3a
Undo comment removal [skip ci]
jdetle Jul 31, 2016
76ec586
Added tests that cause failures
jdetle Aug 1, 2016
0f62787
Added known failures
jdetle Aug 1, 2016
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
25 changes: 17 additions & 8 deletions package/MDAnalysis/coordinates/dcdtimeseries.pyx
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -69,8 +69,10 @@ 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

if stop == dcd.nsets: stop -= 1;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

add a comment to explain the stop gymnastics

n_frames = ((stop-start) / step) + 1;
if (stop-start) % step == 0: n_frames -= 1;
cdef int numdata
numdata = len(format)
if numdata==0:
Expand Down Expand Up @@ -109,18 +111,25 @@ 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
Copy link
Member

Choose a reason for hiding this comment

The 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.

print n_frames
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

remove print (note: I am commenting on individual commits, if this has been removed later, please ignore)

for i from 0 <= i < n_frames:
if (skip > 1):
if (step > 1 and i > 0):
# Check if we have fixed atoms
# XXX not done
numskip = skip - (dcd.setsread % skip) - 1
numskip = step - 1
if(remaining_frames < numskip):
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What does this do?

Can you add a comment to explain, please.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I will add a comment, thanks.

Copy link
Member

@orbeckst orbeckst Jul 31, 2016

Choose a reason for hiding this comment

The 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
dcd.setsread = dcd.setsread + numskip
Copy link
Member

@orbeckst orbeckst Jul 31, 2016

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Are you sure that this should now be one level out, i.e. not inside the if (step > 1 and skip >1) statement any more?

Copy link
Contributor Author

@jdetle jdetle Jul 31, 2016

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes.
Edit: Actually, I am not sure. Let me add some more tests.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Everything is fine, current tests account for this, the logic is sound too.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ok. Btw, can we write this aesthetically more pleasingly as

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
Expand Down
4 changes: 2 additions & 2 deletions package/MDAnalysis/coordinates/src/dcd.c
Original file line number Diff line number Diff line change
Expand Up @@ -467,8 +467,8 @@ __read_timeseries(PyObject *self, PyObject *args)
Py_DECREF(temp);


if (stop == dcd->nsets) { stop = dcd->nsets -1; }
n_frames = ((stop-start) / step)+ 1;
if (stop == dcd->nsets) { stop = dcd->nsets - 1; }
n_frames = ((stop-start) / step) + 1;
if ((stop-start) % step == 0) { n_frames--;}
//n_frames = dcd->nsets / skip;
n_atoms = PyList_Size((PyObject*)atoms);
Expand Down
25 changes: 20 additions & 5 deletions package/MDAnalysis/core/Timeseries.py
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding:utf-8 -*-
# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4
# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4
#
# MDAnalysis --- http://www.MDAnalysis.org
# Copyright (c) 2006-2015 Naveen Michaud-Agrawal, Elizabeth J. Denning, Oliver Beckstein
Expand Down Expand Up @@ -52,6 +52,7 @@
.. autoclass:: WaterDipole

"""
import warnings

from . import AtomGroup

Expand Down Expand Up @@ -104,16 +105,30 @@ def clear(self):
'''clear the timeseries collection'''
self.timeseries = []

def compute(self, trj, start=0, stop=-1, skip=1):
def compute(self, trj, start=None, stop=None, skip=None, step=None):
'''Iterate through the trajectory *trj* and compute the time series.

*trj*
dcd trajectory object (i.e. :attr:`Universe.trajectory`)
*start*
First frame of trajectory to analyse, Default: None becomes 0.
*stop*
Frame index to stop analysis. Default: None becomes
n_frames. Iteration stops *before* this frame number,
which means that the trajectory would be read until the end.
*step*
Step between frames to analyse, Default: None becomes 1.
*Deprecated*
Skip is deprecated in favor of step.

*start, stop, skip*
Frames to calculate parts of the trajectory. It is important to
note that *start* and *stop* are inclusive
'''
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 = trj.check_slice_indices(start, stop, skip)
self.data = trj.correl(self, start, stop, skip)
# Now remap the timeseries data to each timeseries
typestr = "|f8"
Expand Down
9 changes: 6 additions & 3 deletions testsuite/MDAnalysisTests/coordinates/test_dcd.py
Original file line number Diff line number Diff line change
Expand Up @@ -432,9 +432,10 @@ def test_coordinates(self):
class TestDCDCorrel(TestCase):
def setUp(self):
# Note: setUp is executed for *every* test !
super(TestDCDCorrel, self).setUp()
import MDAnalysis.core.Timeseries as TS

self.universe = mda.Universe(PSF, DCD)
self.dcd = self.universe.trajectory
self.ts = self.universe.coord
self.collection = TS.TimeseriesCollection()
C = self.collection
all = self.universe.atoms
Expand All @@ -458,7 +459,9 @@ def setUp(self):

def tearDown(self):
del self.collection
super(TestDCDCorrel, self).tearDown()
del self.universe
del self.dcd
del self.ts

def test_correl(self):
assert_equal(len(self.collection), 10, "Correl: len(collection)")
Expand Down