From 0ca2f0f156f4ffd3e71d4947d45bd04615feb791 Mon Sep 17 00:00:00 2001 From: VOD555 Date: Wed, 10 Oct 2018 09:59:02 -0700 Subject: [PATCH 01/15] Add InterRDF_s --- pmda/rdf.py | 114 ++++++++++++++++++++++++++++++++++++ pmda/test/test_rdf_s.py | 126 ++++++++++++++++++++++++++++++++++++++++ 2 files changed, 240 insertions(+) create mode 100644 pmda/test/test_rdf_s.py diff --git a/pmda/rdf.py b/pmda/rdf.py index 40aebc39..5eb2bca4 100644 --- a/pmda/rdf.py +++ b/pmda/rdf.py @@ -172,3 +172,117 @@ def _reduce(res, result_single_frame): # Add two numpy arrays res += result_single_frame return res + + +class InterRDF_s(ParallelAnalysisBase): + def __init__(self, u, ags, + nbins=75, range=(0.0, 15.0), density=True, **kwargs): + atomgroups = [] + for pair in ags: + atomgroups.append(pair[0]) + atomgroups.append(pair[1]) + super(InterRDF_s, self).__init__(u, atomgroups) + + # List of pairs of AtomGroups + self.ags = ags + self._density = density + self.n = len(ags) + self.nf = u.trajectory.n_frames + self.rdf_settings = {'bins': nbins, + 'range': range} + + def _prepare(self): + # Empty list to store the RDF + count_list = [] + count, edges = np.histogram([-1], **self.rdf_settings) + self.len = len(count) + self.edges = edges + self.bins = 0.5 * (edges[:-1] + edges[1:]) + + # Need to know average volume + self.volume = 0.0 + self._maxrange = self.rdf_settings['range'][1] + + + def _single_frame(self, ts, atomgroups): + ags = [[atomgroups[2*i], atomgroups[2*i+1]] for i in range(self.n)] + count = [np.zeros((ag1.n_atoms, ag2.n_atoms, self.len), dtype=np.float64) + for ag1, ag2 in ags] + for i, (ag1, ag2) in enumerate(ags): + u = ag1.universe + pairs, dist = distances.capped_distance(ag1.positions, + ag2.positions, + self._maxrange, + box=u.dimensions) + + for j, (idx1, idx2) in enumerate(pairs): + count[i][idx1, idx2, :] = np.histogram(dist[j], + **self.rdf_settings)[0] + + volume = u.trajectory.ts.volume + + return np.array([np.array(count), np.array(volume, dtype=np.float64)]) + + def _conclude(self): + self.count = np.sum(self._results[:, 0]) + self.volume = np.sum(self._results[:, 1]) + # Volume in each radial shell + vol = np.power(self.edges[1:], 3) - np.power(self.edges[:-1], 3) + vol *= 4/3.0 * np.pi + + # Empty lists to restore indices, RDF + indices = [] + rdf = [] + + for i, (ag1, ag2) in enumerate(self.ags): + # Number of each selection + nA = len(ag1) + nB = len(ag2) + N = nA * nB + indices.append([ag1.indices, ag2.indices]) + + # Average number density + box_vol = self.volume / self.nf + density = N / box_vol + + if self._density: + rdf.append(self.count[i] / (density * vol * self.nf)) + else: + rdf.append(self.count[i] / (vol * self.nf)) + + self.rdf = rdf + self.indices = indices + + def get_cdf(self): + """Calculate the cumulative distribution functions (CDF) for all sites. + Note that this is the actual count within a given radius, i.e., + :math:`N(r)`. + Returns + ------- + cdf : list + list of arrays with the same structure as :attr:`rdf` + """ + # Calculate cumulative distribution function + # Empty list to restore CDF + cdf = [] + + for count in self.count: + cdf.append(np.cumsum(count, axis=2) / self.nf) + + # Results stored in self.cdf + # self.cdf is a list of cdf between pairs of AtomGroups in ags + self.cdf = cdf + + return cdf + + @staticmethod + def _reduce(res, result_single_frame): + """ 'add' action for an accumulator""" + if isinstance(res, list) and len(res) == 0: + # Convert res from an empty list to a numpy array + # which has the same shape as the single frame result + res = result_single_frame + else: + # Add two numpy arrays + res += result_single_frame + return res diff --git a/pmda/test/test_rdf_s.py b/pmda/test/test_rdf_s.py new file mode 100644 index 00000000..1b711e22 --- /dev/null +++ b/pmda/test/test_rdf_s.py @@ -0,0 +1,126 @@ +# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding:utf-8 -*- +# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 fileencoding=utf-8 +# +# MDAnalysis --- https://www.mdanalysis.org +# Copyright (c) 2006-2018 The MDAnalysis Development Team and contributors +# (see the file AUTHORS for the full list of names) +# +# Released under the GNU Public Licence, v2 or any higher version +# +# Please cite your use of MDAnalysis in published work: +# +# R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy, M. N. Melo, S. L. Seyler, +# D. L. Dotson, J. Domanski, S. Buchoux, I. M. Kenney, and O. Beckstein. +# MDAnalysis: A Python package for the rapid analysis of molecular dynamics +# simulations. In S. Benthall and S. Rostrup editors, Proceedings of the 15th +# Python in Science Conference, pages 102-109, Austin, TX, 2016. SciPy. +# +# N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein. +# MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations. +# J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787 +# +from __future__ import absolute_import, division + +import pytest + +from numpy.testing import assert_almost_equal + +import MDAnalysis as mda +import numpy as np +from pmda.rdf import InterRDF_s +from MDAnalysis.analysis import rdf + +from MDAnalysisTests.datafiles import GRO_MEMPROT, XTC_MEMPROT + + +@pytest.fixture(scope='module') +def u(): + return mda.Universe(GRO_MEMPROT, XTC_MEMPROT) + + +@pytest.fixture(scope='module') +def sels(u): + s1 = u.select_atoms('name ZND and resid 289') + s2 = u.select_atoms('(name OD1 or name OD2) and resid 51 and sphzone 5.0 (resid 289)') + s3 = u.select_atoms('name ZND and (resid 291 or resid 292)') + s4 = u.select_atoms('(name OD1 or name OD2) and sphzone 5.0 (resid 291)') + ags = [[s1, s2], [s3, s4]] + return ags + +@pytest.fixture(scope='module') +def rdf_s(u, sels): + return InterRDF_s(u, sels).run() + +def test_nbins(u): + ags = sels(u) + rdf = InterRDF_s(u, ags, nbins=412).run() + + assert len(rdf.bins) == 412 + + +def test_range(u): + ags = sels(u) + rmin, rmax = 1.0, 13.0 + rdf = InterRDF_s(u, ags, range=(rmin, rmax)).run() + + assert rdf.edges[0] == rmin + assert rdf.edges[-1] == rmax + + +def test_count_size(rdf_s): + # ZND vs OD1 & OD2 + # should see 2 elements in rdf.count + # 1 element in rdf.count[0] + # 2 elements in rdf.count[0][0] + # 2 elements in rdf.count[1] + # 2 elements in rdf.count[1][0] + # 2 elements in rdf.count[1][1] + assert len(rdf_s.count) == 2 + assert len(rdf_s.count[0]) == 1 + assert len(rdf_s.count[0][0]) == 2 + assert len(rdf_s.count[1]) == 2 + assert len(rdf_s.count[1][0]) == 2 + assert len(rdf_s.count[1][1]) == 2 + + +def test_count(rdf_s): + # should see one distance with 5 counts in count[0][0][1] + # should see one distance with 3 counts in count[1][1][0] + assert len(rdf_s.count[0][0][1][rdf_s.count[0][0][1] == 5]) == 1 + assert len(rdf_s.count[1][1][0][rdf_s.count[1][1][0] == 3]) == 1 + + +def test_double_run(rdf_s): + # running rdf twice should give the same result + assert len(rdf_s.count[0][0][1][rdf_s.count[0][0][1] == 5]) == 1 + assert len(rdf_s.count[1][1][0][rdf_s.count[1][1][0] == 3]) == 1 + +def test_cdf(rdf_s): + rdf_s.get_cdf() + assert rdf_s.cdf[0][0][0][-1] == rdf_s.count[0][0][0].sum()/rdf_s.nf + + +def test_reduce(rdf_s): + # should see numpy.array addtion + res = [] + single_frame = np.array([np.array([1, 2]), np.array([3])]) + res = rdf_s._reduce(res, single_frame) + res = rdf_s._reduce(res, single_frame) + assert_almost_equal(res[0], np.array([2, 4])) + assert_almost_equal(res[1], np.array([6])) + + +def test_same_result(u, sels): + # should see same results from analysis.rdf.InterRDF_s and pmda.rdf.InterRDF_s + nrdf = rdf.InterRDF_s(u, sels).run() + prdf = InterRDF_s(u, sels).run(n_jobs=3) + assert_almost_equal(nrdf.count[0][0][0], prdf.count[0][0][0]) + assert_almost_equal(nrdf.rdf[0][0][0], prdf.rdf[0][0][0]) + + +@pytest.mark.parametrize("density, value", [ + (True, 13275.775528444701), + (False, 0.021915460340071267)]) +def test_density(u, sels, density, value): + rdf = InterRDF_s(u, sels, density=density).run() + assert_almost_equal(max(rdf.rdf[0][0][0]), value) From 3582cba2f77cca1813ff46d1822b3145b2219f81 Mon Sep 17 00:00:00 2001 From: VOD555 Date: Wed, 10 Oct 2018 12:29:23 -0700 Subject: [PATCH 02/15] add docs and fix some pep8 failures --- pmda/rdf.py | 69 +++++++++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 64 insertions(+), 5 deletions(-) diff --git a/pmda/rdf.py b/pmda/rdf.py index 5eb2bca4..86bba826 100644 --- a/pmda/rdf.py +++ b/pmda/rdf.py @@ -175,8 +175,69 @@ def _reduce(res, result_single_frame): class InterRDF_s(ParallelAnalysisBase): + """Site-specific intermolecular pair distribution function + + Arguments + --------- + u : Universe + a Universe that contains atoms in `ags` + ags : list + a list of pairs of :class:`~MDAnalysis.core.groups.AtomGroup` + instances + nbins : int (optional) + Number of bins in the histogram [75] + range : tuple or list (optional) + The size of the RDF [0.0, 15.0] + + Example + ------- + + First create the :class:`InterRDF_s` object, by supplying one Universe and + one list of pairs of AtomGroups, then use the :meth:`~InterRDF_s.run` + method:: + + from MDAnalysisTests.datafiles import GRO_MEMPROT, XTC_MEMPROT + u = mda.Universe(GRO_MEMPROT, XTC_MEMPROT) + + s1 = u.select_atoms('name ZND and resid 289') + s2 = u.select_atoms('(name OD1 or name OD2) and resid 51 and sphzone 5.0 + (resid 289)') + s3 = u.select_atoms('name ZND and (resid 291 or resid 292)') + s4 = u.select_atoms('(name OD1 or name OD2) and sphzone 5.0 (resid 291)') + ags = [[s1, s2], [s3, s4]] + + rdf = InterRDF_s(u, ags) + rdf.run() + + Results are available through the :attr:`bins` and :attr:`rdf` attributes:: + + plt.plot(rdf.bins, rdf.rdf[0][0][0]) + + (Which plots the rdf between the first atom in ``s1`` and the first atom in + ``s2``) + + To generate the *cumulative distribution function* (cdf), use the + :meth:`~InterRDF_s.get_cdf` method :: + + cdf = rdf.get_cdf() + + Results are available through the :attr:'cdf' attribute:: + + plt.plot(rdf.bins, rdf.cdf[0][0][0]) + + (Which plots the cdf between the first atom in ``s1`` and the first atom in + ``s2``) + + See Also + -------- + MDAnalysis.analysis.rdf.InterRDF_s + + + .. versionadded:: 0.2.0 + """ + def __init__(self, u, ags, - nbins=75, range=(0.0, 15.0), density=True, **kwargs): + nbins=75, range=(0.0, 15.0), density=True): atomgroups = [] for pair in ags: atomgroups.append(pair[0]) @@ -193,7 +254,6 @@ def __init__(self, u, ags, def _prepare(self): # Empty list to store the RDF - count_list = [] count, edges = np.histogram([-1], **self.rdf_settings) self.len = len(count) self.edges = edges @@ -203,11 +263,10 @@ def _prepare(self): self.volume = 0.0 self._maxrange = self.rdf_settings['range'][1] - def _single_frame(self, ts, atomgroups): ags = [[atomgroups[2*i], atomgroups[2*i+1]] for i in range(self.n)] - count = [np.zeros((ag1.n_atoms, ag2.n_atoms, self.len), dtype=np.float64) - for ag1, ag2 in ags] + count = [np.zeros((ag1.n_atoms, ag2.n_atoms, self.len), + dtype=np.float64) for ag1, ag2 in ags] for i, (ag1, ag2) in enumerate(ags): u = ag1.universe pairs, dist = distances.capped_distance(ag1.positions, From 396b7582eaefe10bb6841239a0a37908c1aba463 Mon Sep 17 00:00:00 2001 From: VOD555 Date: Thu, 25 Oct 2018 17:43:22 -0700 Subject: [PATCH 03/15] fix pylint error --- pmda/rdf.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/pmda/rdf.py b/pmda/rdf.py index 86bba826..63842796 100644 --- a/pmda/rdf.py +++ b/pmda/rdf.py @@ -236,6 +236,9 @@ class InterRDF_s(ParallelAnalysisBase): .. versionadded:: 0.2.0 """ + # pylint: disable=redefined-builtin + # continue to use 'range' as long as MDAnalysis uses it so that + # the user interface remains consistent def __init__(self, u, ags, nbins=75, range=(0.0, 15.0), density=True): atomgroups = [] @@ -252,6 +255,8 @@ def __init__(self, u, ags, self.rdf_settings = {'bins': nbins, 'range': range} + # pylint: enable=redefined-builtin + def _prepare(self): # Empty list to store the RDF count, edges = np.histogram([-1], **self.rdf_settings) From e56c5f77e7b19d8a7614bff39118fd50e4fc00fd Mon Sep 17 00:00:00 2001 From: VOD555 Date: Thu, 25 Oct 2018 17:43:53 -0700 Subject: [PATCH 04/15] fix pep-8 --- pmda/test/test_rdf_s.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/pmda/test/test_rdf_s.py b/pmda/test/test_rdf_s.py index 1b711e22..df1845cb 100644 --- a/pmda/test/test_rdf_s.py +++ b/pmda/test/test_rdf_s.py @@ -41,16 +41,19 @@ def u(): @pytest.fixture(scope='module') def sels(u): s1 = u.select_atoms('name ZND and resid 289') - s2 = u.select_atoms('(name OD1 or name OD2) and resid 51 and sphzone 5.0 (resid 289)') + s2 = u.select_atoms('(name OD1 or name OD2) and resid 51 and + sphzone 5.0 (resid 289)') s3 = u.select_atoms('name ZND and (resid 291 or resid 292)') s4 = u.select_atoms('(name OD1 or name OD2) and sphzone 5.0 (resid 291)') ags = [[s1, s2], [s3, s4]] return ags + @pytest.fixture(scope='module') def rdf_s(u, sels): return InterRDF_s(u, sels).run() + def test_nbins(u): ags = sels(u) rdf = InterRDF_s(u, ags, nbins=412).run() @@ -95,6 +98,7 @@ def test_double_run(rdf_s): assert len(rdf_s.count[0][0][1][rdf_s.count[0][0][1] == 5]) == 1 assert len(rdf_s.count[1][1][0][rdf_s.count[1][1][0] == 3]) == 1 + def test_cdf(rdf_s): rdf_s.get_cdf() assert rdf_s.cdf[0][0][0][-1] == rdf_s.count[0][0][0].sum()/rdf_s.nf @@ -111,7 +115,8 @@ def test_reduce(rdf_s): def test_same_result(u, sels): - # should see same results from analysis.rdf.InterRDF_s and pmda.rdf.InterRDF_s + # should see same results from analysis.rdf.InterRDF_s + # and pmda.rdf.InterRDF_s nrdf = rdf.InterRDF_s(u, sels).run() prdf = InterRDF_s(u, sels).run(n_jobs=3) assert_almost_equal(nrdf.count[0][0][0], prdf.count[0][0][0]) From b17d73705dced39b9333d269f520f5db99c0b8ec Mon Sep 17 00:00:00 2001 From: VOD555 Date: Fri, 26 Oct 2018 18:58:12 -0700 Subject: [PATCH 05/15] fix pep8 error in test_rdf_s.py --- pmda/test/test_rdf_s.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pmda/test/test_rdf_s.py b/pmda/test/test_rdf_s.py index df1845cb..8fa00734 100644 --- a/pmda/test/test_rdf_s.py +++ b/pmda/test/test_rdf_s.py @@ -41,8 +41,8 @@ def u(): @pytest.fixture(scope='module') def sels(u): s1 = u.select_atoms('name ZND and resid 289') - s2 = u.select_atoms('(name OD1 or name OD2) and resid 51 and - sphzone 5.0 (resid 289)') + s2 = u.select_atoms( + '(name OD1 or name OD2) and resid 51 and sphzone 5.0 (resid 289)') s3 = u.select_atoms('name ZND and (resid 291 or resid 292)') s4 = u.select_atoms('(name OD1 or name OD2) and sphzone 5.0 (resid 291)') ags = [[s1, s2], [s3, s4]] From ea65bc4eba11ebd81a15d8ad8d102c224f7fe568 Mon Sep 17 00:00:00 2001 From: VOD555 Date: Fri, 26 Oct 2018 20:00:22 -0700 Subject: [PATCH 06/15] change docs --- pmda/test/test_rdf.py | 19 +++---------------- pmda/test/test_rdf_s.py | 19 +++---------------- 2 files changed, 6 insertions(+), 32 deletions(-) diff --git a/pmda/test/test_rdf.py b/pmda/test/test_rdf.py index 2dfec024..889af140 100644 --- a/pmda/test/test_rdf.py +++ b/pmda/test/test_rdf.py @@ -1,24 +1,11 @@ # -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding:utf-8 -*- -# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 fileencoding=utf-8 +# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 # -# MDAnalysis --- https://www.mdanalysis.org -# Copyright (c) 2006-2017 The MDAnalysis Development Team and contributors +# PMDA +# Copyright (c) 2017 The MDAnalysis Development Team and contributors # (see the file AUTHORS for the full list of names) # # Released under the GNU Public Licence, v2 or any higher version -# -# Please cite your use of MDAnalysis in published work: -# -# R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy, M. N. Melo, S. L. Seyler, -# D. L. Dotson, J. Domanski, S. Buchoux, I. M. Kenney, and O. Beckstein. -# MDAnalysis: A Python package for the rapid analysis of molecular dynamics -# simulations. In S. Benthall and S. Rostrup editors, Proceedings of the 15th -# Python in Science Conference, pages 102-109, Austin, TX, 2016. SciPy. -# -# N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein. -# MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations. -# J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787 -# from __future__ import absolute_import import pytest diff --git a/pmda/test/test_rdf_s.py b/pmda/test/test_rdf_s.py index 8fa00734..f052e64f 100644 --- a/pmda/test/test_rdf_s.py +++ b/pmda/test/test_rdf_s.py @@ -1,24 +1,11 @@ # -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding:utf-8 -*- -# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 fileencoding=utf-8 +# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 # -# MDAnalysis --- https://www.mdanalysis.org -# Copyright (c) 2006-2018 The MDAnalysis Development Team and contributors +# PMDA +# Copyright (c) 2017 The MDAnalysis Development Team and contributors # (see the file AUTHORS for the full list of names) # # Released under the GNU Public Licence, v2 or any higher version -# -# Please cite your use of MDAnalysis in published work: -# -# R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy, M. N. Melo, S. L. Seyler, -# D. L. Dotson, J. Domanski, S. Buchoux, I. M. Kenney, and O. Beckstein. -# MDAnalysis: A Python package for the rapid analysis of molecular dynamics -# simulations. In S. Benthall and S. Rostrup editors, Proceedings of the 15th -# Python in Science Conference, pages 102-109, Austin, TX, 2016. SciPy. -# -# N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein. -# MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations. -# J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787 -# from __future__ import absolute_import, division import pytest From 90242231624f2b287047868719510397decaa145 Mon Sep 17 00:00:00 2001 From: VOD555 Date: Sun, 28 Oct 2018 21:45:43 -0700 Subject: [PATCH 07/15] edit format --- pmda/rdf.py | 1 + 1 file changed, 1 insertion(+) diff --git a/pmda/rdf.py b/pmda/rdf.py index 63842796..65a9142a 100644 --- a/pmda/rdf.py +++ b/pmda/rdf.py @@ -321,6 +321,7 @@ def get_cdf(self): """Calculate the cumulative distribution functions (CDF) for all sites. Note that this is the actual count within a given radius, i.e., :math:`N(r)`. + Returns ------- cdf : list From df874380a98c6f5a813ea1f4c7b81aa0bf094d44 Mon Sep 17 00:00:00 2001 From: VOD555 Date: Sun, 28 Oct 2018 22:10:55 -0700 Subject: [PATCH 08/15] add test for different n_blocks --- pmda/test/test_rdf_s.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/pmda/test/test_rdf_s.py b/pmda/test/test_rdf_s.py index f052e64f..a98db78a 100644 --- a/pmda/test/test_rdf_s.py +++ b/pmda/test/test_rdf_s.py @@ -101,11 +101,12 @@ def test_reduce(rdf_s): assert_almost_equal(res[1], np.array([6])) -def test_same_result(u, sels): +@pytest.mark.parametrize("n_blocks", [1, 2]) +def test_same_result(u, sels, n_blocks): # should see same results from analysis.rdf.InterRDF_s # and pmda.rdf.InterRDF_s nrdf = rdf.InterRDF_s(u, sels).run() - prdf = InterRDF_s(u, sels).run(n_jobs=3) + prdf = InterRDF_s(u, sels).run(n_blocks=n_blocks) assert_almost_equal(nrdf.count[0][0][0], prdf.count[0][0][0]) assert_almost_equal(nrdf.rdf[0][0][0], prdf.rdf[0][0][0]) From 9c4999166be62b0ced7d13e838d2be307810d910 Mon Sep 17 00:00:00 2001 From: VOD555 Date: Sun, 28 Oct 2018 22:24:49 -0700 Subject: [PATCH 09/15] add to CHANGELOG --- CHANGELOG | 3 +++ 1 file changed, 3 insertions(+) diff --git a/CHANGELOG b/CHANGELOG index a7c18fd2..50b572b3 100644 --- a/CHANGELOG +++ b/CHANGELOG @@ -22,6 +22,9 @@ Enhancements * add parallel particle-particle RDF calculation module pmda.rdf (Issue #41) * add readonly_attributes context manager to ParallelAnalysisBase * add parallel implementation of Leaflet Finder (Issue #47) + * add parallel site-specific RDF calculation module pmda.rdf.InterRDF_s + (Issue #60) + Fixes * always distribute frames over blocks so that no empty blocks are From f384984fb42f78a368be42e9840c0102a0f64db5 Mon Sep 17 00:00:00 2001 From: VOD555 Date: Wed, 31 Oct 2018 17:31:05 -0700 Subject: [PATCH 10/15] dask update --- conftest.py | 2 +- pmda/parallel.py | 8 ++++---- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/conftest.py b/conftest.py index a2b4ee6b..6b11f00d 100644 --- a/conftest.py +++ b/conftest.py @@ -29,4 +29,4 @@ def scheduler(request, client): if request.param == 'distributed': return client else: - return multiprocessing + return request.param diff --git a/pmda/parallel.py b/pmda/parallel.py index 06763131..3f92b2c4 100644 --- a/pmda/parallel.py +++ b/pmda/parallel.py @@ -293,13 +293,13 @@ def run(self, """ if scheduler is None: - scheduler = multiprocessing + scheduler = 'multiprocessing' if n_jobs == -1: n_jobs = cpu_count() if n_blocks is None: - if scheduler == multiprocessing: + if scheduler == 'multiprocessing': n_blocks = n_jobs elif isinstance(scheduler, distributed.Client): n_blocks = len(scheduler.ncores()) @@ -308,8 +308,8 @@ def run(self, "Couldn't guess ideal number of blocks from scheduler." "Please provide `n_blocks` in call to method.") - scheduler_kwargs = {'get': scheduler.get} - if scheduler == multiprocessing: + scheduler_kwargs = {'scheduler': scheduler} + if scheduler == 'multiprocessing': scheduler_kwargs['num_workers'] = n_jobs start, stop, step = self._trajectory.check_slice_indices( From 44b5549c744f2512cf236be1c14ca20cc26dd4b0 Mon Sep 17 00:00:00 2001 From: VOD555 Date: Thu, 1 Nov 2018 14:19:46 -0700 Subject: [PATCH 11/15] add tests for different n_blocks --- pmda/test/test_rdf.py | 5 +++-- pmda/test/test_rdf_s.py | 2 +- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/pmda/test/test_rdf.py b/pmda/test/test_rdf.py index 889af140..aeddf3bd 100644 --- a/pmda/test/test_rdf.py +++ b/pmda/test/test_rdf.py @@ -73,11 +73,12 @@ def test_double_run(sels): assert len(rdf.count[rdf.count == 3]) == 7 -def test_same_result(sels): +@pytest.mark.parametrize("n_blocks", [1, 2, 3, 4]) +def test_same_result(sels, n_blocks): # should see same results from analysis.rdf and pmda.rdf s1, s2 = sels nrdf = rdf.InterRDF(s1, s2).run() - prdf = InterRDF(s1, s2).run() + prdf = InterRDF(s1, s2).run(n_blocks=n_blocks) assert_almost_equal(nrdf.count, prdf.count) assert_almost_equal(nrdf.rdf, prdf.rdf) diff --git a/pmda/test/test_rdf_s.py b/pmda/test/test_rdf_s.py index a98db78a..184147c5 100644 --- a/pmda/test/test_rdf_s.py +++ b/pmda/test/test_rdf_s.py @@ -101,7 +101,7 @@ def test_reduce(rdf_s): assert_almost_equal(res[1], np.array([6])) -@pytest.mark.parametrize("n_blocks", [1, 2]) +@pytest.mark.parametrize("n_blocks", [1, 2, 3, 4]) def test_same_result(u, sels, n_blocks): # should see same results from analysis.rdf.InterRDF_s # and pmda.rdf.InterRDF_s From eea692f31489d3e36e68454af8b2853fdefd9448 Mon Sep 17 00:00:00 2001 From: VOD555 Date: Thu, 1 Nov 2018 14:56:34 -0700 Subject: [PATCH 12/15] add scheduler test --- pmda/test/test_rdf_s.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pmda/test/test_rdf_s.py b/pmda/test/test_rdf_s.py index 184147c5..08e2d8f9 100644 --- a/pmda/test/test_rdf_s.py +++ b/pmda/test/test_rdf_s.py @@ -37,7 +37,7 @@ def sels(u): @pytest.fixture(scope='module') -def rdf_s(u, sels): +def rdf_s(u, sels, scheduler): return InterRDF_s(u, sels).run() From d238bb4a83b6e633d344e3176e346e991a227a55 Mon Sep 17 00:00:00 2001 From: VOD555 Date: Thu, 1 Nov 2018 15:02:27 -0700 Subject: [PATCH 13/15] add scheduler test for test_rdf --- pmda/test/test_rdf.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pmda/test/test_rdf.py b/pmda/test/test_rdf.py index aeddf3bd..7d33f58e 100644 --- a/pmda/test/test_rdf.py +++ b/pmda/test/test_rdf.py @@ -50,7 +50,7 @@ def test_range(u): assert rdf.edges[-1] == rmax -def test_count_sum(sels): +def test_count_sum(sels, scheduler): # OD1 vs OD2 # should see 577 comparisons in count s1, s2 = sels From f09b024f3b11729fa5d0ba385314e26e9285ca29 Mon Sep 17 00:00:00 2001 From: VOD555 Date: Thu, 1 Nov 2018 15:41:36 -0700 Subject: [PATCH 14/15] change method to pass information of ags --- pmda/rdf.py | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/pmda/rdf.py b/pmda/rdf.py index 65a9142a..6eccc543 100644 --- a/pmda/rdf.py +++ b/pmda/rdf.py @@ -248,12 +248,19 @@ def __init__(self, u, ags, super(InterRDF_s, self).__init__(u, atomgroups) # List of pairs of AtomGroups - self.ags = ags + #self.ags = ags self._density = density self.n = len(ags) self.nf = u.trajectory.n_frames self.rdf_settings = {'bins': nbins, 'range': range} + ag_shape = [] + indices = [] + for (ag1, ag2) in ags: + indices.append([ag1.indices, ag2.indices]) + ag_shape.append([len(ag1), len(ag2)]) + self.indices = indices + self.ag_shape =ag_shape # pylint: enable=redefined-builtin @@ -295,15 +302,13 @@ def _conclude(self): vol *= 4/3.0 * np.pi # Empty lists to restore indices, RDF - indices = [] rdf = [] - for i, (ag1, ag2) in enumerate(self.ags): + for i, (nA, nB) in enumerate(self.ag_shape): # Number of each selection - nA = len(ag1) - nB = len(ag2) + #nA = len(ag1) + #nB = len(ag2) N = nA * nB - indices.append([ag1.indices, ag2.indices]) # Average number density box_vol = self.volume / self.nf @@ -315,7 +320,6 @@ def _conclude(self): rdf.append(self.count[i] / (vol * self.nf)) self.rdf = rdf - self.indices = indices def get_cdf(self): """Calculate the cumulative distribution functions (CDF) for all sites. From 60a000d8eefbd564e6db00804dc09bea68b5a715 Mon Sep 17 00:00:00 2001 From: VOD555 Date: Thu, 1 Nov 2018 15:59:08 -0700 Subject: [PATCH 15/15] fix pep8 error --- pmda/rdf.py | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/pmda/rdf.py b/pmda/rdf.py index 6eccc543..c6f85bd4 100644 --- a/pmda/rdf.py +++ b/pmda/rdf.py @@ -248,7 +248,6 @@ def __init__(self, u, ags, super(InterRDF_s, self).__init__(u, atomgroups) # List of pairs of AtomGroups - #self.ags = ags self._density = density self.n = len(ags) self.nf = u.trajectory.n_frames @@ -260,7 +259,7 @@ def __init__(self, u, ags, indices.append([ag1.indices, ag2.indices]) ag_shape.append([len(ag1), len(ag2)]) self.indices = indices - self.ag_shape =ag_shape + self.ag_shape = ag_shape # pylint: enable=redefined-builtin @@ -306,8 +305,6 @@ def _conclude(self): for i, (nA, nB) in enumerate(self.ag_shape): # Number of each selection - #nA = len(ag1) - #nB = len(ag2) N = nA * nB # Average number density