From 01bf8b2c2347a8a5e046efb950213d2dc9e3140b Mon Sep 17 00:00:00 2001 From: Max Linke Date: Thu, 22 Dec 2016 20:57:12 +0100 Subject: [PATCH] convert RMSF to AnalysisBase --- package/MDAnalysis/analysis/rms.py | 77 +++++-------------- .../MDAnalysisTests/analysis/test_rms.py | 10 ++- 2 files changed, 26 insertions(+), 61 deletions(-) diff --git a/package/MDAnalysis/analysis/rms.py b/package/MDAnalysis/analysis/rms.py index 866fc688825..b27eb81f979 100644 --- a/package/MDAnalysis/analysis/rms.py +++ b/package/MDAnalysis/analysis/rms.py @@ -537,32 +537,8 @@ def save(self, filename=None): return filename -class RMSF(object): - """Class to perform RMSF analysis on a set of atoms across a trajectory. - - Run the analysis with :meth:`RMSF.run`, which stores the results - in the array :attr:`RMSF.rmsf`. - - This class performs no coordinate transforms; RMSFs are obtained from atom - coordinates as-is. - - .. versionadded:: 0.11.0 - """ - - def __init__(self, atomgroup): - """Calculate RMSF of given atoms across a trajectory. - - Parameters - ---------- - atomgroup : mda.AtomGroup - AtomGroup to obtain RMSF for - - """ - self.atomgroup = atomgroup - self._rmsf = None - - def run(self, start=None, stop=None, step=None, progout=10, quiet=False): - """Calculate RMSF of given atoms across a trajectory. +class RMSF(AnalysisBase): + """Calculate RMSF of given atoms across a trajectory. This method implements an algorithm for computing sums of squares while avoiding overflows and underflows [Welford1962]_. @@ -589,40 +565,27 @@ def run(self, start=None, stop=None, step=None, progout=10, quiet=False): .. [Welford1962] B. P. Welford (1962). "Note on a Method for Calculating Corrected Sums of Squares and Products." Technometrics 4(3):419-420. - """ - traj = self.atomgroup.universe.trajectory - start, stop, step = traj.check_slice_indices(start, stop, step) - sumsquares = np.zeros((self.atomgroup.n_atoms, 3)) - means = np.array(sumsquares) - - if quiet: - progout = None - - # set up progress output - if progout: - percentage = ProgressMeter(self.atomgroup.universe.trajectory.n_frames, - interval=progout) - else: - percentage = ProgressMeter(self.atomgroup.universe.trajectory.n_frames, - quiet=True) + """ + def __init__(self, atomgroup, weights=None, **kwargs): + super(RMSF, self).__init__(atomgroup.universe.trajectory, **kwargs) + self.atomgroup = atomgroup + if weights == 'mass': + weights = self.atomgroup.masses + self.weights = weights - for k, ts in enumerate(self.atomgroup.universe.trajectory[start:stop:step]): - sumsquares += (k/(k + 1.0)) * (self.atomgroup.positions - means)**2 - means = (k * means + self.atomgroup.positions)/(k + 1) + def _prepare(self): + self.sumsquares = np.zeros((self.atomgroup.n_atoms, 3)) + self.mean = self.sumsquares.copy() - percentage.echo(ts.frame) + def _single_frame(self): + k = self._frame_index + self.sumsquares += (k / (k+1.0)) * (self.atomgroup.positions - self.mean) ** 2 + self.mean = (k * self.mean + self.atomgroup.positions) / (k + 1) - rmsf = np.sqrt(sumsquares.sum(axis=1)/(k + 1)) + def _conclude(self): + k = self._frame_index + self.rmsf = np.sqrt(self.sumsquares.sum(axis=1) / (k + 1)) - if not (rmsf >= 0).all(): + if not (self.rmsf >= 0).all(): raise ValueError("Some RMSF values negative; overflow " + "or underflow occurred") - - self._rmsf = rmsf - - @property - def rmsf(self): - """RMSF data; only available after using :meth:`RMSF.run` - - """ - return self._rmsf diff --git a/testsuite/MDAnalysisTests/analysis/test_rms.py b/testsuite/MDAnalysisTests/analysis/test_rms.py index 6499f97e5e4..49ea4dc9670 100644 --- a/testsuite/MDAnalysisTests/analysis/test_rms.py +++ b/testsuite/MDAnalysisTests/analysis/test_rms.py @@ -140,6 +140,7 @@ def test_with_superposition_equal(self): rmsd_superposition = rms.rmsd(A, B, center=True, superposition=True) assert_almost_equal(rmsd, rmsd_superposition) + class TestRMSD(object): @dec.skipif(parser_not_found('DCD'), 'DCD parser not available. Are you using python 3?') @@ -227,6 +228,7 @@ def test_save_before_run(self): RMSD = MDAnalysis.analysis.rms.RMSD(self.universe) RMSD.save('blah') + class TestRMSF(TestCase): def setUp(self): self.universe = MDAnalysis.Universe(GRO, XTC) @@ -239,7 +241,7 @@ def tearDown(self): def test_rmsf(self): rmsfs = MDAnalysis.analysis.rms.RMSF(self.universe.select_atoms('name CA')) - rmsfs.run(quiet=True) + rmsfs.run() test_rmsfs = np.load(rmsfArray) assert_almost_equal(rmsfs.rmsf, test_rmsfs, 5, @@ -247,8 +249,8 @@ def test_rmsf(self): "values") def test_rmsf_single_frame(self): - rmsfs = MDAnalysis.analysis.rms.RMSF(self.universe.select_atoms('name CA')) - rmsfs.run(start=5, stop=6, quiet=True) + rmsfs = MDAnalysis.analysis.rms.RMSF(self.universe.select_atoms('name CA'), start=5, stop=6) + rmsfs.run() assert_almost_equal(rmsfs.rmsf, 0, 5, err_msg="error: rmsfs should all be zero") @@ -261,7 +263,7 @@ def test_rmsf_identical_frames(self): self.universe = MDAnalysis.Universe(GRO, self.outfile) rmsfs = MDAnalysis.analysis.rms.RMSF(self.universe.select_atoms('name CA')) - rmsfs.run(quiet=True) + rmsfs.run() assert_almost_equal(rmsfs.rmsf, 0, 5, err_msg="error: rmsfs should all be 0")