diff --git a/package/CHANGELOG b/package/CHANGELOG index 30a9ed4a2f5..59c390d8b45 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -36,7 +36,8 @@ Fixes * Display of Deprecation warnings doesn't affect other modules anymore (Issue #754) * Changed nframes to n_frames in analysis modules for consistency (Issue #890) * fixed incompatibility with newer matplotlib in analysis.hole - + * Fixed modules that improperly documented and/or used frame slicing + defaults (#914) Changes * Added protected variable _frame_index to to keep track of frame iteration number in AnalysisBase diff --git a/package/MDAnalysis/analysis/contacts.py b/package/MDAnalysis/analysis/contacts.py index b8b5bbf00ff..5a1f8cfb656 100644 --- a/package/MDAnalysis/analysis/contacts.py +++ b/package/MDAnalysis/analysis/contacts.py @@ -396,11 +396,13 @@ def __init__(self, u, selection, refgroup, method="hard_cut", radius=4.5, dictionary of additional kwargs passed to `method`. Check respective functions for reasonable values. start : int, optional - First frame of trajectory to analyse, Default: 0 + First frame of trajectory to analyse, Default: None becomes 0. stop : int, optional - Last frame of trajectory to analyse, Default: -1 + 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 : int, optional - Step between frames to analyse, Default: 1 + Step between frames to analyse, Default: None becomes 1. """ self.u = u diff --git a/package/MDAnalysis/analysis/diffusionmap.py b/package/MDAnalysis/analysis/diffusionmap.py index a4337b50fcc..a6300212cc0 100644 --- a/package/MDAnalysis/analysis/diffusionmap.py +++ b/package/MDAnalysis/analysis/diffusionmap.py @@ -227,11 +227,13 @@ def __init__(self, u, select='all', metric=rmsd, cutoff=1E0-5, weights : array, optional Weights to be given to coordinates for metric calculation start : int, optional - First frame of trajectory to analyse, Default: 0 + First frame of trajectory to analyse, Default: None becomes 0. stop : int, optional - Last frame of trajectory to analyse, Default: -1 + 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 : int, optional - Step between frames to analyse, Default: 1 + Step between frames to analyse, Default: None becomes 1. """ self._u = u traj = self._u.trajectory diff --git a/package/MDAnalysis/analysis/polymer.py b/package/MDAnalysis/analysis/polymer.py index 53e3044f40d..ce826f80b2d 100644 --- a/package/MDAnalysis/analysis/polymer.py +++ b/package/MDAnalysis/analysis/polymer.py @@ -59,11 +59,13 @@ def __init__(self, atomgroups, **kwargs): List of atomgroups. Each atomgroup should represent a single polymer chain, ordered in the correct order. start : int, optional - First frame of trajectory to analyse, Default: 0 + First frame of trajectory to analyse, Default: None becomes 0. stop : int, optional - Last frame of trajectory to analyse, Default: -1 + Last frame of trajectory to analyse, Default: None becomes + n_frames. step : int, optional - Step between frames to analyse, Default: 1 + Frame index to stop analysis. Default: None becomes + n_frames. Iteration stops *before* this frame number. """ super(PersistenceLength, self).__init__( atomgroups[0].universe.trajectory, **kwargs) diff --git a/package/MDAnalysis/analysis/rms.py b/package/MDAnalysis/analysis/rms.py index f644bf63a6e..16ed9ae815b 100644 --- a/package/MDAnalysis/analysis/rms.py +++ b/package/MDAnalysis/analysis/rms.py @@ -556,7 +556,7 @@ def __init__(self, atomgroup): self.atomgroup = atomgroup self._rmsf = None - def run(self, start=0, stop=-1, step=1, progout=10, quiet=False): + def run(self, start=None, stop=None, step=None, progout=10, quiet=False): """Calculate RMSF of given atoms across a trajectory. This method implements an algorithm for computing sums of squares while @@ -565,11 +565,13 @@ def run(self, start=0, stop=-1, step=1, progout=10, quiet=False): Parameters ---------- start : int (optional) - starting frame [0] + starting frame, default None becomes 0. stop : int (optional) - stopping frame [-1] + 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 : int (optional) - step between frames [1] + step between frames, default None becomes 1. progout : int (optional) number of frames to iterate through between updates to progress output; ``None`` for no updates [10] @@ -583,6 +585,8 @@ def run(self, start=0, stop=-1, step=1, progout=10, quiet=False): Calculating Corrected Sums of Squares and Products." Technometrics 4(3):419-420. """ + traj = self.atomgroup.universe.trajectory + start, stop, step = traj.check_slice_indices(start, stop, step) sumsquares = np.zeros((self.atomgroup.n_atoms, 3)) means = np.array(sumsquares) diff --git a/package/MDAnalysis/coordinates/DCD.py b/package/MDAnalysis/coordinates/DCD.py index a61c2db42e8..e121a098295 100644 --- a/package/MDAnalysis/coordinates/DCD.py +++ b/package/MDAnalysis/coordinates/DCD.py @@ -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,18 +535,28 @@ 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): """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() @@ -541,7 +564,7 @@ def correl(self, timeseries, start=0, stop=-1, skip=1): 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: diff --git a/package/MDAnalysis/coordinates/base.py b/package/MDAnalysis/coordinates/base.py index 6723237098d..16eb08bfba8 100644 --- a/package/MDAnalysis/coordinates/base.py +++ b/package/MDAnalysis/coordinates/base.py @@ -1230,7 +1230,7 @@ def check_slice_indices(self, start, stop, step): ---------- start, stop, step : int or None Values representing the slice indices. - Can use `None` to use defaults of (0, -1, and 1) + Can use `None` to use defaults of (0, n_frames, and 1) respectively. Returns diff --git a/package/MDAnalysis/coordinates/dcdtimeseries.pyx b/package/MDAnalysis/coordinates/dcdtimeseries.pyx index c60e20589c8..1f1d4f585e8 100644 --- a/package/MDAnalysis/coordinates/dcdtimeseries.pyx +++ b/package/MDAnalysis/coordinates/dcdtimeseries.pyx @@ -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 = 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: - 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. + if(remaining_frames < numskip): + 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 diff --git a/package/MDAnalysis/coordinates/src/dcd.c b/package/MDAnalysis/coordinates/src/dcd.c index 4b8b761a4af..d0d3956dbe9 100644 --- a/package/MDAnalysis/coordinates/src/dcd.c +++ b/package/MDAnalysis/coordinates/src/dcd.c @@ -63,12 +63,12 @@ __write_dcd_header(PyObject *self, PyObject *args) charmm |= DCD_HAS_EXTRA_BLOCK; if (! self) { - /* we were in fact called as a module function, try to retrieve + /* we were in fact called as a module function, try to retrieve a matching object from args */ if( !PyArg_ParseTuple(args, "Oi|iids", &self, &natoms, &istart, &nsavc, &delta, &remarks) ) return NULL; } else { - /* we were obviously called as an object method so args should + /* we were obviously called as an object method so args should only have the int value. */ if( !PyArg_ParseTuple(args, "i|iids", &natoms, &istart, &nsavc, &delta, &remarks) ) return NULL; @@ -80,23 +80,23 @@ __write_dcd_header(PyObject *self, PyObject *args) PyErr_SetString(PyExc_AttributeError, "dcdfile is not an attribute"); return NULL; } - + if ((temp = PyObject_GetAttrString(self, "dcdfile")) == NULL) { // This gives me a New Reference // Raise exception PyErr_SetString(PyExc_AttributeError, "dcdfile is not an attribute"); return NULL; } - + if (!PyFile_CheckExact(temp)) { // Raise exception PyErr_SetString(PyExc_TypeError, "dcdfile does not refer to a file object"); Py_DECREF(temp); return NULL; } - fd = fileno(PyFile_AsFile(temp)); + fd = fileno(PyFile_AsFile(temp)); // No longer need the reference to temp Py_DECREF(temp); - + dcd = (dcdhandle *)malloc(sizeof(dcdhandle)); memset(dcd, 0, sizeof(dcdhandle)); dcd->fd = fd; @@ -122,7 +122,7 @@ __write_dcd_header(PyObject *self, PyObject *args) return NULL; } Py_DECREF(temp); - + // For debugging purposes temp = PyBuffer_FromMemory(dcd, sizeof(dcdhandle)); // Creates a New Reference if (PyObject_SetAttrString(self, "_dcd_C_str", temp) == -1) { @@ -132,7 +132,7 @@ __write_dcd_header(PyObject *self, PyObject *args) return NULL; } Py_DECREF(temp); - + Py_INCREF(Py_None); return Py_None; } @@ -146,14 +146,14 @@ __write_next_frame(PyObject *self, PyObject *args) int rc, curstep; float* uc_array; double unitcell[6]; - + if (!self) { - /* we were in fact called as a module function, try to retrieve + /* we were in fact called as a module function, try to retrieve a matching object from args */ if( !PyArg_ParseTuple(args, "OO!O!O!O!", &self, &PyArray_Type, &x, &PyArray_Type, &y, &PyArray_Type, &z, &PyArray_Type, &uc) ) return NULL; } else { - /* we were obviously called as an object method so args should + /* we were obviously called as an object method so args should only have the int value. */ if( !PyArg_ParseTuple(args, "O!O!O!O!", &PyArray_Type, &x, &PyArray_Type, &y, &PyArray_Type, &z, &PyArray_Type, &uc) ) return NULL; @@ -198,29 +198,29 @@ __finish_dcd_write(PyObject *self, PyObject *args) //PyObject* temp; //dcdhandle *dcd; - if (! self) { - /* we were in fact called as a module function, try to retrieve + if (! self) { + /* we were in fact called as a module function, try to retrieve a matching object from args */ if( !PyArg_ParseTuple(args, "O", &self) ) return NULL; - } else { - /* we were obviously called as an object method so args should + } else { + /* we were obviously called as an object method so args should only have the int value. */ if( !PyArg_ParseTuple(args, "") ) - return NULL; - } + return NULL; + } /*if ( !PyObject_HasAttrString(self, "_dcd_C_ptr") ) { // Raise exception PyErr_SetString(PyExc_AttributeError, "_dcd_C_ptr is not an attribute"); return NULL; } - + if ((temp = PyObject_GetAttrString(self, "_dcd_C_ptr")) == NULL) { // This gives me a New Reference // Raise exception PyErr_SetString(PyExc_AttributeError, "_dcd_C_ptr is not an attribute"); return NULL; - } + } dcd = (dcdhandle*)PyCObject_AsVoidPtr(temp); free(dcd); Py_DECREF(temp);*/ @@ -250,16 +250,16 @@ __read_dcd_header(PyObject *self, PyObject *args) dcdhandle *dcd = NULL; if (! self) { - /* we were in fact called as a module function, try to retrieve - a matching object from args */ - if( !PyArg_ParseTuple(args, "O", &self) ) - return NULL; - } else { - /* we were obviously called as an object method so args should - only have the int value. */ - if( !PyArg_ParseTuple(args, "") ) - return NULL; - } + /* we were in fact called as a module function, try to retrieve + a matching object from args */ + if( !PyArg_ParseTuple(args, "O", &self) ) + return NULL; + } else { + /* we were obviously called as an object method so args should + only have the int value. */ + if( !PyArg_ParseTuple(args, "") ) + return NULL; + } // Get the file object from the class if (!PyObject_HasAttrString(self, "dcdfile")) { @@ -273,7 +273,7 @@ __read_dcd_header(PyObject *self, PyObject *args) PyErr_SetString(PyExc_AttributeError, "dcdfile is not an attribute"); return NULL; } - + if (!PyFile_CheckExact(temp)) { // Raise exception PyErr_SetString(PyExc_TypeError, "dcdfile does not refer to a file object"); @@ -282,7 +282,7 @@ __read_dcd_header(PyObject *self, PyObject *args) fd = fileno(PyFile_AsFile(temp)); // No longer need the reference to temp Py_DECREF(temp); - + dcd = (dcdhandle *)malloc(sizeof(dcdhandle)); memset(dcd, 0, sizeof(dcdhandle)); dcd->fd = fd; @@ -295,7 +295,7 @@ __read_dcd_header(PyObject *self, PyObject *args) goto error; } - /* + /* * Check that the file is big enough to really hold all the frames it claims to have */ { @@ -319,8 +319,8 @@ __read_dcd_header(PyObject *self, PyObject *args) // Size of header dcd->header_size = fio_ftell(dcd->fd); - /* - * It's safe to use ftell, even though ftell returns a long, because the + /* + * It's safe to use ftell, even though ftell returns a long, because the * header size is < 4GB. */ filesize = stbuf.st_size - fio_ftell(dcd->fd) - firstframesize; @@ -384,7 +384,7 @@ __read_dcd_header(PyObject *self, PyObject *args) goto error; } Py_DECREF(temp); - + temp = PyCObject_FromVoidPtr(dcd, NULL); if (temp == NULL) goto error; if (PyObject_SetAttrString(self, "_dcd_C_ptr", temp) == -1) { @@ -404,7 +404,7 @@ __read_dcd_header(PyObject *self, PyObject *args) goto error; } Py_DECREF(temp); - + // For debugging purposes temp = PyBuffer_FromMemory(dcd, sizeof(dcdhandle)); if (temp == NULL) goto error; @@ -427,214 +427,211 @@ __read_dcd_header(PyObject *self, PyObject *args) static PyObject * __read_timeseries(PyObject *self, PyObject *args) { - PyObject *temp = NULL; - PyArrayObject *coord = NULL; - PyListObject *atoms = NULL; - int lowerb = 0, upperb = 0, range=0; - float *tempX = NULL, *tempY = NULL, *tempZ = NULL; - int rc; - int i, j, index; - int n_atoms = 0, n_frames = 0; - int start = 0, stop = -1, skip = 1, numskip = 0; - dcdhandle *dcd = NULL; - int *atomlist = NULL; - npy_intp dimensions[3]; - float unitcell[6]; - const char* format = "afc"; - - if (!self) { - /* we were in fact called as a module function, try to retrieve - a matching object from args */ - if( !PyArg_ParseTuple(args, "OO!|iiis", &self, &PyList_Type, &atoms, &start, &stop, &skip, &format) ) - return NULL; - } else { - /* we were obviously called as an object method so args should - only have the int value. */ - if( !PyArg_ParseTuple(args, "O!|iiis", &PyList_Type, &atoms, &start, &stop, &skip, &format) ) + PyObject *temp = NULL; + PyArrayObject *coord = NULL; + PyListObject *atoms = NULL; + int lowerb = 0, upperb = 0, range=0; + float *tempX = NULL, *tempY = NULL, *tempZ = NULL; + int rc; + int i, j, index; + int n_atoms = 0, n_frames = 0; + + /* Stop = -1 is incorrect and causes an error, look for a fix of this in line + 469 */ + int start = 0, stop = -1, step = 1, numskip = 0, remaining_frames=0; + dcdhandle *dcd = NULL; + int *atomlist = NULL; + npy_intp dimensions[3]; + float unitcell[6]; + const char* format = "afc"; + + if (!self) { + /* we were in fact called as a module function, try to retrieve + a matching object from args */ + if( !PyArg_ParseTuple(args, "OO!|iiis", &self, &PyList_Type, &atoms, &start, &stop, &step, &format) ) return NULL; - } + } else { + /* we were obviously called as an object method so args should + only have the int value. */ + if( !PyArg_ParseTuple(args, "O!|iiis", &PyList_Type, &atoms, &start, &stop, &step, &format) ) + return NULL; + } - if ((temp = PyObject_GetAttrString(self, "_dcd_C_ptr")) == NULL) { // This gives me a New Reference - // Raise exception - PyErr_SetString(PyExc_AttributeError, "_dcd_C_ptr is not an attribute"); - return NULL; - } + if ((temp = PyObject_GetAttrString(self, "_dcd_C_ptr")) == NULL) { // This gives me a New Reference + // Raise exception + PyErr_SetString(PyExc_AttributeError, "_dcd_C_ptr is not an attribute"); + return NULL; + } - dcd = (dcdhandle*)PyCObject_AsVoidPtr(temp); - Py_DECREF(temp); + dcd = (dcdhandle*)PyCObject_AsVoidPtr(temp); + Py_DECREF(temp); - // Assume that start and stop are valid - if (stop == -1) { stop = dcd->nsets; } - n_frames = (stop-start+1) / skip; - //n_frames = dcd->nsets / skip; - n_atoms = PyList_Size((PyObject*)atoms); - if (n_atoms == 0) { - PyErr_SetString(PyExc_Exception, "No atoms passed into _read_timeseries function"); - return NULL; - } - atomlist = (int*)malloc(sizeof(int)*n_atoms); - memset(atomlist, 0, sizeof(int)*n_atoms); - - // Get the atom indexes - for (i=0;i 0) { n_frames++; } + //n_frames = dcd->nsets / skip; + n_atoms = PyList_Size((PyObject*)atoms); + if (n_atoms == 0) { + PyErr_SetString(PyExc_Exception, "No atoms passed into _read_timeseries function"); + return NULL; + } + atomlist = (int*)malloc(sizeof(int)*n_atoms); + memset(atomlist, 0, sizeof(int)*n_atoms); + + // Get the atom indexes + for (i=0;ifd, dcd->header_size, FIO_SEEK_SET); - dcd->setsread = 0; - dcd->first = 1; + lowerb = atomlist[0]; upperb = atomlist[n_atoms-1]; + range = upperb-lowerb+1; - // Jump to starting frame - jump_to_frame(dcd, start); - - // Now read through trajectory and get the atom pos - tempX = (float*)malloc(sizeof(float)*range); - tempY = (float*)malloc(sizeof(float)*range); - tempZ = (float*)malloc(sizeof(float)*range); - if ((tempX == NULL) | (tempY == NULL) || (tempZ == NULL)) { - PyErr_SetString(PyExc_MemoryError, "Can't allocate temporary space for coordinate arrays"); - goto error; - } - - for (i=0;i 1) { - /*if (dcd->first && dcd->nfixed) { - 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; - if (rc < 0) { - // return an exception - PyErr_SetString(PyExc_IOError, "Error reading first frame from DCD file"); - //fprintf(stderr, "read_dcdstep returned %d\n", rc); - return NULL; - //return MOLFILE_ERROR; - } - dcd->setsread++; - }*/ - // Figure out how many steps to skip - numskip = skip - (dcd->setsread % skip) - 1; - rc = skip_dcdstep(dcd->fd, dcd->natoms, dcd->nfixed, dcd->charmm, numskip); - if (rc < 0) { - // return an exception - PyErr_SetString(PyExc_IOError, "Error skipping frame from DCD file"); - goto error; - } - dcd->setsread+=numskip; + // Figure out the format string + if (strncasecmp(format, "afc", 3) == 0) { + dimensions[0] = n_atoms; dimensions[1] = n_frames; dimensions[2] = 3; + } else if (strncasecmp(format, "acf", 3) == 0) { + dimensions[0] = n_atoms; dimensions[1] = 3; dimensions[2] = n_frames; + } else if (strncasecmp(format, "fac", 3) == 0) { + dimensions[0] = n_frames; dimensions[1] = n_atoms; dimensions[2] = 3; + } else if (strncasecmp(format, "fca", 3) == 0) { + dimensions[0] = n_frames; dimensions[1] = 3; dimensions[2] = n_atoms; + } else if (strncasecmp(format, "caf", 3) == 0) { + dimensions[0] = 3; dimensions[1] = n_atoms; dimensions[2] = n_frames; + } else if (strncasecmp(format, "cfa", 3) == 0) { + dimensions[0] = 3; dimensions[1] = n_frames; dimensions[2] = n_atoms; + } + + coord = (PyArrayObject*) PyArray_SimpleNew(3, dimensions, NPY_DOUBLE); + if (coord == NULL) goto error; + + // Reset trajectory + rc = fio_fseek(dcd->fd, dcd->header_size, FIO_SEEK_SET); + dcd->setsread = 0; + dcd->first = 1; + + // Jump to starting frame + jump_to_frame(dcd, start); + + // Now read through trajectory and get the atom pos + tempX = (float*)malloc(sizeof(float)*range); + tempY = (float*)malloc(sizeof(float)*range); + tempZ = (float*)malloc(sizeof(float)*range); + if ((tempX == NULL) | (tempY == NULL) || (tempZ == NULL)) { + PyErr_SetString(PyExc_MemoryError, "Can't allocate temporary space for coordinate arrays"); + goto error; + } + + remaining_frames = stop-start; + for (i=0;i 1 && i>0) { + // Check if we have fixed atoms + // XXX not done + /* 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. + */ + if(remaining_frames < numskip){ + numskip = 1; + } + rc = skip_dcdstep(dcd->fd, dcd->natoms, dcd->nfixed, dcd->charmm, numskip); + if (rc < 0) { + // return an exception + PyErr_SetString(PyExc_IOError, "Error skipping frame from DCD file"); + goto error; + } } + // on first iteration, numskip == 0, first set is always read. + dcd->setsread += numskip; + //now read from subset 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); + unitcell, dcd->nfixed, dcd->first, dcd->freeind, dcd->fixedcoords, + dcd->reverse, dcd->charmm); dcd->first = 0; dcd->setsread++; + remaining_frames = stop - dcd->setsread; if (rc < 0) { - // return an exception - PyErr_SetString(PyExc_IOError, "Error reading frame from DCD file"); - goto error; - } + // return an exception + PyErr_SetString(PyExc_IOError, "Error reading frame from DCD file"); + goto error; + + } + + // Copy into Numeric array only those atoms we are interested in for (j=0;jdata + a*coord->strides[0] + b*coord->strides[1] + c*coord->strides[2]) - */ - if (strncasecmp(format, "afc", 3) == 0) { - *(double*)(coord->data + j*coord->strides[0] + i*coord->strides[1] + 0*coord->strides[2]) = tempX[index]; - *(double*)(coord->data + j*coord->strides[0] + i*coord->strides[1] + 1*coord->strides[2]) = tempY[index]; - *(double*)(coord->data + j*coord->strides[0] + i*coord->strides[1] + 2*coord->strides[2]) = tempZ[index]; - } else - if (strncasecmp(format, "acf", 3) == 0) { - *(double*)(coord->data + j*coord->strides[0] + 0*coord->strides[1] + i*coord->strides[2]) = tempX[index]; - *(double*)(coord->data + j*coord->strides[0] + 1*coord->strides[1] + i*coord->strides[2]) = tempY[index]; - *(double*)(coord->data + j*coord->strides[0] + 2*coord->strides[1] + i*coord->strides[2]) = tempZ[index]; - } else - if (strncasecmp(format, "fac", 3) == 0) { - *(double*)(coord->data + i*coord->strides[0] + j*coord->strides[1] + 0*coord->strides[2]) = tempX[index]; - *(double*)(coord->data + i*coord->strides[0] + j*coord->strides[1] + 1*coord->strides[2]) = tempY[index]; - *(double*)(coord->data + i*coord->strides[0] + j*coord->strides[1] + 2*coord->strides[2]) = tempZ[index]; - } else - if (strncasecmp(format, "fca", 3) == 0) { - *(double*)(coord->data + i*coord->strides[0] + 0*coord->strides[1] + j*coord->strides[2]) = tempX[index]; - *(double*)(coord->data + i*coord->strides[0] + 1*coord->strides[1] + j*coord->strides[2]) = tempY[index]; - *(double*)(coord->data + i*coord->strides[0] + 2*coord->strides[1] + j*coord->strides[2]) = tempZ[index]; - } else - if (strncasecmp(format, "caf", 3) == 0) { - *(double*)(coord->data + 0*coord->strides[0] + j*coord->strides[1] + i*coord->strides[2]) = tempX[index]; - *(double*)(coord->data + 1*coord->strides[0] + j*coord->strides[1] + i*coord->strides[2]) = tempY[index]; - *(double*)(coord->data + 2*coord->strides[0] + j*coord->strides[1] + i*coord->strides[2]) = tempZ[index]; - } else - if (strncasecmp(format, "cfa", 3) == 0) { - *(double*)(coord->data + 0*coord->strides[0] + i*coord->strides[1] + j*coord->strides[2]) = tempX[index]; - *(double*)(coord->data + 1*coord->strides[0] + i*coord->strides[1] + j*coord->strides[2]) = tempY[index]; - *(double*)(coord->data + 2*coord->strides[0] + i*coord->strides[1] + j*coord->strides[2]) = tempZ[index]; - } + index = atomlist[j]-lowerb; + + /* + * coord[a][b][c] = *(float*)(coord->data + a*coord->strides[0] + b*coord->strides[1] + c*coord->strides[2]) + */ + if (strncasecmp(format, "afc", 3) == 0) { + *(double*)(coord->data + j*coord->strides[0] + i*coord->strides[1] + 0*coord->strides[2]) = tempX[index]; + *(double*)(coord->data + j*coord->strides[0] + i*coord->strides[1] + 1*coord->strides[2]) = tempY[index]; + *(double*)(coord->data + j*coord->strides[0] + i*coord->strides[1] + 2*coord->strides[2]) = tempZ[index]; + } else if (strncasecmp(format, "acf", 3) == 0) { + *(double*)(coord->data + j*coord->strides[0] + 0*coord->strides[1] + i*coord->strides[2]) = tempX[index]; + *(double*)(coord->data + j*coord->strides[0] + 1*coord->strides[1] + i*coord->strides[2]) = tempY[index]; + *(double*)(coord->data + j*coord->strides[0] + 2*coord->strides[1] + i*coord->strides[2]) = tempZ[index]; + } else if (strncasecmp(format, "fac", 3) == 0) { + *(double*)(coord->data + i*coord->strides[0] + j*coord->strides[1] + 0*coord->strides[2]) = tempX[index]; + *(double*)(coord->data + i*coord->strides[0] + j*coord->strides[1] + 1*coord->strides[2]) = tempY[index]; + *(double*)(coord->data + i*coord->strides[0] + j*coord->strides[1] + 2*coord->strides[2]) = tempZ[index]; + } else if (strncasecmp(format, "fca", 3) == 0) { + *(double*)(coord->data + i*coord->strides[0] + 0*coord->strides[1] + j*coord->strides[2]) = tempX[index]; + *(double*)(coord->data + i*coord->strides[0] + 1*coord->strides[1] + j*coord->strides[2]) = tempY[index]; + *(double*)(coord->data + i*coord->strides[0] + 2*coord->strides[1] + j*coord->strides[2]) = tempZ[index]; + } else if (strncasecmp(format, "caf", 3) == 0) { + *(double*)(coord->data + 0*coord->strides[0] + j*coord->strides[1] + i*coord->strides[2]) = tempX[index]; + *(double*)(coord->data + 1*coord->strides[0] + j*coord->strides[1] + i*coord->strides[2]) = tempY[index]; + *(double*)(coord->data + 2*coord->strides[0] + j*coord->strides[1] + i*coord->strides[2]) = tempZ[index]; + } else if (strncasecmp(format, "cfa", 3) == 0) { + *(double*)(coord->data + 0*coord->strides[0] + i*coord->strides[1] + j*coord->strides[2]) = tempX[index]; + *(double*)(coord->data + 1*coord->strides[0] + i*coord->strides[1] + j*coord->strides[2]) = tempY[index]; + *(double*)(coord->data + 2*coord->strides[0] + i*coord->strides[1] + j*coord->strides[2]) = tempZ[index]; + } } // Check if we've been interupted by the user if (PyErr_CheckSignals() == 1) goto error; - } - - // Reset trajectory - rc = fio_fseek(dcd->fd, dcd->header_size, FIO_SEEK_SET); - dcd->setsread = 0; - dcd->first = 1; - free(atomlist); - free(tempX); - free(tempY); - free(tempZ); - return PyArray_Return(coord); - - error: - // Reset trajectory - rc = fio_fseek(dcd->fd, dcd->header_size, FIO_SEEK_SET); - dcd->setsread = 0; - dcd->first = 1; - Py_XDECREF(coord); - if (atomlist != NULL) free(atomlist); - if (tempX != NULL) free(tempX); - if (tempY != NULL) free(tempY); - if (tempZ != NULL) free(tempZ); - return NULL; + } + + // Reset trajectory + rc = fio_fseek(dcd->fd, dcd->header_size, FIO_SEEK_SET); + dcd->setsread = 0; + dcd->first = 1; + free(atomlist); + free(tempX); + free(tempY); + free(tempZ); + return PyArray_Return(coord); + + error: + // Reset trajectory + rc = fio_fseek(dcd->fd, dcd->header_size, FIO_SEEK_SET); + dcd->setsread = 0; + dcd->first = 1; + Py_XDECREF(coord); + if (atomlist != NULL) free(atomlist); + if (tempX != NULL) free(tempX); + if (tempY != NULL) free(tempY); + if (tempZ != NULL) free(tempZ); + return NULL; } + static PyObject * __read_next_frame(PyObject *self, PyObject *args) { @@ -645,16 +642,16 @@ __read_next_frame(PyObject *self, PyObject *args) int rc,numskip; float* unitcell; float alpha, beta, gamma; - + if (!self) { - /* we were in fact called as a module function, try to retrieve + /* we were in fact called as a module function, try to retrieve a matching object from args */ if( !PyArg_ParseTuple(args, "OO!O!O!O!|i", &self, &PyArray_Type, &x, &PyArray_Type, &y, &PyArray_Type, &z, &PyArray_Type, &uc, &skip) ) - return NULL; + return NULL; } else { - /* we were obviously called as an object method so args should + /* we were obviously called as an object method so args should only have the int value. */ - if( !PyArg_ParseTuple(args, "O!O!O!O!|i", &PyArray_Type, &x, &PyArray_Type, &y, &PyArray_Type, &z, &PyArray_Type, &uc, &skip) ) + if( !PyArg_ParseTuple(args, "O!O!O!O!|i", &PyArray_Type, &x, &PyArray_Type, &y, &PyArray_Type, &z, &PyArray_Type, &uc, &skip) ) return NULL; } @@ -669,7 +666,7 @@ __read_next_frame(PyObject *self, PyObject *args) unitcell = (float*) uc->data; unitcell[0] = unitcell[2] = unitcell[5] = 0.0f; - unitcell[1] = unitcell[3] = unitcell[4] = 90.0f; + unitcell[1] = unitcell[3] = unitcell[4] = 90.0f; /* Check for EOF here; that way all EOF's encountered later must be errors */ if (dcd->setsread == dcd->nsets) { @@ -722,7 +719,7 @@ __read_next_frame(PyObject *self, PyObject *args) //fprintf(stderr, "read_dcdstep returned %d\n", rc); return NULL; //return MOLFILE_ERROR; - } + } if (unitcell[1] >= -1.0 && unitcell[1] <= 1.0 && unitcell[3] >= -1.0 && unitcell[3] <= 1.0 && @@ -746,7 +743,7 @@ __read_next_frame(PyObject *self, PyObject *args) unitcell[4] = alpha; unitcell[3] = beta; unitcell[1] = gamma; - + // Return the frame read temp = Py_BuildValue("i", dcd->setsread); return temp; @@ -759,17 +756,17 @@ __finish_dcd_read(PyObject *self, PyObject *args) dcdhandle *dcd; if (! self) { - /* we were in fact called as a module function, try to retrieve - a matching object from args */ + /* we were in fact called as a module function, try to retrieve + a matching object from args */ if( !PyArg_ParseTuple(args, "O", &self) ) - return NULL; - } else { - /* we were obviously called as an object method so args should - only have the int value. */ + return NULL; + } else { + /* we were obviously called as an object method so args should + only have the int value. */ if( !PyArg_ParseTuple(args, "") ) - return NULL; + return NULL; } - + if ( !PyObject_HasAttrString(self, "_dcd_C_ptr") ) { // Raise exception PyErr_SetString(PyExc_AttributeError, "_dcd_C_ptr is not an attribute"); @@ -783,13 +780,13 @@ __finish_dcd_read(PyObject *self, PyObject *args) } dcd = (dcdhandle*)PyCObject_AsVoidPtr(temp); - + close_dcd_read(dcd->freeind, dcd->fixedcoords); free(dcd); Py_DECREF(temp); Py_INCREF(Py_None); return Py_None; -} +} int jump_to_frame(dcdhandle *dcd, int frame) { int rc; @@ -828,23 +825,23 @@ __jump_to_frame(PyObject *self, PyObject *args) int frame; if (! self) { - /* we were in fact called as a module function, try to retrieve + /* we were in fact called as a module function, try to retrieve a matching object from args */ if( !PyArg_ParseTuple(args, "Oi", &self, &frame) ) return NULL; } else { - /* we were obviously called as an object method so args should + /* we were obviously called as an object method so args should only have the int value. */ if( !PyArg_ParseTuple(args, "i", &frame) ) return NULL; - } - + } + if ( !PyObject_HasAttrString(self, "_dcd_C_ptr") ) { // Raise exception PyErr_SetString(PyExc_AttributeError, "_dcd_C_ptr is not an attribute"); return NULL; - } - + } + if ((temp = PyObject_GetAttrString(self, "_dcd_C_ptr")) == NULL) { // This gives me a New Reference // Raise exception PyErr_SetString(PyExc_AttributeError, "_dcd_C_ptr is not an attribute"); @@ -852,7 +849,7 @@ __jump_to_frame(PyObject *self, PyObject *args) } dcd = (dcdhandle*)PyCObject_AsVoidPtr(temp); Py_DECREF(temp); - + /*if (frame > dcd->nsets) { PyErr_SetString(PyExc_IndexError, "Invalid frame"); return NULL; @@ -865,7 +862,7 @@ __jump_to_frame(PyObject *self, PyObject *args) extrablocksize = dcd->charmm & DCD_HAS_EXTRA_BLOCK ? 48 + 8 : 0; ndims = dcd->charmm & DCD_HAS_4DIMS ? 4 : 3; firstframesize = (dcd->natoms+2) * ndims * sizeof(float) + extrablocksize; - framesize = (dcd->natoms-dcd->nfixed+2) * ndims * sizeof(float) + framesize = (dcd->natoms-dcd->nfixed+2) * ndims * sizeof(float) + extrablocksize; // Use zero indexing if (frame == 0) { @@ -890,18 +887,18 @@ __reset_dcd_read(PyObject *self, PyObject *args) PyObject* temp; dcdhandle *dcd; int rc; - + if (! self) { - /* we were in fact called as a module function, try to retrieve + /* we were in fact called as a module function, try to retrieve a matching object from args */ if( !PyArg_ParseTuple(args, "O", &self) ) return NULL; } else { - /* we were obviously called as an object method so args should + /* we were obviously called as an object method so args should only have the int value. */ if( !PyArg_ParseTuple(args, "") ) return NULL; - } + } if ( !PyObject_HasAttrString(self, "_dcd_C_ptr") ) { // Raise exception @@ -943,4 +940,3 @@ init_dcdmodule(void) (void) Py_InitModule("_dcdmodule", DCDMethods); import_array(); } - diff --git a/package/MDAnalysis/core/Timeseries.py b/package/MDAnalysis/core/Timeseries.py index 45a3776f341..f73ab6a701b 100644 --- a/package/MDAnalysis/core/Timeseries.py +++ b/package/MDAnalysis/core/Timeseries.py @@ -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 @@ -52,6 +52,7 @@ .. autoclass:: WaterDipole """ +import warnings from . import AtomGroup @@ -104,17 +105,31 @@ 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 ''' - self.data = trj.correl(self, 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 = trj.check_slice_indices(start, stop, step) + self.data = trj.correl(self, start, stop, step) # Now remap the timeseries data to each timeseries typestr = "|f8" start = 0 diff --git a/testsuite/MDAnalysisTests/coordinates/test_dcd.py b/testsuite/MDAnalysisTests/coordinates/test_dcd.py index 0e901ee87fb..f42ec25dc68 100644 --- a/testsuite/MDAnalysisTests/coordinates/test_dcd.py +++ b/testsuite/MDAnalysisTests/coordinates/test_dcd.py @@ -14,7 +14,7 @@ RefNAMDtriclinicDCD) from MDAnalysisTests.coordinates.base import BaseTimestepTest from MDAnalysisTests import module_not_found, tempdir - +from MDAnalysisTests.plugins.knownfailure import knownfailure @attr('issue') def TestDCD_Issue32(): @@ -23,18 +23,6 @@ def TestDCD_Issue32(): assert_raises(IOError, mda.Universe, PSF, DCD_empty) -class _TestDCD(TestCase): - def setUp(self): - self.universe = mda.Universe(PSF, DCD) - self.dcd = self.universe.trajectory - self.ts = self.universe.coord - - def tearDown(self): - del self.universe - del self.dcd - del self.ts - - class TestDCDReaderClass(TestCase): def test_with_statement(self): from MDAnalysis.coordinates.DCD import DCDReader @@ -55,7 +43,17 @@ def test_with_statement(self): err_msg="with_statement: DCDReader does not read all frames") -class TestDCDReader(_TestDCD): +class TestDCDReader(object): + def setUp(self): + self.universe = mda.Universe(PSF, DCD) + self.dcd = self.universe.trajectory + self.ts = self.universe.coord + + def tearDown(self): + del self.universe + del self.dcd + del self.ts + def test_rewind_dcd(self): self.dcd.rewind() assert_equal(self.ts.frame, 0, "rewinding to frame 0") @@ -137,6 +135,33 @@ def test_volume(self): err_msg="wrong volume for unitcell (no unitcell " "in DCD so this should be 0)") + def test_timeseries_slicing(self): + # check that slicing behaves correctly + # should before issue #914 resolved + x = [(0, 1, 1), (1,1,1), (1, 2, 1), (1, 2, 2), (1, 4, 2), (1, 4, 4), + (0, 5, 5), (3, 5, 1), (None, None, None)] + for start, stop, step in x: + yield self._slice_generation_test, start, stop, step + + def test_backwards_stepping(self): + x = [(4, 0, -1), (5, 0, -2), (5, 0, -4)] + for start, stop, step in x: + yield self._failed_slices_test, start, stop, step + + def _slice_generation_test(self, start, stop, step): + self.u = mda.Universe(PSF, DCD) + ts = self.u.trajectory.timeseries(self.u.atoms) + ts_skip = self.u.trajectory.timeseries(self.u.atoms, start, stop, step) + assert_array_almost_equal(ts[:, start:stop:step,:], ts_skip, 5) + + @knownfailure() + def _failed_slices_test(self, start, stop, step): + self.u = mda.Universe(PSF, DCD) + ts = self.u.trajectory.timeseries(self.u.atoms) + ts_skip = self.u.trajectory.timeseries(self.u.atoms, start, stop, step) + assert_array_almost_equal(ts[:, start:stop:step,:], ts_skip, 5) + + def test_DCDReader_set_dt(dt=100., frame=3): u = mda.Universe(PSF, DCD, dt=dt) assert_almost_equal(u.trajectory[frame].time, frame*dt, @@ -416,14 +441,17 @@ def test_coordinates(self): ts_orig.frame)) -class TestDCDCorrel(_TestDCD): +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() + self.collection_slicing = TS.TimeseriesCollection() C = self.collection + C_step = self.collection_slicing all = self.universe.atoms ca = self.universe.s4AKE.CA ca_termini = mda.core.AtomGroup.AtomGroup([ca[0], ca[-1]]) @@ -440,12 +468,19 @@ def setUp(self): C.addTimeseries(TS.CenterOfGeometry(ca)) # 7 C.addTimeseries(TS.CenterOfMass(all)) # 8 C.addTimeseries(TS.CenterOfGeometry(all)) # 9 + + C_step.addTimeseries(TS.Atom('v', ca_termini)) + # cannot test WaterDipole because there's no water in the test dcd C.compute(self.dcd) + C_step.compute(self.dcd, step=10) def tearDown(self): del self.collection - super(TestDCDCorrel, self).tearDown() + del self.collection_slicing + del self.universe + del self.dcd + del self.ts def test_correl(self): assert_equal(len(self.collection), 10, "Correl: len(collection)") @@ -510,6 +545,11 @@ def test_clear(self): C.clear() assert_equal(len(C), 0, "Correl: clear()") + def test_Atom_slicing(self): + assert_equal(self.collection_slicing[0].shape, (2, 3, 10), + "Correl: Atom positions") + assert_array_almost_equal(self.collection[0][:, :, ::10], + self.collection_slicing[0]) # notes: def compute_correl_references(): diff --git a/testsuite/MDAnalysisTests/data/adk_oplsaa_CA_rmsf.npy b/testsuite/MDAnalysisTests/data/adk_oplsaa_CA_rmsf.npy index d28a8aa8fc1..34d754740a8 100644 Binary files a/testsuite/MDAnalysisTests/data/adk_oplsaa_CA_rmsf.npy and b/testsuite/MDAnalysisTests/data/adk_oplsaa_CA_rmsf.npy differ