Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
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)