Skip to content

Commit

Permalink
Moves msd results to use the Results class (#3265)
Browse files Browse the repository at this point in the history
Towards #3261

## Work done in this PR
* Moves msd results to use the Results class

Co-authored-by: Oliver Beckstein <[email protected]>
  • Loading branch information
IAlibay and orbeckst authored May 6, 2021
1 parent 503cf1e commit b6ae72f
Show file tree
Hide file tree
Showing 3 changed files with 33 additions and 25 deletions.
2 changes: 2 additions & 0 deletions package/CHANGELOG
Original file line number Diff line number Diff line change
Expand Up @@ -174,6 +174,8 @@ Enhancements
checking if it can be used in parallel analysis. (Issue #2996, PR #2950)

Changes
* `msd.EinsteinMSD` now uses the `analysis.base.Results` class to store
analysis results (Issue #3261)
* `bat.BAT` now uses the `analysis.base.Results` class to store the `bat`
results attribute (Issue #3261, #3272)
* `contacts.Contacts` now stores data using the `Contacts.results.timeseries`
Expand Down
26 changes: 15 additions & 11 deletions package/MDAnalysis/analysis/msd.py
Original file line number Diff line number Diff line change
Expand Up @@ -90,7 +90,7 @@
.. code-block:: python
msd = MSD.timeseries
msd = MSD.results.timeseries
Visual inspection of the MSD is important, so let's take a look at it with a
simple plot.
Expand Down Expand Up @@ -261,16 +261,19 @@ class EinsteinMSD(AnalysisBase):
----------
dim_fac : int
Dimensionality :math:`d` of the MSD.
timeseries : :class:`numpy.ndarray`
results.timeseries : :class:`numpy.ndarray`
The averaged MSD over all the particles with respect to lag-time.
msds_by_particle : :class:`numpy.ndarray`
results.msds_by_particle : :class:`numpy.ndarray`
The MSD of each individual particle with respect to lag-time.
ag : :class:`AtomGroup`
The :class:`AtomGroup` resulting from your selection
n_frames : int
Number of frames included in the analysis.
n_particles : int
Number of particles MSD was calculated over.
.. versionadded:: 2.0.0
"""

def __init__(self, u, select='all', msd_type='xyz', fft=True, **kwargs):
Expand Down Expand Up @@ -307,16 +310,17 @@ def __init__(self, u, select='all', msd_type='xyz', fft=True, **kwargs):
self._position_array = None

# result
self.msds_by_particle = None
self.timeseries = None
self.results.msds_by_particle = None
self.results.timeseries = None

def _prepare(self):
# self.n_frames only available here
# these need to be zeroed prior to each run() call
self.msds_by_particle = np.zeros((self.n_frames, self.n_particles))
self.results.msds_by_particle = np.zeros((self.n_frames,
self.n_particles))
self._position_array = np.zeros(
(self.n_frames, self.n_particles, self.dim_fac))
# self.timeseries not set here
# self.results.timeseries not set here

def _parse_msd_type(self):
r""" Sets up the desired dimensionality of the MSD.
Expand Down Expand Up @@ -360,8 +364,8 @@ def _conclude_simple(self):
for lag in lagtimes:
disp = positions[:-lag, :, :] - positions[lag:, :, :]
sqdist = np.square(disp).sum(axis=-1)
self.msds_by_particle[lag, :] = np.mean(sqdist, axis=0)
self.timeseries = self.msds_by_particle.mean(axis=1)
self.results.msds_by_particle[lag, :] = np.mean(sqdist, axis=0)
self.results.timeseries = self.results.msds_by_particle.mean(axis=1)

def _conclude_fft(self): # with FFT, np.float64 bit prescision required.
r""" Calculates the MSD via the FCA fast correlation algorithm.
Expand All @@ -382,6 +386,6 @@ def _conclude_fft(self): # with FFT, np.float64 bit prescision required.

positions = self._position_array.astype(np.float64)
for n in range(self.n_particles):
self.msds_by_particle[:, n] = tidynamics.msd(
self.results.msds_by_particle[:, n] = tidynamics.msd(
positions[:, n, :])
self.timeseries = self.msds_by_particle.mean(axis=1)
self.results.timeseries = self.results.msds_by_particle.mean(axis=1)
30 changes: 16 additions & 14 deletions testsuite/MDAnalysisTests/analysis/test_msd.py
Original file line number Diff line number Diff line change
Expand Up @@ -123,7 +123,7 @@ def test_simple_step_traj_all_dims(self, step_traj, NSTEP, dim,
m_simple = MSD(step_traj, 'all', msd_type=dim, fft=False)
m_simple.run()
poly = characteristic_poly(NSTEP, dim_factor)
assert_almost_equal(m_simple.timeseries, poly, decimal=4)
assert_almost_equal(m_simple.results.timeseries, poly, decimal=4)

@pytest.mark.parametrize("dim, dim_factor", [
('xyz', 3), ('xy', 2), ('xz', 2), ('yz', 2), ('x', 1), ('y', 1),
Expand All @@ -137,13 +137,14 @@ def test_simple_start_stop_step_all_dims(self, step_traj, NSTEP, dim,
m_simple.run(start=10, stop=1000, step=10)
poly = characteristic_poly(NSTEP, dim_factor)
# polynomial must take offset start into account
assert_almost_equal(m_simple.timeseries, poly[0:990:10], decimal=4)
assert_almost_equal(m_simple.results.timeseries, poly[0:990:10],
decimal=4)

def test_random_walk_u_simple(self, random_walk_u):
# regress against random_walk test data
msd_rw = MSD(random_walk_u, 'all', msd_type='xyz', fft=False)
msd_rw.run()
norm = np.linalg.norm(msd_rw.timeseries)
norm = np.linalg.norm(msd_rw.results.timeseries)
val = 3932.39927487146
assert_almost_equal(norm, val, decimal=5)

Expand All @@ -161,25 +162,25 @@ def msd_fft(self, u, SELECTION):

def test_fft_vs_simple_default(self, msd, msd_fft):
# testing on the PSF, DCD trajectory
timeseries_simple = msd.timeseries
timeseries_fft = msd_fft.timeseries
timeseries_simple = msd.results.timeseries
timeseries_fft = msd_fft.results.timeseries
assert_almost_equal(timeseries_simple, timeseries_fft, decimal=4)

def test_fft_vs_simple_default_per_particle(self, msd, msd_fft):
# check fft and simple give same result per particle
per_particle_simple = msd.msds_by_particle
per_particle_fft = msd_fft.msds_by_particle
per_particle_simple = msd.results.msds_by_particle
per_particle_fft = msd_fft.results.msds_by_particle
assert_almost_equal(per_particle_simple, per_particle_fft, decimal=4)

@pytest.mark.parametrize("dim", ['xyz', 'xy', 'xz', 'yz', 'x', 'y', 'z'])
def test_fft_vs_simple_all_dims(self, u, SELECTION, dim):
# check fft and simple give same result for each dimensionality
m_simple = MSD(u, SELECTION, msd_type=dim, fft=False)
m_simple.run()
timeseries_simple = m_simple.timeseries
timeseries_simple = m_simple.results.timeseries
m_fft = MSD(u, SELECTION, msd_type=dim, fft=True)
m_fft.run()
timeseries_fft = m_fft.timeseries
timeseries_fft = m_fft.results.timeseries
assert_almost_equal(timeseries_simple, timeseries_fft, decimal=4)

@pytest.mark.parametrize("dim", ['xyz', 'xy', 'xz', 'yz', 'x', 'y', 'z'])
Expand All @@ -188,10 +189,10 @@ def test_fft_vs_simple_all_dims_per_particle(self, u, SELECTION, dim):
# dimension
m_simple = MSD(u, SELECTION, msd_type=dim, fft=False)
m_simple.run()
per_particle_simple = m_simple.msds_by_particle
per_particle_simple = m_simple.results.msds_by_particle
m_fft = MSD(u, SELECTION, msd_type=dim, fft=True)
m_fft.run()
per_particle_fft = m_fft.msds_by_particle
per_particle_fft = m_fft.results.msds_by_particle
assert_almost_equal(per_particle_simple, per_particle_fft, decimal=4)

@pytest.mark.parametrize("dim, dim_factor", [
Expand All @@ -208,7 +209,7 @@ def test_fft_step_traj_all_dims(self, step_traj, NSTEP, dim, dim_factor):
m_simple.run()
poly = characteristic_poly(NSTEP, dim_factor)
# this was relaxed from decimal=4 for numpy=1.13 test
assert_almost_equal(m_simple.timeseries, poly, decimal=3)
assert_almost_equal(m_simple.results.timeseries, poly, decimal=3)

@pytest.mark.parametrize("dim, dim_factor", [(
'xyz', 3), ('xy', 2), ('xz', 2), ('yz', 2), ('x', 1), ('y', 1),
Expand All @@ -222,12 +223,13 @@ def test_fft_start_stop_step_all_dims(self, step_traj, NSTEP, dim,
m_simple.run(start=10, stop=1000, step=10)
poly = characteristic_poly(NSTEP, dim_factor)
# polynomial must take offset start into account
assert_almost_equal(m_simple.timeseries, poly[0:990:10], decimal=3)
assert_almost_equal(m_simple.results.timeseries, poly[0:990:10],
decimal=3)

def test_random_walk_u_fft(self, random_walk_u):
# regress against random_walk test data
msd_rw = MSD(random_walk_u, 'all', msd_type='xyz', fft=True)
msd_rw.run()
norm = np.linalg.norm(msd_rw.timeseries)
norm = np.linalg.norm(msd_rw.results.timeseries)
val = 3932.39927487146
assert_almost_equal(norm, val, decimal=5)

0 comments on commit b6ae72f

Please sign in to comment.