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
5 changes: 5 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
* `contacts.Contacts` now stores data using the `Contacts.results.timeseries`
attribute. `Contacts.timeseries` is now deprecated (Issue #3261)
* `DensityAnalysis` now uses the `results.density` attribute for storing
data. The `DensityAnalysis.density` attribute is now deprecated
(Issue #3261)
Expand Down Expand Up @@ -225,6 +227,9 @@ Changes
* Added OpenMM coordinate and topology converters (Issue #2863, PR #2917)

Deprecations
* The `analysis.Contacts.timeseries` attribute is now deprecated in favour of
`analysis.Contacts.results.timeseries`. It will be removed in 3.0.0
(Issue #3261)
* The `density` attribute of `analysis.density.DensityAnalysis` is now
deprecated in favour of `results.density`. It will be removed in 3.0.0
(Issue #3261)
Expand Down
49 changes: 35 additions & 14 deletions package/MDAnalysis/analysis/contacts.py
Original file line number Diff line number Diff line change
Expand Up @@ -82,11 +82,11 @@
# iterate through trajectory and perform analysis of "native contacts" Q
ca1.run()
# print number of averave contacts
average_contacts = np.mean(ca1.timeseries[:, 1])
average_contacts = np.mean(ca1.results.timeseries[:, 1])
print('average contacts = {}'.format(average_contacts))
# plot time series q(t)
fig, ax = plt.subplots()
ax.plot(ca1.timeseries[:, 0], ca1.timeseries[:, 1])
ax.plot(ca1.results.timeseries[:, 0], ca1.results.timeseries[:, 1])
ax.set(xlabel='frame', ylabel='fraction of native contacts',
title='Native Contacts, average = {:.2f}'.format(average_contacts))
fig.show()
Expand Down Expand Up @@ -123,10 +123,13 @@
q1q2.run()

f, ax = plt.subplots(1, 2, figsize=plt.figaspect(0.5))
ax[0].plot(q1q2.timeseries[:, 0], q1q2.timeseries[:, 1], label='q1')
ax[0].plot(q1q2.timeseries[:, 0], q1q2.timeseries[:, 2], label='q2')
ax[0].plot(q1q2.results.timeseries[:, 0], q1q2.results.timeseries[:, 1],
label='q1')
ax[0].plot(q1q2.results.timeseries[:, 0], q1q2.results.timeseries[:, 2],
label='q2')
ax[0].legend(loc='best')
ax[1].plot(q1q2.timeseries[:, 1], q1q2.timeseries[:, 2], '.-')
ax[1].plot(q1q2.results.timeseries[:, 1],
q1q2.results.timeseries[:, 2], '.-')
f.show()

Compare the resulting pathway to the `MinActionPath result for AdK`_
Expand Down Expand Up @@ -173,8 +176,8 @@ def is_any_closer(r, r0, dist=2.5):
refgroup=(acidic, basic), kwargs={'dist': 2.5})
nc.run()

bound = nc.timeseries[:, 1]
frames = nc.timeseries[:, 0]
bound = nc.results.timeseries[:, 1]
frames = nc.results.timeseries[:, 0]

f, ax = plt.subplots()

Expand Down Expand Up @@ -364,15 +367,25 @@ class Contacts(AnalysisBase):

Attributes
----------
timeseries : list
list containing *Q* for all refgroup pairs and analyzed frames
results.timeseries : numpy.ndarray
2D array containing *Q* for all refgroup pairs and analyzed frames

timeseries : numpy.ndarray
Alias to the :attr:`results.timeseries` attribute.

.. deprecated:: 2.0.0
Will be removed in MDAnalysis 3.0.0. Please use
:attr:`results.timeseries` instead.


.. versionchanged:: 1.0.0
``save()`` method has been removed. Use ``np.savetxt()`` on
:attr:`Contacts.timeseries` instead.
:attr:`Contacts.results.timeseries` instead.
.. versionchanged:: 1.0.0
added ``pbc`` attribute to calculate distances using PBC.

.. versionchanged:: 2.0.0
:attr:`timeseries` results are now stored in a
:class:`MDAnalysis.analysis.base.Results` instance.
"""
def __init__(self, u, select, refgroup, method="hard_cut", radius=4.5,
pbc=True, kwargs=None, **basekwargs):
Expand Down Expand Up @@ -451,10 +464,10 @@ def __init__(self, u, select, refgroup, method="hard_cut", radius=4.5,
self.initial_contacts.append(contact_matrix(self.r0[-1], radius))

def _prepare(self):
self.timeseries = np.empty((self.n_frames, len(self.r0)+1))
self.results.timeseries = np.empty((self.n_frames, len(self.r0)+1))

def _single_frame(self):
self.timeseries[self._frame_index][0] = self._ts.frame
self.results.timeseries[self._frame_index][0] = self._ts.frame

# compute distance array for a frame
d = distance_array(self.grA.positions, self.grB.positions,
Expand All @@ -466,7 +479,15 @@ def _single_frame(self):
r = d[initial_contacts]
r0 = r0[initial_contacts]
q = self.fraction_contacts(r, r0, **self.fraction_kwargs)
self.timeseries[self._frame_index][i] = q
self.results.timeseries[self._frame_index][i] = q

@property
def timeseries(self):
wmsg = ("The `timeseries` attribute was deprecated in MDAnalysis "
"2.0.0 and will be removed in MDAnalysis 3.0.0. Please use "
"`results.timeseries` instead")
warnings.warn(wmsg, DeprecationWarning)
return self.results.timeseries


def _new_selections(u_orig, selections, frame):
Expand Down
35 changes: 23 additions & 12 deletions testsuite/MDAnalysisTests/analysis/test_contacts.py
Original file line number Diff line number Diff line change
Expand Up @@ -20,13 +20,15 @@
# MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations.
# J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787
#
import warnings
import MDAnalysis as mda
import pytest
from MDAnalysis.analysis import contacts
from MDAnalysis.analysis.distances import distance_array

from numpy.testing import (
assert_almost_equal,
assert_equal,
assert_array_equal,
assert_array_almost_equal
)
Expand Down Expand Up @@ -185,18 +187,18 @@ def test_startframe(self, universe):

"""
CA1 = self._run_Contacts(universe)
assert len(CA1.timeseries) == universe.trajectory.n_frames
assert len(CA1.results.timeseries) == universe.trajectory.n_frames

def test_end_zero(self, universe):
"""test_end_zero: TestContactAnalysis1: stop frame 0 is not ignored"""
CA1 = self._run_Contacts(universe, stop=0)
assert len(CA1.timeseries) == 0
assert len(CA1.results.timeseries) == 0

def test_slicing(self, universe):
start, stop, step = 10, 30, 5
CA1 = self._run_Contacts(universe, start=start, stop=stop, step=step)
frames = np.arange(universe.trajectory.n_frames)[start:stop:step]
assert len(CA1.timeseries) == len(frames)
assert len(CA1.results.timeseries) == len(frames)

def test_villin_folded(self):
# one folded, one unfolded
Expand All @@ -213,7 +215,7 @@ def test_villin_folded(self):
q.run()

results = soft_cut(f, u, sel, sel)
assert_almost_equal(q.timeseries[:, 1], results[:, 1])
assert_almost_equal(q.results.timeseries[:, 1], results[:, 1])

def test_villin_unfolded(self):
# both folded
Expand All @@ -230,7 +232,7 @@ def test_villin_unfolded(self):
q.run()

results = soft_cut(f, u, sel, sel)
assert_almost_equal(q.timeseries[:, 1], results[:, 1])
assert_almost_equal(q.results.timeseries[:, 1], results[:, 1])

def test_hard_cut_method(self, universe):
ca = self._run_Contacts(universe)
Expand All @@ -254,8 +256,8 @@ def test_hard_cut_method(self, universe):
0.48543689, 0.44660194, 0.4368932, 0.40776699, 0.41747573,
0.48543689, 0.45631068, 0.46601942, 0.47572816, 0.51456311,
0.45631068, 0.37864078, 0.42718447]
assert len(ca.timeseries) == len(expected)
assert_array_almost_equal(ca.timeseries[:, 1], expected)
assert len(ca.results.timeseries) == len(expected)
assert_array_almost_equal(ca.results.timeseries[:, 1], expected)

def test_radius_cut_method(self, universe):
acidic = universe.select_atoms(self.sel_acidic)
Expand All @@ -268,7 +270,7 @@ def test_radius_cut_method(self, universe):
expected.append(contacts.radius_cut_q(r[initial_contacts], None, radius=6.0))

ca = self._run_Contacts(universe, method='radius_cut')
assert_array_equal(ca.timeseries[:, 1], expected)
assert_array_equal(ca.results.timeseries[:, 1], expected)

@staticmethod
def _is_any_closer(r, r0, dist=2.5):
Expand All @@ -285,7 +287,7 @@ def test_own_method(self, universe):
1., 0., 1., 1., 1., 1., 1., 1., 0., 1., 1., 0., 1.,
0., 0., 1., 1., 0., 0., 1., 1., 1., 0., 1., 0., 0.,
1., 0., 1., 1., 1., 1., 1.]
assert_array_equal(ca.timeseries[:, 1], bound_expected)
assert_array_equal(ca.results.timeseries[:, 1], bound_expected)

@staticmethod
def _weird_own_method(r, r0):
Expand Down Expand Up @@ -315,7 +317,16 @@ def test_distance_box(self, pbc, expected):
r = contacts.Contacts(u, select=(sel_acidic, sel_basic),
refgroup=(acidic, basic), radius=6.0, pbc=pbc)
r.run()
assert_array_almost_equal(r.timeseries[:, 1], expected)
assert_array_almost_equal(r.results.timeseries[:, 1], expected)

def test_warn_deprecated_attr(self, universe):
"""Test for warning message emitted on using deprecated `timeseries`
attribute"""
CA1 = self._run_Contacts(universe, stop=1)
wmsg = "The `timeseries` attribute was deprecated in MDAnalysis"
with pytest.warns(DeprecationWarning, match=wmsg):
assert_equal(CA1.timeseries, CA1.results.timeseries)
Comment on lines +322 to +328
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

good test: check for message and equality



def test_q1q2():
u = mda.Universe(PSF, DCD)
Expand All @@ -342,7 +353,7 @@ def test_q1q2():
0.93097184, 0.93006358, 0.93188011, 0.93278837, 0.93006358,
0.92915531, 0.92824705, 0.92733878, 0.92643052, 0.93188011,
0.93006358, 0.9346049, 0.93188011]
assert_array_almost_equal(q1q2.timeseries[:, 1], q1_expected)
assert_array_almost_equal(q1q2.results.timeseries[:, 1], q1_expected)

q2_expected = [0.94649446, 0.94926199, 0.95295203, 0.95110701, 0.94833948,
0.95479705, 0.94926199, 0.9501845, 0.94926199, 0.95387454,
Expand All @@ -364,4 +375,4 @@ def test_q1q2():
0.97140221, 0.97601476, 0.97693727, 0.98154982, 0.98431734,
0.97601476, 0.9797048, 0.98154982, 0.98062731, 0.98431734,
0.98616236, 0.9898524, 1.]
assert_array_almost_equal(q1q2.timeseries[:, 2], q2_expected)
assert_array_almost_equal(q1q2.results.timeseries[:, 2], q2_expected)