From 1480ac7a062573099fcc786dab977e9dd86f72fa Mon Sep 17 00:00:00 2001 From: Max Linke Date: Wed, 6 Jul 2016 18:00:49 +0200 Subject: [PATCH 01/10] PEP8 Changes --- package/MDAnalysis/lib/log.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/package/MDAnalysis/lib/log.py b/package/MDAnalysis/lib/log.py index 3f8d6eba708..d4b8638e852 100644 --- a/package/MDAnalysis/lib/log.py +++ b/package/MDAnalysis/lib/log.py @@ -2,8 +2,8 @@ # vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 # # MDAnalysis --- http://www.MDAnalysis.org -# Copyright (c) 2006-2015 Naveen Michaud-Agrawal, Elizabeth J. Denning, Oliver Beckstein -# and contributors (see AUTHORS for the full list) +# Copyright (c) 2006-2015 Naveen Michaud-Agrawal, Elizabeth J. Denning, Oliver +# Beckstein and contributors (see AUTHORS for the full list) # # Released under the GNU Public Licence, v2 or any higher version # @@ -78,6 +78,7 @@ from .. import version + def start_logging(logfile="MDAnalysis.log", version=version.__version__): """Start logging of messages to file and console. @@ -96,7 +97,6 @@ def stop_logging(): clear_handlers(logger) # this _should_ do the job... - def create(logger_name="MDAnalysis", logfile="MDAnalysis.log"): """Create a top level logger. @@ -209,7 +209,8 @@ class ProgressMeter(object): """ - def __init__(self, numsteps, format=None, interval=10, offset=1, quiet=False): + def __init__(self, numsteps, format=None, interval=10, offset=1, + quiet=False): """Set up the ProgressMeter :Arguments: From 095f9b0ad0fee4d8d61a1f9eaca9b3d389253165 Mon Sep 17 00:00:00 2001 From: Max Linke Date: Wed, 6 Jul 2016 18:00:55 +0200 Subject: [PATCH 02/10] Add progressmeter to AnalysisBase This also adds an initialization function to AnalysisBase. This is done now because we need to do more then just to setup the iteration. --- package/MDAnalysis/analysis/base.py | 120 +++++++++++++----- .../MDAnalysisTests/analysis/test_base.py | 43 ++++--- 2 files changed, 108 insertions(+), 55 deletions(-) diff --git a/package/MDAnalysis/analysis/base.py b/package/MDAnalysis/analysis/base.py index 5c6a3dad112..cf99e2bb669 100644 --- a/package/MDAnalysis/analysis/base.py +++ b/package/MDAnalysis/analysis/base.py @@ -2,8 +2,8 @@ # vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 # # MDAnalysis --- http://www.MDAnalysis.org -# Copyright (c) 2006-2015 Naveen Michaud-Agrawal, Elizabeth J. Denning, Oliver Beckstein -# and contributors (see AUTHORS for the full list) +# Copyright (c) 2006-2015 Naveen Michaud-Agrawal, Elizabeth J. Denning, Oliver +# Beckstein and contributors (see AUTHORS for the full list) # # Released under the GNU Public Licence, v2 or any higher version # @@ -24,7 +24,6 @@ """ from six.moves import range -import numpy as np import logging from MDAnalysis.lib.log import ProgressMeter @@ -32,39 +31,90 @@ class AnalysisBase(object): - """Base class for defining multi frame analysis - - The analysis base class is designed as a template for creating - multiframe analysis. - - The class implements the following methods: - - _setup_frames(trajectory, start=None, stop=None, step=None) - Pass a Reader object and define the desired iteration pattern - through the trajectory - - run - The user facing run method. Calls the analysis methods - defined below - - Your analysis can implement the following methods, which are - called from run: - - _prepare - Called before iteration on the trajectory has begun. - Data structures can be set up at this time, however most - error checking should be done in the __init__ - - _single_frame - Called after the trajectory is moved onto each new frame. - - _conclude - Called once iteration on the trajectory is finished. - Apply normalisation and averaging to results here. + """Base class for defining multi frame analysis, it is designed as a + template for creating multiframe analysis. This class will automatically + take care of setting up the trajectory reader for iterating and offers to + show a progress meter. + + To define a new Analysis, `AnalysisBase` needs to be subclassed + `_single_frame` must be defined. It is also possible to define + `_prepare` and `_conclude` for pre and post processing. See the example + below. + + .. code-block:: python + class NewAnalysis(AnalysisBase): + def __init__(self, atomgroup, parameter, **kwargs): + super(NewAnalysis, self).__init__(atomgroup.universe.trajectory, + **kwargs) + self._parameter = parameter + self._ag = atomgroup + + def _prepare(self): + # OPTIONAL + # Called before iteration on the trajectory has begun. + # Data structures can be set up at this time + self.result = [] + + def _single_frame(self): + # REQUIRED + # Called after the trajectory is moved onto each new frame. + # store result of `some_function` for a single frame + self.result.append(some_function(self._ag, self._parameter)) + + def _conclude(self): + # OPTIONAL + # Called once iteration on the trajectory is finished. + # Apply normalisation and averaging to results here. + self.result = np.asarray(self.result) / np.sum(self.result) + + Afterwards the new analysis can be run like this. + + .. code-block:: python + na = NewAnalysis(u.select_atoms('name CA'), 35).run() + print(na.result) """ + def __init__(self, trajectory, start=None, + stop=None, step=None, quiet=True): + """ + Parameters + ---------- + trajectory : mda.Reader + A trajectory Reader + start : int, optional + start frame of analysis + stop : int, optional + stop frame of analysis + step : int, optional + number of frames to skip between each analysed frame + quiet : bool, optional + Turn off verbosity + """ + self._setup_frames(trajectory, start, stop, step) + interval = int(self.n_frames // 100) + if interval == 0: + interval = 1 + self._pm = ProgressMeter(self.nframes if self.nframes else 1, + interval=interval, quiet=quiet) + self._quiet = quiet + def _setup_frames(self, trajectory, start=None, stop=None, step=None): + """ + Pass a Reader object and define the desired iteration pattern + through the trajectory + + Parameters + ---------- + trajectory : mda.Reader + A trajectory Reader + start : int, optional + start frame of analysis + stop : int, optional + stop frame of analysis + step : int, optional + number of frames to skip between each analysed frame + """ self._trajectory = trajectory start, stop, step = trajectory.check_slice_indices( start, stop, step) @@ -78,7 +128,7 @@ def _single_frame(self): Don't worry about normalising, just deal with a single frame. """ - pass + raise NotImplementedError("Only implemented in child classes") def _prepare(self): """Set things up before the analysis loop begins""" @@ -91,7 +141,7 @@ def _conclude(self): """ pass - def run(self, **kwargs): + def run(self): """Perform the calculation""" logger.info("Starting preparation") self._prepare() @@ -101,5 +151,7 @@ def run(self, **kwargs): self._ts = ts # logger.info("--> Doing frame {} of {}".format(i+1, self.n_frames)) self._single_frame() + self._pm.echo(self._frame_index) logger.info("Finishing up") self._conclude() + return self diff --git a/testsuite/MDAnalysisTests/analysis/test_base.py b/testsuite/MDAnalysisTests/analysis/test_base.py index 038a2baacc4..94f26a86db1 100644 --- a/testsuite/MDAnalysisTests/analysis/test_base.py +++ b/testsuite/MDAnalysisTests/analysis/test_base.py @@ -2,8 +2,8 @@ # vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 fileencoding=utf-8 # # MDAnalysis --- http://www.MDAnalysis.org -# Copyright (c) 2006-2015 Naveen Michaud-Agrawal, Elizabeth J. Denning, Oliver Beckstein -# and contributors (see AUTHORS for the full list) +# Copyright (c) 2006-2015 Naveen Michaud-Agrawal, Elizabeth J. Denning, Oliver +# Beckstein and contributors (see AUTHORS for the full list) # # Released under the GNU Public Licence, v2 or any higher version # @@ -18,7 +18,7 @@ from six.moves import range from numpy.testing import ( - assert_, dec + assert_, dec, raises ) import MDAnalysis as mda @@ -30,19 +30,20 @@ class FrameAnalysis(AnalysisBase): """Just grabs frame numbers of frames it goes over""" - def __init__(self, reader, start=None, stop=None, step=None): + def __init__(self, reader, **kwargs): + super(FrameAnalysis, self).__init__(reader, **kwargs) self.traj = reader - self._setup_frames(reader, - start=start, - stop=stop, - step=step) - self.frames = [] def _single_frame(self): self.frames.append(self._ts.frame) +class IncompleteAnalysis(AnalysisBase): + def __init__(self, reader, **kwargs): + super(IncompleteAnalysis, self).__init__(reader, **kwargs) + + class TestAnalysisBase(object): @dec.skipif(parser_not_found('DCD'), 'DCD parser not available. Are you using python 3?') @@ -54,29 +55,29 @@ def tearDown(self): del self.u def test_default(self): - an = FrameAnalysis(self.u.trajectory) + an = FrameAnalysis(self.u.trajectory).run() assert_(an.n_frames == len(self.u.trajectory)) - - an.run() assert_(an.frames == list(range(len(self.u.trajectory)))) def test_start(self): - an = FrameAnalysis(self.u.trajectory, start=20) + an = FrameAnalysis(self.u.trajectory, start=20).run() assert_(an.n_frames == len(self.u.trajectory) - 20) - - an.run() assert_(an.frames == list(range(20, len(self.u.trajectory)))) def test_stop(self): - an = FrameAnalysis(self.u.trajectory, stop=20) + an = FrameAnalysis(self.u.trajectory, stop=20).run() assert_(an.n_frames == 20) - - an.run() assert_(an.frames == list(range(20))) def test_step(self): - an = FrameAnalysis(self.u.trajectory, step=20) + an = FrameAnalysis(self.u.trajectory, step=20).run() assert_(an.n_frames == 5) - - an.run() assert_(an.frames == list(range(98))[::20]) + + def test_quiet(self): + a = FrameAnalysis(self.u.trajectory, quiet=False) + assert_(not a._quiet) + + @raises(NotImplementedError) + def test_incomplete_defined_analysis(self): + IncompleteAnalysis(self.u.trajectory).run() From ba8c775a5d4a74e7125b48273d1de17c4d08aa0e Mon Sep 17 00:00:00 2001 From: Max Linke Date: Thu, 7 Jul 2016 20:47:33 +0200 Subject: [PATCH 03/10] Adopt AlignTrj --- package/MDAnalysis/analysis/align.py | 17 +++++++---------- .../MDAnalysisTests/analysis/test_align.py | 3 +-- 2 files changed, 8 insertions(+), 12 deletions(-) diff --git a/package/MDAnalysis/analysis/align.py b/package/MDAnalysis/analysis/align.py index 4b1c2c755a0..c7fe0945282 100644 --- a/package/MDAnalysis/analysis/align.py +++ b/package/MDAnalysis/analysis/align.py @@ -2,8 +2,8 @@ # vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 # # MDAnalysis --- http://www.MDAnalysis.org -# Copyright (c) 2006-2015 Naveen Michaud-Agrawal, Elizabeth J. Denning, Oliver Beckstein -# and contributors (see AUTHORS for the full list) +# Copyright (c) 2006-2015 Naveen Michaud-Agrawal, Elizabeth J. Denning, Oliver +# Beckstein and contributors (see AUTHORS for the full list) # # Released under the GNU Public Licence, v2 or any higher version # @@ -472,8 +472,7 @@ class AlignTraj(AnalysisBase): def __init__(self, mobile, reference, select='all', filename=None, prefix='rmsfit_', mass_weighted=False, tol_mass=0.1, - strict=False, force=True, quiet=False, start=None, stop=None, - step=None, **kwargs): + strict=False, force=True, **kwargs): """Initialization Parameters @@ -516,18 +515,16 @@ def __init__(self, mobile, reference, select='all', filename=None, exception """ - self._quiet = quiet + super(AlignTraj, self).__init__(mobile.trajectory, **kwargs) if self._quiet: logging.disable(logging.WARN) - traj = mobile.trajectory select = rms.process_selection(select) self.ref_atoms = reference.select_atoms(*select['reference']) self.mobile_atoms = mobile.select_atoms(*select['mobile']) - kwargs.setdefault('remarks', 'RMS fitted trajectory to reference') if filename is None: - path, fn = os.path.split(traj.filename) + path, fn = os.path.split(self._trajectory.filename) filename = os.path.join(path, prefix + fn) logger.info('filename of rms_align with no filename given' ': {0}'.format(filename)) @@ -540,7 +537,8 @@ def __init__(self, mobile, reference, select='all', filename=None, natoms = self.mobile_atoms.n_atoms self.ref_atoms, self.mobile_atoms = get_matching_atoms( - self.ref_atoms, self.mobile_atoms, tol_mass=tol_mass, strict=strict) + self.ref_atoms, self.mobile_atoms, tol_mass=tol_mass, + strict=strict) self._writer = mda.Writer(self.filename, natoms) @@ -551,7 +549,6 @@ def __init__(self, mobile, reference, select='all', filename=None, self._weights = None logger.info("RMS-fitting on {0:d} atoms.".format(len(self.ref_atoms))) - self._setup_frames(traj, start, stop, step) def _prepare(self): # reference centre of mass system diff --git a/testsuite/MDAnalysisTests/analysis/test_align.py b/testsuite/MDAnalysisTests/analysis/test_align.py index dddccae33ab..3813cf4df57 100644 --- a/testsuite/MDAnalysisTests/analysis/test_align.py +++ b/testsuite/MDAnalysisTests/analysis/test_align.py @@ -202,8 +202,7 @@ def different_atoms(): assert_raises(SelectionError, different_atoms) - -class TestAlignmentProcessing(TestCase): +class TestAlignmentProcessing(object): def setUp(self): self.seq = FASTA self.tempdir = tempdir.TempDir() From 174042c129ff59eddbb223291048610ebe7db7e4 Mon Sep 17 00:00:00 2001 From: Max Linke Date: Thu, 7 Jul 2016 20:30:15 +0200 Subject: [PATCH 04/10] Adopt Contacts --- package/MDAnalysis/analysis/contacts.py | 9 ++++----- testsuite/MDAnalysisTests/analysis/test_contacts.py | 6 ++---- 2 files changed, 6 insertions(+), 9 deletions(-) diff --git a/package/MDAnalysis/analysis/contacts.py b/package/MDAnalysis/analysis/contacts.py index 2f958a92136..ede234f5b93 100644 --- a/package/MDAnalysis/analysis/contacts.py +++ b/package/MDAnalysis/analysis/contacts.py @@ -376,7 +376,7 @@ class Contacts(AnalysisBase): """ def __init__(self, u, selection, refgroup, method="hard_cut", radius=4.5, - kwargs=None, start=None, stop=None, step=None,): + kwargs=None, **basekwargs): """Initialization Parameters @@ -404,6 +404,9 @@ def __init__(self, u, selection, refgroup, method="hard_cut", radius=4.5, Step between frames to analyse, Default: 1 """ + self.u = u + super(Contacts, self).__init__(self.u.trajectory, **basekwargs) + if method == 'hard_cut': self.fraction_contacts = hard_cut_q elif method == 'soft_cut': @@ -413,10 +416,6 @@ def __init__(self, u, selection, refgroup, method="hard_cut", radius=4.5, raise ValueError("method has to be callable") self.fraction_contacts = method - # setup boilerplate - self.u = u - self._setup_frames(self.u.trajectory, start, stop, step) - self.selection = selection self.grA = u.select_atoms(selection[0]) self.grB = u.select_atoms(selection[1]) diff --git a/testsuite/MDAnalysisTests/analysis/test_contacts.py b/testsuite/MDAnalysisTests/analysis/test_contacts.py index 712b5cf6ac4..dbcdad587b7 100644 --- a/testsuite/MDAnalysisTests/analysis/test_contacts.py +++ b/testsuite/MDAnalysisTests/analysis/test_contacts.py @@ -162,14 +162,12 @@ def tearDown(self): def _run_Contacts(self, **kwargs): acidic = self.universe.select_atoms(self.sel_acidic) basic = self.universe.select_atoms(self.sel_basic) - Contacts = contacts.Contacts( + return contacts.Contacts( self.universe, selection=(self.sel_acidic, self.sel_basic), refgroup=(acidic, basic), radius=6.0, - **kwargs) - Contacts.run() - return Contacts + **kwargs).run() def test_startframe(self): """test_startframe: TestContactAnalysis1: start frame set to 0 (resolution of From 4c9017126d13bc2874bc9e34c39e74a095854a8e Mon Sep 17 00:00:00 2001 From: Max Linke Date: Thu, 7 Jul 2016 20:36:26 +0200 Subject: [PATCH 05/10] Addopt Persistencelength --- package/MDAnalysis/analysis/polymer.py | 18 +++++++----------- 1 file changed, 7 insertions(+), 11 deletions(-) diff --git a/package/MDAnalysis/analysis/polymer.py b/package/MDAnalysis/analysis/polymer.py index 9e54492b114..754c6c37c89 100644 --- a/package/MDAnalysis/analysis/polymer.py +++ b/package/MDAnalysis/analysis/polymer.py @@ -2,8 +2,8 @@ # vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 # # MDAnalysis --- http://www.MDAnalysis.org -# Copyright (c) 2006-2015 Naveen Michaud-Agrawal, Elizabeth J. Denning, Oliver Beckstein -# and contributors (see AUTHORS for the full list) +# Copyright (c) 2006-2015 Naveen Michaud-Agrawal, Elizabeth J. Denning, Oliver +# Beckstein and contributors (see AUTHORS for the full list) # # Released under the GNU Public Licence, v2 or any higher version # @@ -32,7 +32,6 @@ import logging from .. import NoDataError -from ..lib.util import blocks_of from ..lib.distances import calc_bonds from .base import AnalysisBase @@ -51,8 +50,7 @@ class PersistenceLength(AnalysisBase): .. versionadded:: 0.13.0 """ - def __init__(self, atomgroups, - start=None, stop=None, step=None): + def __init__(self, atomgroups, **kwargs): """Calculate the persistence length for polymer chains Parameters @@ -67,17 +65,16 @@ def __init__(self, atomgroups, step : int, optional Step between frames to analyse, Default: 1 """ + super(PersistenceLength, self).__init__( + atomgroups[0].universe.trajectory, **kwargs) self._atomgroups = atomgroups # Check that all chains are the same length lens = [len(ag) for ag in atomgroups] chainlength = len(atomgroups[0]) - if not all( l == chainlength for l in lens): + if not all(l == chainlength for l in lens): raise ValueError("Not all AtomGroups were the same size") - self._setup_frames(atomgroups[0].universe.trajectory, - start, stop, step) - self._results = np.zeros(chainlength - 1, dtype=np.float32) def _single_frame(self): @@ -117,10 +114,9 @@ def _calc_bond_length(self): def perform_fit(self): """Fit the results to an exponential decay""" - from scipy.optimize import curve_fit try: - results = self.results + self.results except AttributeError: raise NoDataError("Use the run method first") self.x = np.arange(len(self.results)) * self.lb From d55a4a0a0cce1c196f8d631ebb01341b89091533 Mon Sep 17 00:00:00 2001 From: Max Linke Date: Thu, 7 Jul 2016 20:39:02 +0200 Subject: [PATCH 06/10] Adopt InterRDF --- package/MDAnalysis/analysis/rdf.py | 16 ++++++---------- 1 file changed, 6 insertions(+), 10 deletions(-) diff --git a/package/MDAnalysis/analysis/rdf.py b/package/MDAnalysis/analysis/rdf.py index a4b54de230e..48f8a225aa2 100644 --- a/package/MDAnalysis/analysis/rdf.py +++ b/package/MDAnalysis/analysis/rdf.py @@ -2,8 +2,8 @@ # vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 # # MDAnalysis --- http://www.MDAnalysis.org -# Copyright (c) 2006-2015 Naveen Michaud-Agrawal, Elizabeth J. Denning, Oliver Beckstein -# and contributors (see AUTHORS for the full list) +# Copyright (c) 2006-2015 Naveen Michaud-Agrawal, Elizabeth J. Denning, Oliver +# Beckstein and contributors (see AUTHORS for the full list) # # Released under the GNU Public Licence, v2 or any higher version # @@ -79,18 +79,14 @@ class InterRDF(AnalysisBase): """ def __init__(self, g1, g2, nbins=75, range=(0.0, 15.0), exclusion_block=None, - start=None, stop=None, step=None): + **kwargs): + super(InterRDF, self).__init__(g1.universe.trajectory, **kwargs) self.g1 = g1 self.g2 = g2 self.u = g1.universe - self._setup_frames(self.u.trajectory, - start=start, - stop=stop, - step=step) - - self.rdf_settings = {'bins':nbins, - 'range':range} + self.rdf_settings = {'bins': nbins, + 'range': range} # Empty histogram to store the RDF count, edges = np.histogram([-1], **self.rdf_settings) From ba662ae5386f4996076eae7f9759a2a4428f6306 Mon Sep 17 00:00:00 2001 From: Max Linke Date: Thu, 7 Jul 2016 20:43:47 +0200 Subject: [PATCH 07/10] Adopt LinearDensity --- package/MDAnalysis/analysis/lineardensity.py | 91 +++++++++++-------- .../analysis/test_lineardensity.py | 16 ++-- 2 files changed, 59 insertions(+), 48 deletions(-) diff --git a/package/MDAnalysis/analysis/lineardensity.py b/package/MDAnalysis/analysis/lineardensity.py index 69df9fbc28b..d0c54bd784b 100644 --- a/package/MDAnalysis/analysis/lineardensity.py +++ b/package/MDAnalysis/analysis/lineardensity.py @@ -2,8 +2,8 @@ # vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 # # MDAnalysis --- http://www.MDAnalysis.org -# Copyright (c) 2006-2015 Naveen Michaud-Agrawal, Elizabeth J. Denning, Oliver Beckstein -# and contributors (see AUTHORS for the full list) +# Copyright (c) 2006-2015 Naveen Michaud-Agrawal, Elizabeth J. Denning, Oliver +# Beckstein and contributors (see AUTHORS for the full list) # # Released under the GNU Public Licence, v2 or any higher version # @@ -13,7 +13,6 @@ # MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations. # J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787 # - """ Linear Density --- :mod:`MDAnalysis.analysis.lineardensity` =========================================================== @@ -29,6 +28,7 @@ from MDAnalysis.analysis.base import AnalysisBase + class LinearDensity(AnalysisBase): """Linear density profile LinearDensity(selection, grouping='atoms', binsize=0.25) @@ -71,11 +71,13 @@ class LinearDensity(AnalysisBase): .. versionadded:: 0.14.0 """ - def __init__(self, selection, grouping='atoms', binsize=0.25, - start=None, stop=None, step=None): - self._ags = [selection] # allows use of run(parallel=True) + + def __init__(self, selection, grouping='atoms', binsize=0.25, **kwargs): + super(LinearDensity, self).__init__(selection.universe.trajectory, + **kwargs) + # allows use of run(parallel=True) + self._ags = [selection] self._universe = selection.universe - self._setup_frames(self._universe.trajectory, start, stop, step) self.binsize = binsize @@ -88,12 +90,13 @@ def __init__(self, selection, grouping='atoms', binsize=0.25, # Box sides self.dimensions = self._universe.dimensions[:3] self.volume = np.prod(self.dimensions) - bins = (self.dimensions // self.binsize).astype(int) # number of bins + # number of bins + bins = (self.dimensions // self.binsize).astype(int) # Here we choose a number of bins of the largest cell side so that # x, y and z values can use the same "coord" column in the output file self.nbins = bins.max() - slices_vol = self.volume/bins + slices_vol = self.volume / bins self.keys = ['pos', 'pos_std', 'char', 'char_std'] @@ -115,10 +118,10 @@ def _prepare(self): group = getattr(self._ags[0], self.grouping) # Get masses and charges for the selection - try: # in case it's not an atom + try: # in case it's not an atom self.masses = np.array([elem.total_mass() for elem in group]) self.charges = np.array([elem.total_charge() for elem in group]) - except AttributeError: # much much faster for atoms + except AttributeError: # much much faster for atoms self.masses = self._ags[0].masses self.charges = self._ags[0].charges @@ -130,9 +133,10 @@ def _single_frame(self): # Find position of atom/group of atoms if self.grouping == 'atoms': - positions = self._ags[0].positions # faster for atoms + positions = self._ags[0].positions # faster for atoms else: - positions = np.array([elem.centroid() for elem in self.group]) # COM for res/frag/etc + # COM for res/frag/etc + positions = np.array([elem.centroid() for elem in self.group]) for dim in ['x', 'y', 'z']: idx = self.results[dim]['dim'] @@ -140,8 +144,10 @@ def _single_frame(self): key = 'pos' key_std = 'pos_std' # histogram for positions weighted on masses - hist, _ = np.histogram(positions[:, idx], weights=self.masses, - bins=self.nbins, range=(0.0, max(self.dimensions))) + hist, _ = np.histogram(positions[:, idx], + weights=self.masses, + bins=self.nbins, + range=(0.0, max(self.dimensions))) self.results[dim][key] += hist self.results[dim][key_std] += np.square(hist) @@ -149,33 +155,35 @@ def _single_frame(self): key = 'char' key_std = 'char_std' # histogram for positions weighted on charges - hist, _ = np.histogram(positions[:, idx], weights=self.charges, - bins=self.nbins, range=(0.0, max(self.dimensions))) + hist, _ = np.histogram(positions[:, idx], + weights=self.charges, + bins=self.nbins, + range=(0.0, max(self.dimensions))) self.results[dim][key] += hist self.results[dim][key_std] += np.square(hist) def _conclude(self): - k = 6.022e-1 # divide by avodagro and convert from A3 to cm3 + k = 6.022e-1 # divide by avodagro and convert from A3 to cm3 # Average results over the number of configurations for dim in ['x', 'y', 'z']: for key in ['pos', 'pos_std', 'char', 'char_std']: self.results[dim][key] /= self.n_frames # Compute standard deviation for the error - self.results[dim]['pos_std'] = np.sqrt(self.results[dim]['pos_std'] - - np.square(self.results[dim]['pos'])) - self.results[dim]['char_std'] = np.sqrt(self.results[dim]['char_std'] - - np.square(self.results[dim]['char'])) + self.results[dim]['pos_std'] = np.sqrt(self.results[dim][ + 'pos_std'] - np.square(self.results[dim]['pos'])) + self.results[dim]['char_std'] = np.sqrt(self.results[dim][ + 'char_std'] - np.square(self.results[dim]['char'])) for dim in ['x', 'y', 'z']: - self.results[dim]['pos'] /= (self.results[dim]['slice volume']*k) - self.results[dim]['char'] /= (self.results[dim]['slice volume']*k) - self.results[dim]['pos_std'] /= (self.results[dim]['slice volume']*k) - self.results[dim]['char_std'] /= (self.results[dim]['slice volume']*k) + self.results[dim]['pos'] /= self.results[dim]['slice volume'] * k + self.results[dim]['char'] /= self.results[dim]['slice volume'] * k + self.results[dim]['pos_std'] /= self.results[dim]['slice volume'] * k + self.results[dim]['char_std'] /= self.results[dim]['slice volume'] * k def save(self, description='', form='txt'): - """ Save density profile to file + """Save density profile to file Allows to save the density profile to either a ASCII txt file or a binary numpy npz file. Output file has extension 'ldens' and begins with the name of the trajectory file. @@ -185,22 +193,24 @@ def save(self, description='', form='txt'): description : str An arbitrary description added to the output filename. Can be useful form : str {'txt', 'npz'} - Format of the output. 'txt' will generate an ASCII text file while 'npz' - will produce a numpy binary file. + Format of the output. 'txt' will generate an ASCII text file while + 'npz' will produce a numpy binary file. Example ------- - After initializing and running a `LinearDensity` object, results can be written - to file as follows: + After initializing and running a `LinearDensity` object, results can be + written to file as follows: ldens.save(description='mydensprof', form='txt') which will output the linear density profiles in a file named .mydensprof_.ldens + """ # Take root of trajectory filename for output file naming - trajname = path.splitext(path.basename(self._universe.trajectory.filename))[0] + trajname = path.splitext(path.basename( + self._universe.trajectory.filename))[0] # additional string for naming the output file description = description + "_" + str(self.grouping) - filename = trajname+"." + description + ".ldens" + filename = trajname + "." + description + ".ldens" if form is 'txt': self._savetxt(filename) @@ -223,11 +233,14 @@ def _savetxt(self, filename): output.append(self.results[dim]['char']) output.append(self.results[dim]['char_std']) - density = self.totalmass/self.volume - header = "1 coord [Ang] 2-7 mass density (x,sx,y,sz,z,sz) [g/cm^3]" + \ - "8-13 charge density (x,sx,y,sz,z,sz) [e/A^3]\n Average density: "\ - + str(density) + " g/cm3" - np.savetxt(filename, np.column_stack(output), fmt='%10.5f', header=header) + density = self.totalmass / self.volume + header = ("1 coord [Ang] 2-7 mass density (x,sx,y,sz,z,sz) [g/cm^3]" + "8-13 charge density (x,sx,y,sz,z,sz) [e/A^3]\n Average " + "density: {} g/cm3".format(density)) + np.savetxt(filename, + np.column_stack(output), + fmt='%10.5f', + header=header) def _savez(self, filename): bins = np.linspace(0.0, max(self.dimensions), num=self.nbins) @@ -237,7 +250,7 @@ def _savez(self, filename): self.results[dim].pop('dim') self.results[dim].pop('slice volume') for key in self.results[dim]: - dictionary[dim+"_"+key] = self.results[dim][key] + dictionary[dim + "_" + key] = self.results[dim][key] np.savez(filename, **dictionary) diff --git a/testsuite/MDAnalysisTests/analysis/test_lineardensity.py b/testsuite/MDAnalysisTests/analysis/test_lineardensity.py index 6c861518d95..42ed115038a 100644 --- a/testsuite/MDAnalysisTests/analysis/test_lineardensity.py +++ b/testsuite/MDAnalysisTests/analysis/test_lineardensity.py @@ -2,8 +2,8 @@ # vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 fileencoding=utf-8 # # MDAnalysis --- http://www.MDAnalysis.org -# Copyright (c) 2006-2015 Naveen Michaud-Agrawal, Elizabeth J. Denning, Oliver Beckstein -# and contributors (see AUTHORS for the full list) +# Copyright (c) 2006-2015 Naveen Michaud-Agrawal, Elizabeth J. Denning, Oliver +# Beckstein and contributors (see AUTHORS for the full list) # # Released under the GNU Public Licence, v2 or any higher version # @@ -20,23 +20,21 @@ from MDAnalysis.analysis.lineardensity import LinearDensity from numpy.testing import TestCase, assert_allclose + class TestLinearDensity(TestCase): def setUp(self): self.universe = mda.Universe(waterPSF, waterDCD) self.sel_string = 'all' self.selection = self.universe.select_atoms(self.sel_string) - - self.xpos = np.array([ 0., 0., 0., 0.0072334, 0.00473299, 0., - 0., 0., 0., 0. ]) + + self.xpos = np.array([0., 0., 0., 0.0072334, 0.00473299, 0., + 0., 0., 0., 0.]) def test_serial(self): - ld = LinearDensity(self.selection, binsize=5) - ld.run() + ld = LinearDensity(self.selection, binsize=5).run() assert_allclose(self.xpos, ld.results['x']['pos'], rtol=1e-6, atol=0) # def test_parallel(self): # ld = LinearDensity(self.universe, self.selection, binsize=5) # ld.run(parallel=True) # assert_equal(self.xpos, ld.results['x']['pos']) - - From 9755a9ffe875b85a4f35fb3576b559592eb39f8b Mon Sep 17 00:00:00 2001 From: Max Linke Date: Thu, 7 Jul 2016 20:34:56 +0200 Subject: [PATCH 08/10] Follow best practices for matplotlib Returning an axis object is considered to be better since the user is able to adjust the plot by himself if he wishes to do so. --- package/MDAnalysis/analysis/polymer.py | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/package/MDAnalysis/analysis/polymer.py b/package/MDAnalysis/analysis/polymer.py index 754c6c37c89..53e3044f40d 100644 --- a/package/MDAnalysis/analysis/polymer.py +++ b/package/MDAnalysis/analysis/polymer.py @@ -125,15 +125,16 @@ def perform_fit(self): self.fit = np.exp(-self.x/self.lp) - def plot(self): + def plot(self, ax=None): """Oooh fancy""" import matplotlib.pyplot as plt - plt.ylabel('C(x)') - plt.xlabel('x') - plt.xlim([0.0, 40 * self.lb]) - plt.plot(self.x, self.results, 'ro') - plt.plot(self.x, self.fit) - plt.show() + if ax is None: + ax = plt.gca() + ax.plot(self.x, self.results, 'ro', label='Result') + ax.plot(self.x, self.fit, label='Fit') + ax.set(xlabel='x', ylabel='C(x)', xlim=[0.0, 40 * self.lb]) + ax.legend(loc='best') + return ax def fit_exponential_decay(x, y): From 9edaa24923d653bd70260441ddf5c83b30dfdbe4 Mon Sep 17 00:00:00 2001 From: Max Linke Date: Fri, 8 Jul 2016 08:02:51 +0200 Subject: [PATCH 09/10] update changelog --- package/CHANGELOG | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/package/CHANGELOG b/package/CHANGELOG index c4966abaf02..b6433cfce90 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -18,6 +18,13 @@ The rules for this file: * 0.15.1 +Enhancements + + * All classes derived from 'AnalysisBase' now display a ProgressMeter + with 'quiet=False' + * The 'run' method from all 'AnalysisBase' derived classes return the + class itself. + Fixes * GROWriter resids now truncated properly (Issue #886) * reading/writing lambda value in trr files (Issue #859) From ffdf47931e98e5db74d4b33c4c09b99a3980a0e4 Mon Sep 17 00:00:00 2001 From: Max Linke Date: Fri, 8 Jul 2016 10:23:29 +0200 Subject: [PATCH 10/10] Ensure old API still works Since user code might already exist we keep make sure the old API still works as well. --- package/MDAnalysis/analysis/base.py | 17 +++++++++++------ testsuite/MDAnalysisTests/analysis/test_base.py | 12 ++++++++++++ 2 files changed, 23 insertions(+), 6 deletions(-) diff --git a/package/MDAnalysis/analysis/base.py b/package/MDAnalysis/analysis/base.py index cf99e2bb669..3d7d5efe79f 100644 --- a/package/MDAnalysis/analysis/base.py +++ b/package/MDAnalysis/analysis/base.py @@ -90,13 +90,8 @@ def __init__(self, trajectory, start=None, quiet : bool, optional Turn off verbosity """ - self._setup_frames(trajectory, start, stop, step) - interval = int(self.n_frames // 100) - if interval == 0: - interval = 1 - self._pm = ProgressMeter(self.nframes if self.nframes else 1, - interval=interval, quiet=quiet) self._quiet = quiet + self._setup_frames(trajectory, start, stop, step) def _setup_frames(self, trajectory, start=None, stop=None, step=None): @@ -122,6 +117,16 @@ def _setup_frames(self, trajectory, start=None, self.stop = stop self.step = step self.n_frames = len(range(start, stop, step)) + interval = int(self.n_frames // 100) + if interval == 0: + interval = 1 + + # ensure _quiet is set when __init__ wasn't called, this is to not + # break pre 0.16.0 API usage of AnalysisBase + if not hasattr(self, '_quiet'): + self._quiet = True + self._pm = ProgressMeter(self.n_frames if self.n_frames else 1, + interval=interval, quiet=self._quiet) def _single_frame(self): """Calculate data from a single frame of trajectory diff --git a/testsuite/MDAnalysisTests/analysis/test_base.py b/testsuite/MDAnalysisTests/analysis/test_base.py index 94f26a86db1..ee1a1e3d5c8 100644 --- a/testsuite/MDAnalysisTests/analysis/test_base.py +++ b/testsuite/MDAnalysisTests/analysis/test_base.py @@ -44,6 +44,15 @@ def __init__(self, reader, **kwargs): super(IncompleteAnalysis, self).__init__(reader, **kwargs) +class OldAPIAnalysis(AnalysisBase): + """for version 0.15.0""" + def __init__(self, reader, **kwargs): + self._setup_frames(reader, **kwargs) + + def _single_frame(self): + pass + + class TestAnalysisBase(object): @dec.skipif(parser_not_found('DCD'), 'DCD parser not available. Are you using python 3?') @@ -81,3 +90,6 @@ def test_quiet(self): @raises(NotImplementedError) def test_incomplete_defined_analysis(self): IncompleteAnalysis(self.u.trajectory).run() + + def test_old_api(self): + OldAPIAnalysis(self.u.trajectory).run()