DensityAnalysis and RMSF start/stop/step Issue Fix (#111)

* Fix #108 
* Corrected start/stop/step issue in DensityAnalysis and RMSF
* Changes:
  - Added self.n_frames to parallel.py
  - Updated DensityAnalysis conclude() to read self.n_frames instead of trajectory length
  - Updated RMSF to read self.n_frames in _conclude()
  - Added test_n_frames() to both DensityAnalysis and RMSF tests
  - Added test_accuracy() to test_rmsf.py to verify subselecting trajectories doesn't affect accuracy
  - Added new datafiles to test_density.py to use for testing timestep sub-selections
  - Consolidated two tests in test_density.py into one test that verifies both the array values and the sum of their values
* Fixed commented out metadata in conclude()
* PEP8 fix

Author: nawtrey

pmda/density.py

@@ -251,7 +251,6 @@ def __init__(self, atomgroup, delta=1.0, atomselection=None,
         self._ydim = ydim
         self._zdim = zdim
         self._trajectory = u.trajectory
-        self._n_frames = u.trajectory.n_frames
         if updating and atomselection is None:
             raise ValueError("updating=True requires a atomselection string")
         elif not updating and atomselection is not None:
@@ -300,12 +299,12 @@ def _single_frame(self, ts, atomgroups):
 
     def _conclude(self):
         self._grid = self._results[:].sum(axis=0)
-        self._grid /= float(self._n_frames)
+        self._grid /= float(self.n_frames)
         metadata = self._metadata if self._metadata is not None else {}
         metadata['psf'] = self._atomgroup.universe.filename
         metadata['dcd'] = self._trajectory.filename
         metadata['atomselection'] = self._atomselection
-        metadata['n_frames'] = self._n_frames
+        metadata['n_frames'] = self.n_frames
         metadata['totaltime'] = self._atomgroup.universe.trajectory.totaltime
         metadata['dt'] = self._trajectory.dt
         metadata['time_unit'] = mda.core.flags['time_unit']

pmda/parallel.py

@@ -352,6 +352,8 @@ def run(self,
                                                                  stop, step)
         n_frames = len(range(start, stop, step))
 
+        self.n_frames = n_frames
+
         if n_frames == 0:
             warnings.warn("run() analyses no frames: check start/stop/step")
         if n_frames < n_blocks:

pmda/rms/rmsf.py

@@ -193,10 +193,9 @@ def _conclude(self):
         # serial case
         if n_blocks == 1:
             # get length of trajectory slice
-            k = len(self._blocks[0])
             self.mean = self._results[0, 0]
             self.sumsquares = self._results[0, 1]
-            self.rmsf = np.sqrt(self.sumsquares.sum(axis=1) / k)
+            self.rmsf = np.sqrt(self.sumsquares.sum(axis=1) / self.n_frames)
         # parallel case
         else:
             mean = self._results[:, 0]
@@ -207,10 +206,9 @@ def _conclude(self):
                 vals.append((len(self._blocks[i]), mean[i], sos[i]))
             # combine block results using fold method
             results = fold_second_order_moments(vals)
-            self.totalts = results[0]
             self.mean = results[1]
             self.sumsquares = results[2]
-            self.rmsf = np.sqrt(self.sumsquares.sum(axis=1) / self.totalts)
+            self.rmsf = np.sqrt(self.sumsquares.sum(axis=1) / self.n_frames)
             self._negative_rmsf(self.rmsf)
 
     @staticmethod

pmda/test/test_density.py

@@ -1,76 +1,77 @@
 import pytest
 import MDAnalysis as mda
 import numpy as np
-import pmda.density
+from pmda.density import DensityAnalysis
 from MDAnalysis.analysis.density import density_from_Universe
 from numpy.testing import assert_almost_equal
 from MDAnalysisTests.datafiles import GRO_MEMPROT, XTC_MEMPROT
+from MDAnalysisTests.datafiles import PSF_TRICLINIC, DCD_TRICLINIC
 
 
 @pytest.fixture(scope='module')
-def u():
+def u1():
     return mda.Universe(GRO_MEMPROT, XTC_MEMPROT)
 
 
-def test_density_sum(u):
-    pdensity = pmda.density.DensityAnalysis(u.atoms, atomselection='name OD1',
-                                            updating=True)
-    core_number = 2
-    pdensity.run(n_blocks=core_number, n_jobs=core_number)
-    dens = mda.analysis.density.density_from_Universe(u,
-                                                      atomselection='name OD1',
-                                                      update_selection=True)
-    assert np.sum(dens.grid) == np.sum(pdensity.density.grid)
+@pytest.fixture(scope='module')
+def u2():
+    return mda.Universe(PSF_TRICLINIC, DCD_TRICLINIC)
 
 
-@pytest.mark.parametrize("n_blocks", [1, 2, 3, 4])
-def test_density_grid(u, n_blocks):
-    pdensity = pmda.density.DensityAnalysis(u.atoms, atomselection='name OD1',
-                                            updating=True)
-    core_number = n_blocks
-    pdensity.run(n_blocks=core_number, n_jobs=core_number)
-    dens = mda.analysis.density.density_from_Universe(u,
-                                                      atomselection='name OD1',
-                                                      update_selection=True)
-    assert_almost_equal(dens.grid, pdensity.density.grid)
+@pytest.mark.parametrize("n_blocks", [1, 2])
+@pytest.mark.parametrize("stop", [4, 5])
+@pytest.mark.parametrize("step", [1, 2])
+def test_density_values(u1, n_blocks, stop, step):
+    parallel = DensityAnalysis(u1.atoms, atomselection='name OD1',
+                               updating=True)
+    parallel.run(n_blocks=n_blocks, n_jobs=n_blocks, start=0, stop=stop,
+                 step=step)
+    serial = density_from_Universe(u1, atomselection='name OD1',
+                                   update_selection=True, start=0, stop=stop,
+                                   step=step)
+    assert np.sum(serial.grid) == np.sum(parallel.density.grid)
+    assert_almost_equal(serial.grid, parallel.density.grid, decimal=17)
 
 
-def test_updating(u):
+def test_updating(u1):
     with pytest.raises(ValueError):
-        pdensity = pmda.density.DensityAnalysis(u.atoms,
-                                                updating=True)
+        pdensity = DensityAnalysis(u1.atoms, updating=True)
 
 
-def test_atomselection(u):
+def test_atomselection(u1):
     with pytest.raises(ValueError):
-        pdensity = pmda.density.DensityAnalysis(u.atoms,
-                                                atomselection='name OD1')
+        pdensity = DensityAnalysis(u1.atoms, atomselection='name OD1')
 
 
-def test_gridcenter(u):
-    aselect = 'name OD1'
+def test_gridcenter(u1):
     gridcenter = np.array([10, 10, 10])
     xdim = 190
     ydim = 200
     zdim = 210
-    dens = mda.analysis.density.density_from_Universe(u,
-                                                      atomselection=aselect,
-                                                      update_selection=True,
-                                                      gridcenter=gridcenter,
-                                                      xdim=xdim,
-                                                      ydim=ydim,
-                                                      zdim=zdim)
-    pdens = pmda.density.DensityAnalysis(u.atoms,
-                                         atomselection=aselect,
-                                         updating=True,
-                                         gridcenter=gridcenter,
-                                         xdim=xdim,
-                                         ydim=ydim,
-                                         zdim=zdim)
+    serial = density_from_Universe(u1, atomselection='name OD1',
+                                   update_selection=True,
+                                   gridcenter=gridcenter,
+                                   xdim=xdim,
+                                   ydim=ydim,
+                                   zdim=zdim)
+    parallel = DensityAnalysis(u1.atoms, atomselection='name OD1',
+                               updating=True,
+                               gridcenter=gridcenter,
+                               xdim=xdim,
+                               ydim=ydim,
+                               zdim=zdim)
     core_number = 4
-    pdens.run(n_blocks=core_number, n_jobs=core_number)
-    assert_almost_equal(dens.grid, pdens.density.grid)
-    assert_almost_equal(pdens._gridcenter, gridcenter)
-    assert len(pdens.density.edges[0]) == xdim + 1
-    assert len(pdens.density.edges[1]) == ydim + 1
-    assert len(pdens.density.edges[2]) == zdim + 1
+    parallel.run(n_blocks=core_number, n_jobs=core_number)
+    assert_almost_equal(serial.grid, parallel.density.grid)
+    assert_almost_equal(parallel._gridcenter, gridcenter)
+    assert len(parallel.density.edges[0]) == xdim + 1
+    assert len(parallel.density.edges[1]) == ydim + 1
+    assert len(parallel.density.edges[2]) == zdim + 1
+
+
+@pytest.mark.parametrize("start", [0, 1])
+@pytest.mark.parametrize("stop", [5, 7, 10])
+@pytest.mark.parametrize("step", [1, 2, 3])
+def test_n_frames(u2, start, stop, step):
+    pdensity = DensityAnalysis(u2.atoms).run(start, stop, step)
+    assert pdensity.n_frames == len(range(start, stop, step))

pmda/test/test_rmsf.py

@@ -23,8 +23,22 @@ def u():
     return mda.Universe(PSF, DCD)
 
 
+@pytest.mark.parametrize('n_blocks', (2, 3, 4))
+@pytest.mark.parametrize('decimal', (7, 10, 11))
+@pytest.mark.parametrize("step", [1, 2, 3, 5])
+def test_rmsf_accuracy(u, n_blocks, decimal, step):
+    PMDA_vals = RMSF(u.atoms).run(step=step,
+                                  n_blocks=n_blocks,
+                                  n_jobs=n_blocks)
+    MDA_vals = mda.analysis.rms.RMSF(u.atoms).run(step=step)
+    assert_almost_equal(MDA_vals.mean, PMDA_vals.mean, decimal=decimal)
+    assert_almost_equal(MDA_vals.sumsquares, PMDA_vals.sumsquares,
+                        decimal=decimal)
+    assert_almost_equal(MDA_vals.rmsf, PMDA_vals.rmsf, decimal=decimal)
+
+
 @pytest.mark.parametrize('n_blocks', (1, 2, 3, 4, 5, 10))
-@pytest.mark.parametrize('n_frames', (10, 50, 100))
+@pytest.mark.parametrize('n_frames', (10, 45, 90))
 def test_RMSF_values(u, n_blocks, n_frames):
     PMDA_vals = RMSF(u.atoms).run(stop=n_frames,
                                   n_blocks=n_blocks,
@@ -49,3 +63,11 @@ def test_RMSF_n_jobs(u, n_blocks):
 def test_negative_RMSF_raises_ValueError():
     with pytest.raises(ValueError):
         RMSF._negative_rmsf(np.array([-1, -1, -1]))
+
+
+@pytest.mark.parametrize("stop", [10, 33, 45, 90])
+@pytest.mark.parametrize("step", [1, 2, 3, 5])
+def test_n_frames(u, stop, step):
+    F = RMSF(u.atoms)
+    F.run(stop=stop, step=step)
+    assert F.n_frames == len(range(0, stop, step))