diff --git a/pmda/density.py b/pmda/density.py index 4623f7d7..1b158ea3 100644 --- a/pmda/density.py +++ b/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'] diff --git a/pmda/parallel.py b/pmda/parallel.py index 18faaf16..bc41ef33 100644 --- a/pmda/parallel.py +++ b/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: diff --git a/pmda/rms/rmsf.py b/pmda/rms/rmsf.py index 70f7013d..a8bc2969 100644 --- a/pmda/rms/rmsf.py +++ b/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 diff --git a/pmda/test/test_density.py b/pmda/test/test_density.py index 02e9f4bb..2e66e708 100644 --- a/pmda/test/test_density.py +++ b/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)) diff --git a/pmda/test/test_rmsf.py b/pmda/test/test_rmsf.py index 63eca813..58e62cd6 100644 --- a/pmda/test/test_rmsf.py +++ b/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))