From 93b52b00a65d6acb24caf105ffa1881ea5c5afab Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Fri, 29 Jul 2016 10:53:31 -0700 Subject: [PATCH 01/28] Reformatting commit --- package/MDAnalysis/coordinates/src/dcd.c | 490 +++++++++++------------ 1 file changed, 235 insertions(+), 255 deletions(-) diff --git a/package/MDAnalysis/coordinates/src/dcd.c b/package/MDAnalysis/coordinates/src/dcd.c index 4b8b761a4af..dea24b6884a 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,212 +427,193 @@ __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; + 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; - } - - 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; - } + } 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) ) + return NULL; + } - dcd = (dcdhandle*)PyCObject_AsVoidPtr(temp); - Py_DECREF(temp); + 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); + + // 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;insets; } - 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;ifd, dcd->header_size, FIO_SEEK_SET); - dcd->setsread = 0; - dcd->first = 1; + } 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; + } - // 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 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; } 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++; 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 * @@ -645,16 +626,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 +650,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 +703,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 +727,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 +740,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 +764,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 +809,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 +833,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 +846,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 +871,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 +924,3 @@ init_dcdmodule(void) (void) Py_InitModule("_dcdmodule", DCDMethods); import_array(); } - From 9a58b1b802e7411618370c31158cd185ef7bb0e7 Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Tue, 26 Jul 2016 17:13:27 -0700 Subject: [PATCH 02/28] RMSF fix --- package/MDAnalysis/analysis/rms.py | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/package/MDAnalysis/analysis/rms.py b/package/MDAnalysis/analysis/rms.py index f644bf63a6e..c6170938d1e 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,11 @@ 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] + stopping frame, default None becomes n_frames. 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 +583,9 @@ 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. """ + start, stop, step = self.atomgroup.universe.check_slice_indices(start, + stop, + step) sumsquares = np.zeros((self.atomgroup.n_atoms, 3)) means = np.array(sumsquares) From 814542109f200125d3055f775103b84e203ee9a4 Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Tue, 26 Jul 2016 17:24:12 -0700 Subject: [PATCH 03/28] Fixed DCD defaults and documentation. --- package/MDAnalysis/coordinates/DCD.py | 18 ++++++++++-------- package/MDAnalysis/coordinates/base.py | 2 +- 2 files changed, 11 insertions(+), 9 deletions(-) diff --git a/package/MDAnalysis/coordinates/DCD.py b/package/MDAnalysis/coordinates/DCD.py index a61c2db42e8..47acfd1a740 100644 --- a/package/MDAnalysis/coordinates/DCD.py +++ b/package/MDAnalysis/coordinates/DCD.py @@ -498,14 +498,15 @@ 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, 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 @@ -513,7 +514,7 @@ def timeseries(self, asel, start=0, stop=-1, skip=1, format='afc'): where the shape is (frame, number of atoms, coordinates) """ - start, stop, skip = self.check_slice_indices(start, stop, skip) + start, stop, skip = 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']: @@ -524,16 +525,17 @@ def timeseries(self, asel, start=0, stop=-1, skip=1, format='afc'): # XXX needs to be implemented return self._read_timeseries(atom_numbers, start, stop, skip, format) - def correl(self, timeseries, start=0, stop=-1, skip=1): + def correl(self, timeseries, start=None, stop=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. """ - start, stop, skip = self.check_slice_indices(start, stop, skip) + start, stop, step = self.check_slice_indices(start, stop, step) atomlist = timeseries._getAtomList() format = timeseries._getFormat() lowerb, upperb = timeseries._getBounds() 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 From 0649737a6780e9725c7f1d71297d275cc2bbeb1f Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Tue, 26 Jul 2016 17:35:06 -0700 Subject: [PATCH 04/28] Fixed analysis docs --- package/MDAnalysis/analysis/contacts.py | 7 ++++--- package/MDAnalysis/analysis/diffusionmap.py | 7 ++++--- package/MDAnalysis/analysis/polymer.py | 8 +++++--- package/MDAnalysis/analysis/rms.py | 3 ++- 4 files changed, 15 insertions(+), 10 deletions(-) diff --git a/package/MDAnalysis/analysis/contacts.py b/package/MDAnalysis/analysis/contacts.py index b8b5bbf00ff..6c03381ebcb 100644 --- a/package/MDAnalysis/analysis/contacts.py +++ b/package/MDAnalysis/analysis/contacts.py @@ -396,11 +396,12 @@ 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. 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..b42be7f044f 100644 --- a/package/MDAnalysis/analysis/diffusionmap.py +++ b/package/MDAnalysis/analysis/diffusionmap.py @@ -227,11 +227,12 @@ 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. 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 c6170938d1e..f7227b620a8 100644 --- a/package/MDAnalysis/analysis/rms.py +++ b/package/MDAnalysis/analysis/rms.py @@ -567,7 +567,8 @@ def run(self, start=None, stop=None, step=None, progout=10, quiet=False): start : int (optional) starting frame, default None becomes 0. stop : int (optional) - stopping frame, default None becomes n_frames. + Frame index to stop analysis. Default: None becomes + n_frames. Iteration stops *before* this frame number. step : int (optional) step between frames, default None becomes 1. progout : int (optional) From 9e7d439bde51cdb42190db70ae6cd95b1531ec4e Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Tue, 26 Jul 2016 17:40:25 -0700 Subject: [PATCH 05/28] Check_slice index fix --- package/MDAnalysis/analysis/rms.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/package/MDAnalysis/analysis/rms.py b/package/MDAnalysis/analysis/rms.py index f7227b620a8..26be0ba8e13 100644 --- a/package/MDAnalysis/analysis/rms.py +++ b/package/MDAnalysis/analysis/rms.py @@ -584,9 +584,8 @@ def run(self, start=None, stop=None, step=None, progout=10, quiet=False): Calculating Corrected Sums of Squares and Products." Technometrics 4(3):419-420. """ - start, stop, step = self.atomgroup.universe.check_slice_indices(start, - stop, - step) + 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) From e9b65862e4b17d10ec037fe22130a4d0a6475d62 Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Tue, 26 Jul 2016 17:49:04 -0700 Subject: [PATCH 06/28] DCD fix --- package/MDAnalysis/coordinates/DCD.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/package/MDAnalysis/coordinates/DCD.py b/package/MDAnalysis/coordinates/DCD.py index 47acfd1a740..98bfdb6c8a7 100644 --- a/package/MDAnalysis/coordinates/DCD.py +++ b/package/MDAnalysis/coordinates/DCD.py @@ -523,9 +523,9 @@ def timeseries(self, asel, start=None, stop=None, step=None, 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=None, stop=None, skip=None): + def correl(self, timeseries, start=None, stop=None, step=None): """Populate a TimeseriesCollection object with timeseries computed from the trajectory :Arguments: @@ -543,7 +543,7 @@ def correl(self, timeseries, start=None, stop=None, skip=None): 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: From 04dc85d1d8be6c8ccd9875fb396453b25738b7c9 Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Tue, 26 Jul 2016 18:07:22 -0700 Subject: [PATCH 07/28] Changed RMSFarray to reflect correct data --- .../data/adk_oplsaa_CA_rmsf.npy | Bin 1844 -> 1792 bytes 1 file changed, 0 insertions(+), 0 deletions(-) diff --git a/testsuite/MDAnalysisTests/data/adk_oplsaa_CA_rmsf.npy b/testsuite/MDAnalysisTests/data/adk_oplsaa_CA_rmsf.npy index d28a8aa8fc182e78ca521ce41752051b7cb6db6d..34d754740a80aafbdcd41d2b04bb8d2685039a4d 100644 GIT binary patch literal 1792 zcmX}qc|6pK8wYSR-(zMBLx-bQk)y%7L$363+q7Ek3JohELUR3*a_vOg&6TrkD#Dg^ zCYqu{DbcpIEZecdIvUDuesZM4&hMY!^WXD%Ua#km&+DDmsl)uB;Ml+z zzLgr^W{(A5PmRB4e@t9VK&1cvn4sYQ<<0>SvB8o&HZ&kASh8!I7@O+p=&2>C{ePq{ z8Vba9{RF;_+BE-RIe~4J=J}?7FhN~@T(RWFM1ET|S!_n&bmds&#|i@XTCE0o1{5x2 z{q(u~O%l$bmVfbQ2#jR+oRw?HAR{KUpWRSHe|^$~s31Y}i2 z@9OI)Sovxh4Meg~|Ip0FJ%z@JlkQ<&ktDn}i1SpXNL1fS(6vh^L07x&?7l|e{j-)V z%d;dHyTq?^_EC5irM?jOgu=IeZE-_2BpQn!o(Wh>!Pd&FASsGMkB3!RjthmFntr})l6y^m_-tw(z^k-;z{3%T%Q8lX1Cy<4y3PYa!4i3KB`+(>_ zCk;-rRF>x>E-HOj8viDcMM9YJa!+k}82#3?ui<?LO#SGO%g{+26JVV)~scIxl} zcODOsjg=wut7R||o_Ql?BNxgF9IHP&WWbpZcv$c)52`Kesyz~Aa6Wx-ip%FB;$mG= z&~+X>hL_cBNRfj<#Lp8a&17)Dr}y>T_p-S1x8UH}QhD4ieBn7&B@d4p1%_I`93Bu~ z`j#D+#gWSDq{3q|=a9NNe0?i=}9A7S^R8z=fGVjS(qH}*;_Ux15Kgu>i9QYbiYpG&C0VO(og-p zWIr2?&Q%VVm1(S-(aG`Cpb+c0*2pc2#6YGp2EzzE^|GFOcbti@>qoY5NhZ_|namK) z3{d0abD9+lT-9|CXj+Rr*h0W3R?%{q2LaF-HmHUMVu;hZ86f{=i?YMc_qOW2c6QiPgQ2 z*Lu5>2$^}5mQzQ-*;Sx?aEd^ccmJ`3Z3LKW{<8lhPawg1TsNN2MEQ2p)S*lQ zmUKn2pDl%Q)z5oV>}g!ex;4G&FdLsq6Fx8R*~0w}J?-(r9>Yn-qV^JdoKkYC5_{WW ztTyyv?p1pfk9nBKx^0B%Vth>ff-Iu7ko$X=G(_6(dk9|+s`7qWC|}P)lUwwU)7xn{ z#BR8&T1mszRmdOdWI;-l?$VmeLFzyId^X=_Hn-WR5w)ad4s49G4roX}-JcaPEwrgH4@11l|LJ(M)^PY+a?!FEvNs zOdMEP~@N!?C_l12dW8dW zxF1Be`kP@`Dh|_EVc|`DcUR^g6xJNrDKz{RmeH_N z(&&C=N5QPa)G|;&A)sY$tp7U-=k=T|nlmU=Z~E@Oa3zUCmOwb5C~=Z^`{wWh0lR4P z^6+pHkK;yOmFSR&bbRB|BXMA)xnBC$Bnlg(Z-yt0l2Gg2K4vgTNc!TLXsZ)Ads_2@ z@jp!L6pfa)s!IInYBqPBBrrSKKJclKz@e$+GtJKlM7nBUQ)nWgSk!xNww6TlF>d0Z zBZ*#vK2@71B&IKAT;8fpA^J$I%8d{Tsj^E6=hY}Mn>%hh8j^70i-w#ejtRm09;F>5 zK7^$oAG=2&Hk&L;ZmAY>oL%0}7_ib{j=RzwlF~vup175Fo`?WzDHZe!Nx3#j6 z%YW{Y*3LqUU66iF42{mAh{I7TG~U+VRh-SFVX>-BKn##rKD_PO7ay6h3+(^UK@ga< L4eU7dhJ}9tQ9#c6 literal 1844 zcmXYxeK?f)8pr2(e$Tu;W)Pba5~<^~9oeEtHD9x`F)d|DayEI(a*7#En7kz0WP4$$ zbgb0n#AS&&X-PS>Nlii&Q5;e!)=S!6u-mIM*Y3Ihxc|Aoe|+!j{@$NdIy8FsPPxi1 zG$uCOZs+bD%E;i@*kBbK6dN8I6TM3ryE|0LDyTnIY;@S4Qw1&KsQsvrbRPw=l~J;L z>i-MElqz{Rt6==-?M^a|jM53Xwld;Jt@85n@|@PyA5F$lh8*;*;Bh#&qn-XTvzbPt z9eX&rnZnT0h_fbg3inS-)rx~01ZBT+A1R;^7?`xQVKoisO~!FQaR?OZpq3qm1lblV z7~^3YnVKtV4@E>2&nC-#h zxI|C zZ5JOAD$%0+Tt03U6$L!M!{F)qOV$279%^#R9J1T_@SG@#ZM)3_-1F@wEqJKbUp~w5 zB_B5r_fsbq2w)m{F+-ls$Lh3uO~UyC_`44u%C8G~_^fZsrRSG<*#Fgu7k2p!DvDbF z6?KAOSB865-~^4OE05m)dVs>ix~to~Z*ma!qf7k291dJR_?fwZgHpQWn}S6YhE;+V zr}rFGR@ASloI#=g=hA@S2ntEdj%*WI(3m&YotgTa#wAtvUH>u)ts|84iU$07) z=yI+~oY`!t%fS~Wdn1x(&~VtSDioK{_;vch%HSUeDi_Z-j(JK@%eN`nDPa)Sadb#+ z#z1gD;#j(bUioX!{og}<&#W)doCoJv`LbE7xG zg*L6C@;3tG3+rUh^cfT%&MmmqPjJvNBsS<7gHPopQ4Z-0Y=Sn7n(imCyK>4qyOQ9) ztF}rmj}s_c49gua>;5im;9PA4jSYSt8`kRf@-1K8g>()&wyP~peo5o=}+fwoDyH!=edxZQ=eUkcBysjrW*QuO8rPrJ5C5Wer{j>38=HYQd{ z`H4<&7&t_1Zmxn za`})Xl=+9r1#siBHw-3pIWk+su~;s|6}!Jz81zapQj`(8d{6=<>#H7oCPB>)7q z^D*#`W{z|=i(-pbz0@8Fl4Rr7;tDAsKcMkBD#fU_KYs6T`k2(lKOA~00#~DGZ8H>M zZ(7~_zCb5-nEHuX`FP57?kL^ ze5h7itg6=WIq$L)XBUI(LWS}9nOvl#G}%k*2)y5T1`IyuqQ2ebTJtd;cNaT$xl-P(q!cPx6+#;m;-3J{Wcux7=kJmcXM2LTFCuv7ifPO`DO+o7FMu+{?teqq>s Date: Tue, 26 Jul 2016 18:30:55 -0700 Subject: [PATCH 08/28] Updated changelog [skip ci] --- package/CHANGELOG | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/package/CHANGELOG b/package/CHANGELOG index 9a1bb52719d..bd287842752 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -34,7 +34,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 From cb5852e9cea9fa1affdd04297e00a06b51a66d7e Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Wed, 27 Jul 2016 15:56:26 -0700 Subject: [PATCH 09/28] Added test for slicing, slowly changing DCD source for timeseries slicing. --- package/MDAnalysis/coordinates/src/dcd.c | 49 ++++++++++++------- .../MDAnalysisTests/coordinates/test_dcd.py | 9 ++++ 2 files changed, 41 insertions(+), 17 deletions(-) diff --git a/package/MDAnalysis/coordinates/src/dcd.c b/package/MDAnalysis/coordinates/src/dcd.c index dea24b6884a..ebdbfaddcd1 100644 --- a/package/MDAnalysis/coordinates/src/dcd.c +++ b/package/MDAnalysis/coordinates/src/dcd.c @@ -435,7 +435,13 @@ __read_timeseries(PyObject *self, PyObject *args) int rc; int i, j, index; int n_atoms = 0, n_frames = 0; +<<<<<<< HEAD int start = 0, stop = -1, skip = 1, numskip = 0; +======= + /* Stop = -1 is incorrect and causes an error, look for a fix of this in line + 469 */ + int start = 0, stop = -1, skip = 1, numskip = 0, remaining_frames=0; +>>>>>>> Added test for slicing, slowly changing DCD source for timeseries slicing. dcdhandle *dcd = NULL; int *atomlist = NULL; npy_intp dimensions[3]; @@ -464,8 +470,10 @@ __read_timeseries(PyObject *self, PyObject *args) Py_DECREF(temp); // Assume that start and stop are valid - if (stop == -1) { stop = dcd->nsets; } - n_frames = (stop-start+1) / skip; + + if (stop == -1) { stop = dcd->nsets -1; } + n_frames = ((stop-start) / skip)+1; + //n_frames = dcd->nsets / skip; n_atoms = PyList_Size((PyObject*)atoms); if (n_atoms == 0) { @@ -497,20 +505,15 @@ __read_timeseries(PyObject *self, PyObject *args) // 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) { + } 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) { + } 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) { + } 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) { + } 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) { + } else if (strncasecmp(format, "cfa", 3) == 0) { dimensions[0] = 3; dimensions[1] = n_frames; dimensions[2] = n_atoms; } @@ -534,32 +537,43 @@ __read_timeseries(PyObject *self, PyObject *args) goto error; } - for (i=0;i 1) { // Figure out how many steps to skip - numskip = skip - (dcd->setsread % skip) - 1; + numskip = skip-1; + if(remaining_frames < numskip){ + numskip = 1; + } + printf("%d numskip, %d skip %d stop %d setsread %d n_frames %d i %d remaining_frames \n", numskip, skip, stop, dcd->setsread, n_frames, i, remaining_frames); 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; } + //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); dcd->first = 0; - dcd->setsread++; + dcd->setsread += numskip; + remaining_frames = stop - dcd->setsread++; if (rc < 0) { // 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]) */ @@ -616,6 +630,7 @@ __read_timeseries(PyObject *self, PyObject *args) return NULL; } + static PyObject * __read_next_frame(PyObject *self, PyObject *args) { diff --git a/testsuite/MDAnalysisTests/coordinates/test_dcd.py b/testsuite/MDAnalysisTests/coordinates/test_dcd.py index 0e901ee87fb..7868e8c00f6 100644 --- a/testsuite/MDAnalysisTests/coordinates/test_dcd.py +++ b/testsuite/MDAnalysisTests/coordinates/test_dcd.py @@ -137,6 +137,15 @@ 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 fail before issue #914 resolved + self.u = mda.Universe(PSF, DCD) + ts = self.u.trajectory.timeseries(self.u.atoms) + ts_skip = self.u.trajectory.timeseries(self.u.atoms, skip=10) + assert_array_almost_equal(ts[:,::10], 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, From 1b1d4f28b33684277f56e33712a213f0b9f5a4f5 Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Wed, 27 Jul 2016 15:57:04 -0700 Subject: [PATCH 10/28] DCD.py changes --- package/MDAnalysis/coordinates/DCD.py | 11 ++++++++--- package/MDAnalysis/coordinates/src/dcd.c | 10 +++------- 2 files changed, 11 insertions(+), 10 deletions(-) diff --git a/package/MDAnalysis/coordinates/DCD.py b/package/MDAnalysis/coordinates/DCD.py index 98bfdb6c8a7..84bab725130 100644 --- a/package/MDAnalysis/coordinates/DCD.py +++ b/package/MDAnalysis/coordinates/DCD.py @@ -498,7 +498,8 @@ def _read_frame(self, frame): ts.frame = frame return ts - def timeseries(self, asel, start=None, stop=None, step=None, format='afc'): + def timeseries(self, asel, start=None, stop=-1, step=None, skip=None, + format='afc'): """Return a subset of coordinate data for an AtomGroup :Arguments: @@ -514,7 +515,11 @@ def timeseries(self, asel, start=None, stop=None, step=None, format='afc'): where the shape is (frame, number of atoms, coordinates) """ - start, stop, skip = self.check_slice_indices(start, stop, step) + if skip is not None: + step = skip + start, stop, step = self.check_slice_indices(start, stop, step) + # for old cython _read_timeseries code + skip = 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']: @@ -523,7 +528,7 @@ def timeseries(self, asel, start=None, stop=None, step=None, 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, step, format) + return self._read_timeseries(atom_numbers, start, stop, skip, format) def correl(self, timeseries, start=None, stop=None, step=None): """Populate a TimeseriesCollection object with timeseries computed from the trajectory diff --git a/package/MDAnalysis/coordinates/src/dcd.c b/package/MDAnalysis/coordinates/src/dcd.c index ebdbfaddcd1..0e465cd7360 100644 --- a/package/MDAnalysis/coordinates/src/dcd.c +++ b/package/MDAnalysis/coordinates/src/dcd.c @@ -435,13 +435,10 @@ __read_timeseries(PyObject *self, PyObject *args) int rc; int i, j, index; int n_atoms = 0, n_frames = 0; -<<<<<<< HEAD - int start = 0, stop = -1, skip = 1, numskip = 0; -======= + /* Stop = -1 is incorrect and causes an error, look for a fix of this in line 469 */ int start = 0, stop = -1, skip = 1, numskip = 0, remaining_frames=0; ->>>>>>> Added test for slicing, slowly changing DCD source for timeseries slicing. dcdhandle *dcd = NULL; int *atomlist = NULL; npy_intp dimensions[3]; @@ -542,7 +539,7 @@ __read_timeseries(PyObject *self, PyObject *args) for (i=0;i 1) { // Figure out how many steps to skip - numskip = skip-1; + numskip = skip; if(remaining_frames < numskip){ numskip = 1; } @@ -560,7 +557,7 @@ __read_timeseries(PyObject *self, PyObject *args) dcd->reverse, dcd->charmm); dcd->first = 0; dcd->setsread += numskip; - remaining_frames = stop - dcd->setsread++; + remaining_frames = stop - dcd->setsread; if (rc < 0) { // return an exception PyErr_SetString(PyExc_IOError, "Error reading frame from DCD file"); @@ -568,7 +565,6 @@ __read_timeseries(PyObject *self, PyObject *args) } - } // Copy into Numeric array only those atoms we are interested in for (j=0;j Date: Wed, 27 Jul 2016 16:49:59 -0700 Subject: [PATCH 11/28] Everything working. Fixed logic error in code. Removed printing of debug statement. --- package/MDAnalysis/coordinates/src/dcd.c | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/package/MDAnalysis/coordinates/src/dcd.c b/package/MDAnalysis/coordinates/src/dcd.c index 0e465cd7360..4924a9c9b4f 100644 --- a/package/MDAnalysis/coordinates/src/dcd.c +++ b/package/MDAnalysis/coordinates/src/dcd.c @@ -535,15 +535,14 @@ __read_timeseries(PyObject *self, PyObject *args) } /*Should actually be number of frames, not true yet*/ - remaining_frames = stop; + remaining_frames = stop-start; for (i=0;i 1) { + if (skip > 1 && i>0) { // Figure out how many steps to skip - numskip = skip; + numskip = skip -1; if(remaining_frames < numskip){ numskip = 1; } - printf("%d numskip, %d skip %d stop %d setsread %d n_frames %d i %d remaining_frames \n", numskip, skip, stop, dcd->setsread, n_frames, i, remaining_frames); rc = skip_dcdstep(dcd->fd, dcd->natoms, dcd->nfixed, dcd->charmm, numskip); if (rc < 0) { // return an exception @@ -551,12 +550,13 @@ __read_timeseries(PyObject *self, PyObject *args) goto error; } } + 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); dcd->first = 0; - dcd->setsread += numskip; + dcd->setsread += 1; remaining_frames = stop - dcd->setsread; if (rc < 0) { // return an exception From 55910a0a279ca6e76998142355cf061345e14d8c Mon Sep 17 00:00:00 2001 From: jdetle Date: Wed, 27 Jul 2016 21:49:21 -0700 Subject: [PATCH 12/28] Change to ++ Small fixes to make code work better --- package/MDAnalysis/coordinates/DCD.py | 6 ++---- package/MDAnalysis/coordinates/src/dcd.c | 6 +++++- testsuite/MDAnalysisTests/coordinates/test_dcd.py | 2 +- 3 files changed, 8 insertions(+), 6 deletions(-) diff --git a/package/MDAnalysis/coordinates/DCD.py b/package/MDAnalysis/coordinates/DCD.py index 84bab725130..ad1819b9a8e 100644 --- a/package/MDAnalysis/coordinates/DCD.py +++ b/package/MDAnalysis/coordinates/DCD.py @@ -498,7 +498,7 @@ def _read_frame(self, frame): ts.frame = frame return ts - def timeseries(self, asel, start=None, stop=-1, step=None, skip=None, + def timeseries(self, asel, start=None, stop=None, step=None, skip=None, format='afc'): """Return a subset of coordinate data for an AtomGroup @@ -518,8 +518,6 @@ def timeseries(self, asel, start=None, stop=-1, step=None, skip=None, if skip is not None: step = skip start, stop, step = self.check_slice_indices(start, stop, step) - # for old cython _read_timeseries code - skip = 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']: @@ -528,7 +526,7 @@ def timeseries(self, asel, start=None, stop=-1, step=None, skip=None, # 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=None, stop=None, step=None): """Populate a TimeseriesCollection object with timeseries computed from the trajectory diff --git a/package/MDAnalysis/coordinates/src/dcd.c b/package/MDAnalysis/coordinates/src/dcd.c index 4924a9c9b4f..0498b548b04 100644 --- a/package/MDAnalysis/coordinates/src/dcd.c +++ b/package/MDAnalysis/coordinates/src/dcd.c @@ -467,8 +467,12 @@ __read_timeseries(PyObject *self, PyObject *args) Py_DECREF(temp); // Assume that start and stop are valid +<<<<<<< HEAD if (stop == -1) { stop = dcd->nsets -1; } +======= + if (stop == dcd->nsets) { stop = dcd->nsets -1; } +>>>>>>> Change to ++ n_frames = ((stop-start) / skip)+1; //n_frames = dcd->nsets / skip; @@ -556,7 +560,7 @@ __read_timeseries(PyObject *self, PyObject *args) unitcell, dcd->nfixed, dcd->first, dcd->freeind, dcd->fixedcoords, dcd->reverse, dcd->charmm); dcd->first = 0; - dcd->setsread += 1; + dcd->setsread++; remaining_frames = stop - dcd->setsread; if (rc < 0) { // return an exception diff --git a/testsuite/MDAnalysisTests/coordinates/test_dcd.py b/testsuite/MDAnalysisTests/coordinates/test_dcd.py index 7868e8c00f6..d915c3a739f 100644 --- a/testsuite/MDAnalysisTests/coordinates/test_dcd.py +++ b/testsuite/MDAnalysisTests/coordinates/test_dcd.py @@ -142,7 +142,7 @@ def test_timeseries_slicing(self): # should fail before issue #914 resolved self.u = mda.Universe(PSF, DCD) ts = self.u.trajectory.timeseries(self.u.atoms) - ts_skip = self.u.trajectory.timeseries(self.u.atoms, skip=10) + ts_skip = self.u.trajectory.timeseries(self.u.atoms, step=10) assert_array_almost_equal(ts[:,::10], ts_skip, 5) From 83c7e4f962122b684074f9fa770115a724a5922c Mon Sep 17 00:00:00 2001 From: jdetle Date: Wed, 27 Jul 2016 22:40:47 -0700 Subject: [PATCH 13/28] Elaborated on docs --- package/MDAnalysis/analysis/contacts.py | 3 ++- package/MDAnalysis/analysis/diffusionmap.py | 3 ++- package/MDAnalysis/analysis/rms.py | 3 ++- 3 files changed, 6 insertions(+), 3 deletions(-) diff --git a/package/MDAnalysis/analysis/contacts.py b/package/MDAnalysis/analysis/contacts.py index 6c03381ebcb..5a1f8cfb656 100644 --- a/package/MDAnalysis/analysis/contacts.py +++ b/package/MDAnalysis/analysis/contacts.py @@ -399,7 +399,8 @@ def __init__(self, u, selection, refgroup, method="hard_cut", radius=4.5, First frame of trajectory to analyse, Default: None becomes 0. stop : int, optional Frame index to stop analysis. Default: None becomes - n_frames. Iteration stops *before* this frame number. + 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: None becomes 1. diff --git a/package/MDAnalysis/analysis/diffusionmap.py b/package/MDAnalysis/analysis/diffusionmap.py index b42be7f044f..a6300212cc0 100644 --- a/package/MDAnalysis/analysis/diffusionmap.py +++ b/package/MDAnalysis/analysis/diffusionmap.py @@ -230,7 +230,8 @@ def __init__(self, u, select='all', metric=rmsd, cutoff=1E0-5, First frame of trajectory to analyse, Default: None becomes 0. stop : int, optional Frame index to stop analysis. Default: None becomes - n_frames. Iteration stops *before* this frame number. + 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: None becomes 1. """ diff --git a/package/MDAnalysis/analysis/rms.py b/package/MDAnalysis/analysis/rms.py index 26be0ba8e13..16ed9ae815b 100644 --- a/package/MDAnalysis/analysis/rms.py +++ b/package/MDAnalysis/analysis/rms.py @@ -568,7 +568,8 @@ def run(self, start=None, stop=None, step=None, progout=10, quiet=False): starting frame, default None becomes 0. stop : int (optional) Frame index to stop analysis. Default: None becomes - n_frames. Iteration stops *before* this frame number. + 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, default None becomes 1. progout : int (optional) From beffbeaa923c85598fc67a918d8dee8d7f2930b2 Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Thu, 28 Jul 2016 11:02:28 -0700 Subject: [PATCH 14/28] Added generator but it appears to not be working, please help. --- testsuite/MDAnalysisTests/coordinates/test_dcd.py | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/testsuite/MDAnalysisTests/coordinates/test_dcd.py b/testsuite/MDAnalysisTests/coordinates/test_dcd.py index d915c3a739f..0548933e281 100644 --- a/testsuite/MDAnalysisTests/coordinates/test_dcd.py +++ b/testsuite/MDAnalysisTests/coordinates/test_dcd.py @@ -140,10 +140,20 @@ def test_volume(self): def test_timeseries_slicing(self): # check that slicing behaves correctly # should fail before issue #914 resolved + x = [] + for i in range(10): + for j in range(i+1, 11): + for k in range(j, 11): + tuples = (i, j, k) + x.append(tuples) + for start, stop, step in x: + yield self._slice_generation_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, step=10) - assert_array_almost_equal(ts[:,::10], ts_skip, 5) + 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): From 8e24480f247759aeb05216ae5e6823d361e1e9f3 Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Thu, 28 Jul 2016 15:43:12 -0700 Subject: [PATCH 15/28] Deprecated skip for step --- package/MDAnalysis/coordinates/DCD.py | 18 +++++++++++++++++- 1 file changed, 17 insertions(+), 1 deletion(-) diff --git a/package/MDAnalysis/coordinates/DCD.py b/package/MDAnalysis/coordinates/DCD.py index ad1819b9a8e..92c4c711278 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 +from numpy.lib.utils import deprecate from ..core import flags from .. import units as mdaunits # use mdaunits instead of units to avoid a clash @@ -498,6 +499,7 @@ def _read_frame(self, frame): ts.frame = frame return ts + def timeseries(self, asel, start=None, stop=None, step=None, skip=None, format='afc'): """Return a subset of coordinate data for an AtomGroup @@ -514,9 +516,15 @@ def timeseries(self, asel, start=None, stop=None, step=None, skip=None, 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. """ if skip is not None: step = skip + + @deprecate(new_name="step", + message="This parameter will be removed in 0.17") start, stop, step = self.check_slice_indices(start, stop, step) if len(asel) == 0: raise NoDataError("Timeseries requires at least one atom to analyze") @@ -528,7 +536,7 @@ def timeseries(self, asel, start=None, stop=None, step=None, skip=None, # XXX needs to be implemented return self._read_timeseries(atom_numbers, start, stop, step, format) - def correl(self, timeseries, start=None, stop=None, step=None): + def correl(self, timeseries, start=None, stop=None, step=None, skip=None): """Populate a TimeseriesCollection object with timeseries computed from the trajectory :Arguments: @@ -537,7 +545,15 @@ def correl(self, timeseries, start=None, stop=None, step=None): *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. """ + if skip is not None: + step = skip + + @deprecate(new_name="step", + message="This parameter will be removed in 0.17") start, stop, step = self.check_slice_indices(start, stop, step) atomlist = timeseries._getAtomList() format = timeseries._getFormat() From d687cc7492b611c51517569bad1d74b32b95c77b Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Fri, 29 Jul 2016 08:00:41 -0700 Subject: [PATCH 16/28] Changed to warnings.warn() because decorator wasnt working --- package/MDAnalysis/coordinates/DCD.py | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/package/MDAnalysis/coordinates/DCD.py b/package/MDAnalysis/coordinates/DCD.py index 92c4c711278..fce5df80de0 100644 --- a/package/MDAnalysis/coordinates/DCD.py +++ b/package/MDAnalysis/coordinates/DCD.py @@ -73,7 +73,7 @@ import numpy as np import struct import types -from numpy.lib.utils import deprecate +import warnings from ..core import flags from .. import units as mdaunits # use mdaunits instead of units to avoid a clash @@ -86,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 @@ -499,7 +500,6 @@ def _read_frame(self, frame): ts.frame = frame return ts - def timeseries(self, asel, start=None, stop=None, step=None, skip=None, format='afc'): """Return a subset of coordinate data for an AtomGroup @@ -522,9 +522,10 @@ def timeseries(self, asel, start=None, stop=None, step=None, skip=None, """ 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) - @deprecate(new_name="step", - message="This parameter will be removed in 0.17") start, stop, step = self.check_slice_indices(start, stop, step) if len(asel) == 0: raise NoDataError("Timeseries requires at least one atom to analyze") @@ -551,10 +552,10 @@ def correl(self, timeseries, start=None, stop=None, step=None, skip=None): """ 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) - @deprecate(new_name="step", - message="This parameter will be removed in 0.17") - start, stop, step = self.check_slice_indices(start, stop, step) atomlist = timeseries._getAtomList() format = timeseries._getFormat() lowerb, upperb = timeseries._getBounds() From 6893abdb220980432cb51c136036fc526d3376ea Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Fri, 29 Jul 2016 09:01:10 -0700 Subject: [PATCH 17/28] Fixed silly correl bug Simple arg changes --- package/MDAnalysis/coordinates/DCD.py | 1 + package/MDAnalysis/coordinates/src/dcd.c | 14 +++++++++----- 2 files changed, 10 insertions(+), 5 deletions(-) diff --git a/package/MDAnalysis/coordinates/DCD.py b/package/MDAnalysis/coordinates/DCD.py index fce5df80de0..e121a098295 100644 --- a/package/MDAnalysis/coordinates/DCD.py +++ b/package/MDAnalysis/coordinates/DCD.py @@ -556,6 +556,7 @@ def correl(self, timeseries, start=None, stop=None, step=None, skip=None): "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() diff --git a/package/MDAnalysis/coordinates/src/dcd.c b/package/MDAnalysis/coordinates/src/dcd.c index 0498b548b04..827b9411dc7 100644 --- a/package/MDAnalysis/coordinates/src/dcd.c +++ b/package/MDAnalysis/coordinates/src/dcd.c @@ -438,7 +438,7 @@ __read_timeseries(PyObject *self, PyObject *args) /* Stop = -1 is incorrect and causes an error, look for a fix of this in line 469 */ - int start = 0, stop = -1, skip = 1, numskip = 0, remaining_frames=0; + int start = 0, stop = -1, step = 1, numskip = 0, remaining_frames=0; dcdhandle *dcd = NULL; int *atomlist = NULL; npy_intp dimensions[3]; @@ -448,12 +448,12 @@ __read_timeseries(PyObject *self, PyObject *args) 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) ) + 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, &skip, &format) ) + if( !PyArg_ParseTuple(args, "O!|iiis", &PyList_Type, &atoms, &start, &stop, &step, &format) ) return NULL; } @@ -472,9 +472,13 @@ __read_timeseries(PyObject *self, PyObject *args) if (stop == -1) { stop = dcd->nsets -1; } ======= if (stop == dcd->nsets) { stop = dcd->nsets -1; } +<<<<<<< HEAD >>>>>>> Change to ++ n_frames = ((stop-start) / skip)+1; +======= + n_frames = ((stop-start) / step)+1; +>>>>>>> Fixed silly correl bug //n_frames = dcd->nsets / skip; n_atoms = PyList_Size((PyObject*)atoms); if (n_atoms == 0) { @@ -541,9 +545,9 @@ __read_timeseries(PyObject *self, PyObject *args) /*Should actually be number of frames, not true yet*/ remaining_frames = stop-start; for (i=0;i 1 && i>0) { + if (step > 1 && i>0) { // Figure out how many steps to skip - numskip = skip -1; + numskip = step -1; if(remaining_frames < numskip){ numskip = 1; } From 6b275882712faaea80e3ab14b6766281f4c38862 Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Fri, 29 Jul 2016 16:13:14 -0700 Subject: [PATCH 18/28] Fixed issue caused by integer division, when stop-start / nframes has no remainder it should remain 1. --- package/MDAnalysis/coordinates/src/dcd.c | 13 ++----------- .../MDAnalysisTests/coordinates/test_dcd.py | 19 +++++++++++++++++-- 2 files changed, 19 insertions(+), 13 deletions(-) diff --git a/package/MDAnalysis/coordinates/src/dcd.c b/package/MDAnalysis/coordinates/src/dcd.c index 827b9411dc7..772bd2fb513 100644 --- a/package/MDAnalysis/coordinates/src/dcd.c +++ b/package/MDAnalysis/coordinates/src/dcd.c @@ -466,19 +466,10 @@ __read_timeseries(PyObject *self, PyObject *args) dcd = (dcdhandle*)PyCObject_AsVoidPtr(temp); Py_DECREF(temp); - // Assume that start and stop are valid -<<<<<<< HEAD - if (stop == -1) { stop = dcd->nsets -1; } -======= if (stop == dcd->nsets) { stop = dcd->nsets -1; } -<<<<<<< HEAD ->>>>>>> Change to ++ - n_frames = ((stop-start) / skip)+1; - -======= - n_frames = ((stop-start) / step)+1; ->>>>>>> Fixed silly correl bug + n_frames = ((stop-start) / step)+ 1; + if ((stop-start) == step) { n_frames = 1;} //n_frames = dcd->nsets / skip; n_atoms = PyList_Size((PyObject*)atoms); if (n_atoms == 0) { diff --git a/testsuite/MDAnalysisTests/coordinates/test_dcd.py b/testsuite/MDAnalysisTests/coordinates/test_dcd.py index 0548933e281..b336e5be3c1 100644 --- a/testsuite/MDAnalysisTests/coordinates/test_dcd.py +++ b/testsuite/MDAnalysisTests/coordinates/test_dcd.py @@ -55,7 +55,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") @@ -154,7 +164,12 @@ def _slice_generation_test(self, start, stop, step): 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_slice_generation_test(self): + self.u = mda.Universe(PSF, DCD) + ts = self.u.trajectory.timeseries(self.u.atoms) + ts_skip = self.u.trajectory.timeseries(self.u.atoms, 0, 4, 4) + assert_array_almost_equal(ts[:, 0:4:4,:], ts_skip, 5) def test_DCDReader_set_dt(dt=100., frame=3): u = mda.Universe(PSF, DCD, dt=dt) From 32fce067fefbef6b569d7af2e78063472f98a8e8 Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Fri, 29 Jul 2016 21:08:02 -0700 Subject: [PATCH 19/28] Reduced number of tests by a signifcant margin, removed test of test --- testsuite/MDAnalysisTests/coordinates/test_dcd.py | 11 +++-------- 1 file changed, 3 insertions(+), 8 deletions(-) diff --git a/testsuite/MDAnalysisTests/coordinates/test_dcd.py b/testsuite/MDAnalysisTests/coordinates/test_dcd.py index b336e5be3c1..c62ecb91776 100644 --- a/testsuite/MDAnalysisTests/coordinates/test_dcd.py +++ b/testsuite/MDAnalysisTests/coordinates/test_dcd.py @@ -151,9 +151,9 @@ def test_timeseries_slicing(self): # check that slicing behaves correctly # should fail before issue #914 resolved x = [] - for i in range(10): - for j in range(i+1, 11): - for k in range(j, 11): + for i in range(5): + for j in range(i+1, 6): + for k in range(j, 6): tuples = (i, j, k) x.append(tuples) for start, stop, step in x: @@ -165,11 +165,6 @@ def _slice_generation_test(self, start, stop, step): 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_slice_generation_test(self): - self.u = mda.Universe(PSF, DCD) - ts = self.u.trajectory.timeseries(self.u.atoms) - ts_skip = self.u.trajectory.timeseries(self.u.atoms, 0, 4, 4) - assert_array_almost_equal(ts[:, 0:4:4,:], ts_skip, 5) def test_DCDReader_set_dt(dt=100., frame=3): u = mda.Universe(PSF, DCD, dt=dt) From 68ae89786ab022d73ca3b0575a4bc14b1d80c6e0 Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Sat, 30 Jul 2016 13:13:06 -0700 Subject: [PATCH 20/28] Fixed incomplete logic Changed tuples Eliminated _testDCD --- package/MDAnalysis/coordinates/src/dcd.c | 2 +- .../MDAnalysisTests/coordinates/test_dcd.py | 24 ++++--------------- 2 files changed, 5 insertions(+), 21 deletions(-) diff --git a/package/MDAnalysis/coordinates/src/dcd.c b/package/MDAnalysis/coordinates/src/dcd.c index 772bd2fb513..88e09293530 100644 --- a/package/MDAnalysis/coordinates/src/dcd.c +++ b/package/MDAnalysis/coordinates/src/dcd.c @@ -469,7 +469,7 @@ __read_timeseries(PyObject *self, PyObject *args) if (stop == dcd->nsets) { stop = dcd->nsets -1; } n_frames = ((stop-start) / step)+ 1; - if ((stop-start) == step) { n_frames = 1;} + if ((stop-start) % step == 0) { n_frames--;} //n_frames = dcd->nsets / skip; n_atoms = PyList_Size((PyObject*)atoms); if (n_atoms == 0) { diff --git a/testsuite/MDAnalysisTests/coordinates/test_dcd.py b/testsuite/MDAnalysisTests/coordinates/test_dcd.py index c62ecb91776..0000d447058 100644 --- a/testsuite/MDAnalysisTests/coordinates/test_dcd.py +++ b/testsuite/MDAnalysisTests/coordinates/test_dcd.py @@ -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 @@ -150,12 +138,8 @@ def test_volume(self): def test_timeseries_slicing(self): # check that slicing behaves correctly # should fail before issue #914 resolved - x = [] - for i in range(5): - for j in range(i+1, 6): - for k in range(j, 6): - tuples = (i, j, k) - x.append(tuples) + 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)] for start, stop, step in x: yield self._slice_generation_test, start, stop, step @@ -164,7 +148,7 @@ def _slice_generation_test(self, start, stop, step): 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) @@ -445,7 +429,7 @@ 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() From 6f9e57d874524887eee4e9fe96a3d4b6731b924e Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Sat, 30 Jul 2016 18:47:03 -0700 Subject: [PATCH 21/28] Changes to DCD correl functions --- .../MDAnalysis/coordinates/dcdtimeseries.pyx | 25 +++++++++++++------ package/MDAnalysis/coordinates/src/dcd.c | 4 +-- package/MDAnalysis/core/Timeseries.py | 25 +++++++++++++++---- .../MDAnalysisTests/coordinates/test_dcd.py | 9 ++++--- 4 files changed, 45 insertions(+), 18 deletions(-) diff --git a/package/MDAnalysis/coordinates/dcdtimeseries.pyx b/package/MDAnalysis/coordinates/dcdtimeseries.pyx index c60e20589c8..db1c94fd9a4 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,10 @@ 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 + + if stop == dcd.nsets: stop -= 1; + n_frames = ((stop-start) / step) + 1; + if (stop-start) % step == 0: n_frames -= 1; cdef int numdata numdata = len(format) if numdata==0: @@ -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 + print n_frames 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): + 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 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 88e09293530..714e5f73cad 100644 --- a/package/MDAnalysis/coordinates/src/dcd.c +++ b/package/MDAnalysis/coordinates/src/dcd.c @@ -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); diff --git a/package/MDAnalysis/core/Timeseries.py b/package/MDAnalysis/core/Timeseries.py index 45a3776f341..317b62a8cfd 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,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" diff --git a/testsuite/MDAnalysisTests/coordinates/test_dcd.py b/testsuite/MDAnalysisTests/coordinates/test_dcd.py index 0000d447058..0f058456697 100644 --- a/testsuite/MDAnalysisTests/coordinates/test_dcd.py +++ b/testsuite/MDAnalysisTests/coordinates/test_dcd.py @@ -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 @@ -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)") From e476ee5c2d913aabc936c7330f391f666f3fe6c1 Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Sat, 30 Jul 2016 19:48:25 -0700 Subject: [PATCH 22/28] Changed to logic to take advantage of new default, simplify code. --- package/MDAnalysis/coordinates/dcdtimeseries.pyx | 6 ++---- package/MDAnalysis/coordinates/src/dcd.c | 6 ++---- 2 files changed, 4 insertions(+), 8 deletions(-) diff --git a/package/MDAnalysis/coordinates/dcdtimeseries.pyx b/package/MDAnalysis/coordinates/dcdtimeseries.pyx index db1c94fd9a4..d7c2ab4cb48 100644 --- a/package/MDAnalysis/coordinates/dcdtimeseries.pyx +++ b/package/MDAnalysis/coordinates/dcdtimeseries.pyx @@ -70,9 +70,8 @@ def __read_timecorrel(object self, object atoms, object atomcounts, object forma dcd = PyCObject_AsVoidPtr(self._dcd_C_ptr) cdef int n_frames - if stop == dcd.nsets: stop -= 1; - n_frames = ((stop-start) / step) + 1; - if (stop-start) % step == 0: n_frames -= 1; + n_frames = ((stop-start) / step); + if (stop-start) % step > 0 : n_frames += 1; cdef int numdata numdata = len(format) if numdata==0: @@ -113,7 +112,6 @@ def __read_timecorrel(object self, object atoms, object atomcounts, object forma cdef float unitcell[6] cdef int remaining_frames = stop-start numskip = 0 - print n_frames for i from 0 <= i < n_frames: if (step > 1 and i > 0): # Check if we have fixed atoms diff --git a/package/MDAnalysis/coordinates/src/dcd.c b/package/MDAnalysis/coordinates/src/dcd.c index 714e5f73cad..9f3d6fe222f 100644 --- a/package/MDAnalysis/coordinates/src/dcd.c +++ b/package/MDAnalysis/coordinates/src/dcd.c @@ -466,10 +466,8 @@ __read_timeseries(PyObject *self, PyObject *args) dcd = (dcdhandle*)PyCObject_AsVoidPtr(temp); Py_DECREF(temp); - - if (stop == dcd->nsets) { stop = dcd->nsets - 1; } - n_frames = ((stop-start) / step) + 1; - if ((stop-start) % step == 0) { n_frames--;} + n_frames = ((stop-start) / step); + if ((stop-start) % step > 0) { n_frames++; } //n_frames = dcd->nsets / skip; n_atoms = PyList_Size((PyObject*)atoms); if (n_atoms == 0) { From d036859610ca8253a179078e3c0e725e33f757e9 Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Sat, 30 Jul 2016 20:55:44 -0700 Subject: [PATCH 23/28] I made a silly mistake in Timeseries --- package/MDAnalysis/core/Timeseries.py | 4 ++-- testsuite/MDAnalysisTests/coordinates/test_dcd.py | 10 ++++++++++ 2 files changed, 12 insertions(+), 2 deletions(-) diff --git a/package/MDAnalysis/core/Timeseries.py b/package/MDAnalysis/core/Timeseries.py index 317b62a8cfd..f73ab6a701b 100644 --- a/package/MDAnalysis/core/Timeseries.py +++ b/package/MDAnalysis/core/Timeseries.py @@ -128,8 +128,8 @@ def compute(self, trj, start=None, stop=None, skip=None, step=None): "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) + 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 0f058456697..ee928803d35 100644 --- a/testsuite/MDAnalysisTests/coordinates/test_dcd.py +++ b/testsuite/MDAnalysisTests/coordinates/test_dcd.py @@ -437,7 +437,9 @@ def setUp(self): 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]]) @@ -454,11 +456,16 @@ 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 + del self.collection_slicing del self.universe del self.dcd del self.ts @@ -526,6 +533,9 @@ 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") # notes: def compute_correl_references(): From 0c268ed06c1034336a0c2dfc9ab07e0d3894bc62 Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Sat, 30 Jul 2016 21:04:36 -0700 Subject: [PATCH 24/28] More testing --- testsuite/MDAnalysisTests/coordinates/test_dcd.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/testsuite/MDAnalysisTests/coordinates/test_dcd.py b/testsuite/MDAnalysisTests/coordinates/test_dcd.py index ee928803d35..da17411aeb4 100644 --- a/testsuite/MDAnalysisTests/coordinates/test_dcd.py +++ b/testsuite/MDAnalysisTests/coordinates/test_dcd.py @@ -536,6 +536,8 @@ def test_clear(self): 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(): From d3d4dadac56fef66101e461140d0856143198081 Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Sun, 31 Jul 2016 11:50:03 -0700 Subject: [PATCH 25/28] Added comments. --- package/MDAnalysis/coordinates/dcdtimeseries.pyx | 10 ++++++++-- package/MDAnalysis/coordinates/src/dcd.c | 10 ++++++++-- 2 files changed, 16 insertions(+), 4 deletions(-) diff --git a/package/MDAnalysis/coordinates/dcdtimeseries.pyx b/package/MDAnalysis/coordinates/dcdtimeseries.pyx index d7c2ab4cb48..aaec810c54a 100644 --- a/package/MDAnalysis/coordinates/dcdtimeseries.pyx +++ b/package/MDAnalysis/coordinates/dcdtimeseries.pyx @@ -115,15 +115,21 @@ def __read_timecorrel(object self, object atoms, object atomcounts, object forma for i from 0 <= i < n_frames: if (step > 1 and 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): 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 += 1 diff --git a/package/MDAnalysis/coordinates/src/dcd.c b/package/MDAnalysis/coordinates/src/dcd.c index 9f3d6fe222f..6b898e3c897 100644 --- a/package/MDAnalysis/coordinates/src/dcd.c +++ b/package/MDAnalysis/coordinates/src/dcd.c @@ -531,12 +531,17 @@ __read_timeseries(PyObject *self, PyObject *args) goto error; } - /*Should actually be number of frames, not true yet*/ remaining_frames = stop-start; for (i=0;i 1 && i>0) { - // Figure out how many steps to skip + /* 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; } @@ -547,6 +552,7 @@ __read_timeseries(PyObject *self, PyObject *args) 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, From 6833e3a772be12069a9953f667fd395b9893c55a Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Sun, 31 Jul 2016 14:53:17 -0700 Subject: [PATCH 26/28] Undo comment removal [skip ci] --- package/MDAnalysis/coordinates/dcdtimeseries.pyx | 1 + package/MDAnalysis/coordinates/src/dcd.c | 2 ++ 2 files changed, 3 insertions(+) diff --git a/package/MDAnalysis/coordinates/dcdtimeseries.pyx b/package/MDAnalysis/coordinates/dcdtimeseries.pyx index aaec810c54a..1f1d4f585e8 100644 --- a/package/MDAnalysis/coordinates/dcdtimeseries.pyx +++ b/package/MDAnalysis/coordinates/dcdtimeseries.pyx @@ -115,6 +115,7 @@ def __read_timecorrel(object self, object atoms, object atomcounts, object forma for i from 0 <= i < n_frames: if (step > 1 and 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 diff --git a/package/MDAnalysis/coordinates/src/dcd.c b/package/MDAnalysis/coordinates/src/dcd.c index 6b898e3c897..0edc7abd29d 100644 --- a/package/MDAnalysis/coordinates/src/dcd.c +++ b/package/MDAnalysis/coordinates/src/dcd.c @@ -534,6 +534,8 @@ __read_timeseries(PyObject *self, PyObject *args) 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; From 76ec586cb22a520dcf21fa5ac390948deafece6a Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Mon, 1 Aug 2016 10:04:40 -0700 Subject: [PATCH 27/28] Added tests that cause failures --- testsuite/MDAnalysisTests/coordinates/test_dcd.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/testsuite/MDAnalysisTests/coordinates/test_dcd.py b/testsuite/MDAnalysisTests/coordinates/test_dcd.py index da17411aeb4..c7ebb4e523e 100644 --- a/testsuite/MDAnalysisTests/coordinates/test_dcd.py +++ b/testsuite/MDAnalysisTests/coordinates/test_dcd.py @@ -139,7 +139,7 @@ def test_timeseries_slicing(self): # check that slicing behaves correctly # should fail 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)] + (0, 5, 5), (3, 5, 1), (5, 0, -2), (5, 0, -1), (None, None, None)] for start, stop, step in x: yield self._slice_generation_test, start, stop, step From 0f62787ec3a667ff4f77807df03e9ec5c403865c Mon Sep 17 00:00:00 2001 From: John Detlefs Date: Mon, 1 Aug 2016 13:13:46 -0700 Subject: [PATCH 28/28] Added known failures --- package/MDAnalysis/coordinates/src/dcd.c | 2 +- .../MDAnalysisTests/coordinates/test_dcd.py | 18 +++++++++++++++--- 2 files changed, 16 insertions(+), 4 deletions(-) diff --git a/package/MDAnalysis/coordinates/src/dcd.c b/package/MDAnalysis/coordinates/src/dcd.c index 0edc7abd29d..d0d3956dbe9 100644 --- a/package/MDAnalysis/coordinates/src/dcd.c +++ b/package/MDAnalysis/coordinates/src/dcd.c @@ -535,7 +535,7 @@ __read_timeseries(PyObject *self, PyObject *args) for (i=0;i 1 && i>0) { // Check if we have fixed atoms -- // XXX not done + // 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; diff --git a/testsuite/MDAnalysisTests/coordinates/test_dcd.py b/testsuite/MDAnalysisTests/coordinates/test_dcd.py index c7ebb4e523e..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(): @@ -137,18 +137,30 @@ def test_volume(self): def test_timeseries_slicing(self): # check that slicing behaves correctly - # should fail before issue #914 resolved + # 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), (5, 0, -2), (5, 0, -1), (None, None, None)] + (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)