Skip to content

Commit

Permalink
convert RMSF to AnalysisBase
Browse files Browse the repository at this point in the history
  • Loading branch information
kain88-de committed Dec 29, 2016
1 parent 98e0e00 commit 01bf8b2
Show file tree
Hide file tree
Showing 2 changed files with 26 additions and 61 deletions.
77 changes: 20 additions & 57 deletions package/MDAnalysis/analysis/rms.py
Original file line number Diff line number Diff line change
Expand Up @@ -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]_.
Expand All @@ -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
10 changes: 6 additions & 4 deletions testsuite/MDAnalysisTests/analysis/test_rms.py
Original file line number Diff line number Diff line change
Expand Up @@ -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?')
Expand Down Expand Up @@ -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)
Expand All @@ -239,16 +241,16 @@ 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,
err_msg="error: rmsf profile should match test " +
"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")
Expand All @@ -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")

0 comments on commit 01bf8b2

Please sign in to comment.