diff --git a/package/CHANGELOG b/package/CHANGELOG index f0f791f0da5..22e6ada09e6 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -106,6 +106,7 @@ Fixes * Fix syntax warning over comparison of literals using is (Issue #3066) Enhancements + * Added `Results` class for storing analysis results (#3115, PR #3233) * Added intra_bonds, intra_angles, intra_dihedrals etc. to return only the connections involving atoms within the AtomGroup, instead of including atoms outside the AtomGroup (Issue #1264, #2821, PR #3200) @@ -173,6 +174,8 @@ Enhancements checking if it can be used in parallel analysis. (Issue #2996, PR #2950) Changes + * `GNMAnalysis`, `LinearDensity`, `PersistenceLength` and + `AnalysisFromFunction` use the `results` attribute. * Fixed inaccurate docstring inside the RMSD class (Issue #2796, PR #3134) * TPRParser now loads TPR files with `tpr_resid_from_one=True` by default, which starts TPR resid indexing from 1 (instead of 0 as in 1.x) (Issue #2364, PR #3152) diff --git a/package/MDAnalysis/analysis/__init__.py b/package/MDAnalysis/analysis/__init__.py index 1ac4648978a..1c9f78b1880 100644 --- a/package/MDAnalysis/analysis/__init__.py +++ b/package/MDAnalysis/analysis/__init__.py @@ -21,102 +21,6 @@ # J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787 # -""" -:mod:`MDAnalysis.analysis` --- Analysis code based on MDAnalysis -================================================================ - -The :mod:`MDAnalysis.analysis` sub-package contains various recipes and -algorithms that can be used to analyze MD trajectories. - -If you use them please check if the documentation mentions any specific caveats -and also if there are any published papers associated with these algorithms. - -Available analysis modules --------------------------- - -:mod:`~MDAnalysis.analysis.align` - Fitting and aligning of coordinate frames, including the option to - use a sequence alignment to define equivalent atoms to fit on. - -:mod:`~MDAnalysis.analysis.contacts` - Analyse the number of native contacts relative to a reference - state, also known as a "q1-q2" analysis. - -:mod:`~MDAnalysis.analysis.density` - Creating and manipulating densities such as the density ow water - molecules around a protein. Makes use of the external - GridDataFormats_ package. - -:mod:`~MDAnalysis.analysis.distances` - Functions to calculate distances between atoms and selections; it - contains the often-used - :func:`~MDAnalysis.analysis.distances.distance_array` function. - -:mod:`~MDAnalysis.analysis.hbonds` - Analyze hydrogen bonds, including both the per frame results as well - as the dynamic properties and lifetimes. - -:mod:`~MDAnalysis.analysis.helix_analysis` - Analysis of helices with the HELANAL_ algorithm. - -:mod:`~MDAnalysis.analysis.hole2` - Run and process output from the :program:`HOLE` program - to analyze pores, tunnels and cavities in proteins. - -:mod:`~MDAnalysis.analysis.gnm` - Gaussian normal mode analysis of MD trajectories with the - help of an elastic network. - -:mod:`~MDAnalysis.analysis.leaflet` - Find lipids in the upper and lower (or inner and outer) leaflet of - a bilayer; the algorithm can deal with any deformations as long as - the two leaflets are topologically distinct. - -:mod:`~MDAnalysis.analysis.msd` - Implements the calculation of Mean Squared Displacements (MSDs) by the - Einstein relation. - -:mod:`~MDAnalysis.analysis.nuclinfo` - Analyse the nucleic acid for the backbone dihedrals, chi, sugar - pucker, and Watson-Crick distance (minor and major groove - distances). - -:mod:`~MDAnalysis.analysis.psa` - Perform Path Similarity Analysis (PSA) on a set of trajectories to measure - their mutual similarities, including the ability to perform hierarchical - clustering and generate heat map-dendrogram plots. - -:mod:`~MDAnalysis.analysis.rdf` - Calculation of pair distribution functions - -:mod:`~MDAnalysis.analysis.rms` - Calculation of RMSD and RMSF. - -:mod:`~MDAnalysis.analysis.waterdynamics` - Analysis of water. - -:mod:`~MDAnalysis.analysis.legacy.x3dna` - Analysis of helicoidal parameters driven by X3DNA_. (Note that this - module is not fully supported any more and needs to be explicitly - imported from :mod:`MDAnalysis.analysis.legacy`.) - -.. _GridDataFormats: https://github.com/orbeckst/GridDataFormats -.. _HELANAL: http://www.ccrnp.ncifcrf.gov/users/kumarsan/HELANAL/helanal.html -.. _X3DNA: http://x3dna.org/ - -.. versionchanged:: 0.10.0 - The analysis submodules are not automatically imported any more. Manually - import any submodule that you need. - -.. versionchanged:: 0.16.0 - :mod:`~MDAnalysis.analysis.legacy.x3dna` was moved to the - :mod:`MDAnalysis.analysis.legacy` package - -.. versionchanged:: 2.0.0 - Adds MSD, and changes hole for hole2.hole. - -""" - __all__ = [ 'align', 'base', diff --git a/package/MDAnalysis/analysis/base.py b/package/MDAnalysis/analysis/base.py index 7237fcae26e..c7619203581 100644 --- a/package/MDAnalysis/analysis/base.py +++ b/package/MDAnalysis/analysis/base.py @@ -24,10 +24,10 @@ Analysis building blocks --- :mod:`MDAnalysis.analysis.base` ============================================================ -A collection of useful building blocks for creating Analysis -classes. +The building blocks for creating Analysis classes. """ +from collections import UserDict import inspect import logging import itertools @@ -40,20 +40,127 @@ logger = logging.getLogger(__name__) +class Results(UserDict): + r"""Container object for storing results. + + ``Results`` are dictionaries that provide two ways by which values can be + accessed: by dictionary key ``results["value_key"]`` or by object + attribute, ``results.value_key``. ``Results`` stores all results obtained + from an analysis after calling :func:`run()`. + + The implementation is similar to the :class:`sklearn.utils.Bunch` + class in `scikit-learn`_. + + .. _`scikit-learn`: https://scikit-learn.org/ + + Raises + ------ + AttributeError + If an assigned attribute has the same name as a default attribute. + + ValueError + If a key is not of type ``str`` and therefore is not able to be + accessed by attribute. + + Notes + ----- + Pickling of ``Results`` is currently not supported. + + Examples + -------- + >>> from MDAnalysis.analysis.base import Results + >>> results = Results(a=1, b=2) + >>> results['b'] + 2 + >>> results.b + 2 + >>> results.a = 3 + >>> results['a'] + 3 + >>> results.c = [1, 2, 3, 4] + >>> results['c'] + [1, 2, 3, 4] + + + .. versionadded:: 2.0.0 + """ + def _validate_key(self, key): + if key in dir(UserDict) or (key == "data" and self._dict_frozen): + raise AttributeError(f"'{key}' is a protected dictionary " + "attribute") + elif isinstance(key, str) and not key.isidentifier(): + raise ValueError(f"'{key}' is not a valid attribute") + + def __init__(self, **kwargs): + if "data" in kwargs.keys(): + raise AttributeError(f"'data' is a protected dictionary attribute") + + self._dict_frozen = False + for key in kwargs: + self._validate_key(key) + super().__init__(**kwargs) + self._dict_frozen = True + + def __setitem__(self, key, item): + self._validate_key(key) + super().__setitem__(key, item) + + def __setattr__(self, attr, value): + self._validate_key(attr) + super().__setattr__(attr, value) + + # attribute available as key + if self._dict_frozen and attr != "_dict_frozen": + super().__setitem__(attr, value) + + def __getattr__(self, attr): + try: + return self.data[attr] + except KeyError as err: + raise AttributeError("'Results' object has no " + f"attribute '{attr}'") from err + + class AnalysisBase(object): - """Base class for defining multi frame analysis + r"""Base class for defining multi-frame analysis - The class it is designed as a template for creating multiframe analyses. + The class is designed as a template for creating multi-frame analyses. This class will automatically take care of setting up the trajectory reader for iterating, and it offers to show a progress meter. + Computed results are stored inside the :attr:`results` attribute. 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. + and :meth:`_single_frame` must be defined. It is also possible to define + :meth:`_prepare` and :meth:`_conclude` for pre- and post-processing. + All results should be stored as attributes of the :class:`Results` + container. + Parameters + ---------- + trajectory : MDAnalysis.coordinates.base.ReaderBase + A trajectory Reader + verbose : bool, optional + Turn on more logging and debugging + + Attributes + ---------- + times: numpy.ndarray + array of Timestep times. Only exists after calling + :meth:`AnalysisBase.run` + frames: numpy.ndarray + array of Timestep frame indices. Only exists after calling + :meth:`AnalysisBase.run` + results: :class:`Results` + results of calculation are stored after call + to :meth:`AnalysisBase.run` + + + Example + ------- .. code-block:: python + from MDAnalysis.analysis.base import AnalysisBase + class NewAnalysis(AnalysisBase): def __init__(self, atomgroup, parameter, **kwargs): super(NewAnalysis, self).__init__(atomgroup.universe.trajectory, @@ -65,53 +172,51 @@ def _prepare(self): # OPTIONAL # Called before iteration on the trajectory has begun. # Data structures can be set up at this time - self.result = [] + self.results.example_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)) + # store an example_result of `some_function` for a single frame + self.results.example_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) + self.results.example_result = np.asarray(self.example_result) + self.results.example_result /= np.sum(self.result) - Afterwards the new analysis can be run like this. + Afterwards the new analysis can be run like this .. code-block:: python - na = NewAnalysis(u.select_atoms('name CA'), 35).run(start=10, stop=20) - print(na.result) + import MDAnalysis as mda + from MDAnalysisTests.datafiles import PSF, DCD - Attributes - ---------- - times: np.ndarray - array of Timestep times. Only exists after calling run() - frames: np.ndarray - array of Timestep frame indices. Only exists after calling run() + u = mda.Universe(PSF, DCD) - """ + na = NewAnalysis(u.select_atoms('name CA'), 35) + na.run(start=10, stop=20) + print(na.results.example_result) + # results can also be accessed by key + print(na.results["example_result"]) - def __init__(self, trajectory, verbose=False, **kwargs): - """ - Parameters - ---------- - trajectory : mda.Reader - A trajectory Reader - verbose : bool, optional - Turn on more logging and debugging, default ``False`` + .. versionchanged:: 1.0.0 + Support for setting ``start``, ``stop``, and ``step`` has been + removed. These should now be directly passed to + :meth:`AnalysisBase.run`. - .. versionchanged:: 1.0.0 - Support for setting ``start``, ``stop``, and ``step`` has been - removed. These should now be directly passed to - :meth:`AnalysisBase.run`. - """ + .. versionchanged:: 2.0.0 + Added :attr:`results` + """ + + def __init__(self, trajectory, verbose=False, **kwargs): self._trajectory = trajectory self._verbose = verbose + self.results = Results() def _setup_frames(self, trajectory, start=None, stop=None, step=None): """ @@ -132,7 +237,6 @@ def _setup_frames(self, trajectory, start=None, stop=None, step=None): .. versionchanged:: 1.0.0 Added .frames and .times arrays as attributes - """ self._trajectory = trajectory start, stop, step = trajectory.check_slice_indices(start, stop, step) @@ -157,7 +261,7 @@ def _prepare(self): def _conclude(self): """Finalise the results you've gathered. - Called at the end of the run() method to finish everything up. + Called at the end of the :meth:`run` method to finish everything up. """ pass # pylint: disable=unnecessary-pass @@ -198,46 +302,57 @@ def run(self, start=None, stop=None, step=None, verbose=None): class AnalysisFromFunction(AnalysisBase): - """ - Create an analysis from a function working on AtomGroups + r"""Create an :class:`AnalysisBase` from a function working on AtomGroups + + Parameters + ---------- + function : callable + function to evaluate at each frame + trajectory : mda.coordinates.Reader, optional + trajectory to iterate over. If ``None`` the first AtomGroup found in + args and kwargs is used as a source for the trajectory. + *args : list + arguments for ``function`` + **kwargs : dict + arguments for ``function`` and :class:`AnalysisBase` Attributes ---------- - results : ndarray - results of calculation are stored after call to ``run`` + results.frames : numpy.ndarray + simulation frames used in analysis + results.times : numpy.ndarray + simulation times used in analysis + results.timeseries : numpy.ndarray + Results for each frame of the wrapped function, + stored after call to :meth:`AnalysisFromFunction.run`. + + Raises + ------ + ValueError + if ``function`` has the same ``kwargs`` as :class:`AnalysisBase` Example ------- - >>> def rotation_matrix(mobile, ref): - >>> return mda.analysis.align.rotation_matrix(mobile, ref)[0] + .. code-block:: python - >>> rot = AnalysisFromFunction(rotation_matrix, trajectory, mobile, ref).run() - >>> print(rot.results) + def rotation_matrix(mobile, ref): + return mda.analysis.align.rotation_matrix(mobile, ref)[0] - Raises - ------ - ValueError : if ``function`` has the same kwargs as ``BaseAnalysis`` - """ + rot = AnalysisFromFunction(rotation_matrix, trajectory, + mobile, ref).run() + print(rot.results.timeseries) - def __init__(self, function, trajectory=None, *args, **kwargs): - """Parameters - ---------- - function : callable - function to evaluate at each frame - trajectory : mda.coordinates.Reader (optional) - trajectory to iterate over. If ``None`` the first AtomGroup found in - args and kwargs is used as a source for the trajectory. - *args : list - arguments for ``function`` - **kwargs : dict - arguments for ``function`` and ``AnalysisBase`` - .. versionchanged:: 1.0.0 - Support for directly passing the ``start``, ``stop``, and ``step`` - arguments has been removed. These should instead be passed - to :meth:`AnalysisFromFunction.run`. + .. versionchanged:: 1.0.0 + Support for directly passing the ``start``, ``stop``, and ``step`` + arguments has been removed. These should instead be passed + to :meth:`AnalysisFromFunction.run`. - """ + .. versionchanged:: 2.0.0 + Former :attr:`results` are now stored as :attr:`results.timeseries` + """ + + def __init__(self, function, trajectory=None, *args, **kwargs): if (trajectory is not None) and (not isinstance( trajectory, coordinates.base.ProtoReader)): args = (trajectory,) + args @@ -261,33 +376,67 @@ def __init__(self, function, trajectory=None, *args, **kwargs): super(AnalysisFromFunction, self).__init__(trajectory) def _prepare(self): - self.results = [] + self.results.timeseries = [] def _single_frame(self): - self.results.append(self.function(*self.args, **self.kwargs)) + self.results.timeseries.append(self.function(*self.args, + **self.kwargs)) def _conclude(self): - self.results = np.asarray(self.results) + self.results.frames = self.frames + self.results.times = self.times + self.results.timeseries = np.asarray(self.results.timeseries) def analysis_class(function): - """ - Transform a function operating on a single frame to an analysis class + r"""Transform a function operating on a single frame to an + :class:`AnalysisBase` class. + + Parameters + ---------- + function : callable + function to evaluate at each frame + + Attributes + ---------- + results.frames : numpy.ndarray + simulation frames used in analysis + results.times : numpy.ndarray + simulation times used in analysis + results.timeseries : numpy.ndarray + Results for each frame of the wrapped function, + stored after call to :meth:`AnalysisFromFunction.run`. + + Raises + ------ + ValueError + if ``function`` has the same ``kwargs`` as :class:`AnalysisBase` + + Examples + -------- - For an usage in a library we recommend the following style: + For use in a library, we recommend the following style + + .. code-block:: python + + def rotation_matrix(mobile, ref): + return mda.analysis.align.rotation_matrix(mobile, ref)[0] + RotationMatrix = analysis_class(rotation_matrix) + + It can also be used as a decorator + + .. code-block:: python - >>> def rotation_matrix(mobile, ref): - >>> return mda.analysis.align.rotation_matrix(mobile, ref)[0] - >>> RotationMatrix = analysis_class(rotation_matrix) + @analysis_class + def RotationMatrix(mobile, ref): + return mda.analysis.align.rotation_matrix(mobile, ref)[0] - It can also be used as a decorator: + rot = RotationMatrix(u.trajectory, mobile, ref).run(step=2) + print(rot.results.timeseries) - >>> @analysis_class - >>> def RotationMatrix(mobile, ref): - >>> return mda.analysis.align.rotation_matrix(mobile, ref)[0] - >>> rot = RotationMatrix(u.trajectory, mobile, ref).run(step=2) - >>> print(rot.results) + .. versionchanged:: 2.0.0 + Former ``results`` are now stored as ``results.timeseries`` """ class WrapperClass(AnalysisFromFunction): @@ -318,7 +467,8 @@ def _filter_baseanalysis_kwargs(function, kwargs): Raises ------ - ValueError : if ``function`` has the same kwargs as ``BaseAnalysis`` + ValueError + if ``function`` has the same ``kwargs`` as :class:`AnalysisBase` """ try: # pylint: disable=deprecated-method diff --git a/package/MDAnalysis/analysis/gnm.py b/package/MDAnalysis/analysis/gnm.py index d4fd2700f34..1f67b4d45f3 100644 --- a/package/MDAnalysis/analysis/gnm.py +++ b/package/MDAnalysis/analysis/gnm.py @@ -41,7 +41,9 @@ .. _Cookbook: https://github.com/MDAnalysis/MDAnalysisCookbook The basic approach is to pass a trajectory to :class:`GNMAnalysis` and then run -the analysis:: +the analysis: + +.. code-block:: python u = MDAnalysis.Universe(PSF, DCD) C = MDAnalysis.analysis.gnm.GNMAnalysis(u, ReportVector="output.txt") @@ -94,6 +96,8 @@ from .base import AnalysisBase +from MDAnalysis.analysis.base import Results + logger = logging.getLogger('MDAnalysis.analysis.GNM') @@ -204,27 +208,28 @@ class GNMAnalysis(AnalysisBase): universe : Universe Analyze the full trajectory in the universe. select : str (optional) - MDAnalysis selection string, default "protein and name CA" + MDAnalysis selection string cutoff : float (optional) Consider selected atoms within the cutoff as neighbors for the Gaussian network model. ReportVector : str (optional) filename to write eigenvectors to, by default no output is written - (``None``) Bonus_groups : tuple This is a tuple of selection strings that identify additional groups (such as ligands). The center of mass of each group will be added as a single point in the ENM (it is a popular way of treating small ligands such as drugs). You need to ensure that none of the atoms in `Bonus_groups` is contained in `selection` as this could lead to - double counting. No checks are applied. Default is ``None``. + double counting. No checks are applied. Attributes ---------- - results : list - GNM results per frame: - results = [(time,eigenvalues[1],eigenvectors[1]), - (time,eigenvalues[1],eigenvectors[1]), ...] + results.times : numpy.ndarray + simulation times used in analysis + results.eigenvalues : numpy.ndarray + calculated eigenvalues + results.eigenvectors : numpy.ndarray + calculated eigenvectors See Also -------- @@ -238,7 +243,9 @@ class GNMAnalysis(AnalysisBase): Changed `selection` keyword to `select` .. versionchanged:: 2.0.0 - Use :class:`~MDAnalysis.analysis.AnalysisBase` as parent class. + Use :class:`~MDAnalysis.analysis.AnalysisBase` as parent class and + store results as attributes ``times``, ``eigenvalues`` and + ``eigenvectors`` of the ``results`` attribute. """ def __init__(self, @@ -251,16 +258,18 @@ def __init__(self, self.u = universe self.select = select self.cutoff = cutoff - self.results = [] # final result + self.results = Results() + self.results.eigenvalues = [] + self.results.eigenvectors = [] self._timesteps = None # time for each frame self.ReportVector = ReportVector self.Bonus_groups = [self.u.select_atoms(item) for item in Bonus_groups] \ if Bonus_groups else [] self.ca = self.u.select_atoms(self.select) - def _generate_output(self, w, v, outputobject, time, matrix, - nmodes=2, ReportVector=None, counter=0): - """Appends eigenvalues and eigenvectors to results. + def _generate_output(self, w, v, outputobject, + ReportVector=None, counter=0): + """Appends time, eigenvalues and eigenvectors to results. This generates the output by adding eigenvalue and eigenvector data to an appendable object and optionally @@ -275,14 +284,13 @@ def _generate_output(self, w, v, outputobject, time, matrix, print( "", counter, - time, item[0] + 1, w[list_map[1]], item[1], file=oup) - outputobject.append((time, w[list_map[1]], v[list_map[1]])) - # outputobject.append((time, [ w[list_map[i]] for i in range(nmodes) ], - # [ v[list_map[i]] for i in range( nmodes) ] )) + + outputobject.eigenvalues.append(w[list_map[1]]) + outputobject.eigenvectors.append(v[list_map[1]]) def generate_kirchoff(self): """Generate the Kirchhoff matrix of contacts. @@ -318,9 +326,6 @@ def generate_kirchoff(self): return matrix - def _prepare(self): - self.timeseries = [] - def _single_frame(self): matrix = self.generate_kirchoff() try: @@ -336,13 +341,13 @@ def _single_frame(self): w, v, self.results, - self._ts.time, - matrix, ReportVector=self.ReportVector, counter=self._ts.frame) def _conclude(self): - self._timesteps = self.times + self.results.times = self.times + self.results.eigenvalues = np.asarray(self.results.eigenvalues) + self.results.eigenvectors = np.asarray(self.results.eigenvectors) class closeContactGNMAnalysis(GNMAnalysis): @@ -357,19 +362,27 @@ class closeContactGNMAnalysis(GNMAnalysis): universe : Universe Analyze the full trajectory in the universe. select : str (optional) - MDAnalysis selection string, default "protein" + MDAnalysis selection string cutoff : float (optional) Consider selected atoms within the cutoff as neighbors for the - Gaussian network model [4.5 Å]. + Gaussian network model. ReportVector : str (optional) filename to write eigenvectors to, by default no output is written - (``None``) weights : {"size", None} (optional) If set to "size" (the default) then weight the contact by :math:`1/\sqrt{N_i N_j}` where :math:`N_i` and :math:`N_j` are the number of atoms in the residues :math:`i` and :math:`j` that contain the atoms that form a contact. + Attributes + ---------- + results.times : numpy.ndarray + simulation times used in analysis + results.eigenvalues : numpy.ndarray + calculated eigenvalues + results.eigenvectors : numpy.ndarray + calculated eigenvectors + Notes ----- The `MassWeight` option has now been removed. @@ -390,7 +403,9 @@ class closeContactGNMAnalysis(GNMAnalysis): Changed `selection` keyword to `select` .. versionchanged:: 2.0.0 - Use :class:`~MDAnalysis.analysis.AnalysisBase` as parent class. + Use :class:`~MDAnalysis.analysis.AnalysisBase` as parent class and + store results as attributes ``times``, ``eigenvalues`` and + ``eigenvectors`` of the `results` attribute. """ def __init__(self, diff --git a/package/MDAnalysis/analysis/lineardensity.py b/package/MDAnalysis/analysis/lineardensity.py index 0ee87a5d720..78ce2fe4433 100644 --- a/package/MDAnalysis/analysis/lineardensity.py +++ b/package/MDAnalysis/analysis/lineardensity.py @@ -25,14 +25,14 @@ =========================================================== A tool to compute mass and charge density profiles along the three -cartesian axes of the simulation cell. Works only for orthorombic, +cartesian axes [xyz] of the simulation cell. Works only for orthorombic, fixed volume cells (thus for simulations in canonical NVT ensemble). """ import os.path as path import numpy as np -from MDAnalysis.analysis.base import AnalysisBase +from MDAnalysis.analysis.base import AnalysisBase, Results class LinearDensity(AnalysisBase): @@ -49,25 +49,35 @@ class LinearDensity(AnalysisBase): Bin width in Angstrom used to build linear density histograms. Defines the resolution of the resulting density profile (smaller --> higher resolution) - verbose : bool (optional) - Show detailed progress of the calculation if set to ``True``; the - default is ``False``. + verbose : bool, optional + Show detailed progress of the calculation if set to ``True`` Attributes ---------- - results : dict - Keys 'x', 'y', and 'z' for the three directions. Under these - keys, find 'pos', 'pos_std' (mass-weighted density and - standard deviation), 'char', 'char_std' (charge density and - its standard deviation), 'slice_volume' (volume of bin). + results.x.dim : int + index of the [xyz] axes + results.x.pos : numpy.ndarray + mass density in [xyz] direction + results.x.pos_std : numpy.ndarray + standard deviation of the mass density in [xyz] direction + results.x.char : numpy.ndarray + charge density in [xyz] direction + results.x.char_std : numpy.ndarray + standard deviation of the charge density in [xyz] direction + results.x.slice_volume : float + volume of bin in [xyz] direction Example ------- - First create a LinearDensity object by supplying a selection, - then use the :meth:`run` method:: + First create a ``LinearDensity`` object by supplying a selection, + then use the :meth:`run` method. Finally access the results + stored in results, i.e. the mass density in the x direction. - ldens = LinearDensity(selection) - ldens.run() + .. code-block:: python + + ldens = LinearDensity(selection) + ldens.run() + print(ldens.results.x.pos) .. versionadded:: 0.14.0 @@ -81,6 +91,11 @@ class LinearDensity(AnalysisBase): .. versionchanged:: 1.0.0 Changed `selection` keyword to `select` + + .. versionchanged:: 2.0.0 + Results are now instances of + :class:`~MDAnalysis.core.analysis.Results` allowing access + via key and attribute. """ def __init__(self, select, grouping='atoms', binsize=0.25, **kwargs): @@ -96,8 +111,10 @@ def __init__(self, select, grouping='atoms', binsize=0.25, **kwargs): # AtomGroup.wrap()) self.grouping = grouping - # Dictionary containing results - self.results = {'x': {'dim': 0}, 'y': {'dim': 1}, 'z': {'dim': 2}} + # Initiate result instances + self.results["x"] = Results(dim=0) + self.results["y"] = Results(dim=1) + self.results["z"] = Results(dim=2) # Box sides self.dimensions = self._universe.dimensions[:3] self.volume = np.prod(self.dimensions) @@ -114,9 +131,9 @@ def __init__(self, select, grouping='atoms', binsize=0.25, **kwargs): # Initialize results array with zeros for dim in self.results: idx = self.results[dim]['dim'] - self.results[dim].update({'slice volume': slices_vol[idx]}) + self.results[dim]['slice_volume'] = slices_vol[idx] for key in self.keys: - self.results[dim].update({key: np.zeros(self.nbins)}) + self.results[dim][key] = np.zeros(self.nbins) # Variables later defined in _prepare() method self.masses = None @@ -177,7 +194,7 @@ def _single_frame(self): def _conclude(self): k = 6.022e-1 # divide by avodagro and convert from A3 to cm3 - # Average results over the number of configurations + # 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 @@ -188,10 +205,9 @@ def _conclude(self): '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 + norm = k * self.results[dim]['slice_volume'] + for key in self.keys: + self.results[dim][key] /= norm def _add_other_results(self, other): # For parallel analysis diff --git a/package/MDAnalysis/analysis/polymer.py b/package/MDAnalysis/analysis/polymer.py index 78d99131e34..21973db3057 100644 --- a/package/MDAnalysis/analysis/polymer.py +++ b/package/MDAnalysis/analysis/polymer.py @@ -163,12 +163,11 @@ class PersistenceLength(AnalysisBase): List of AtomGroups. Each should represent a single polymer chain, ordered in the correct order. verbose : bool (optional) - Show detailed progress of the calculation if set to ``True``; the - default is ``False``. + Show detailed progress of the calculation if set to ``True``. Attributes ---------- - results : numpy.ndarray + results.bond_autocorrelation : numpy.ndarray the measured bond autocorrelation lb : float the average bond length @@ -187,6 +186,8 @@ class PersistenceLength(AnalysisBase): The run method now automatically performs the exponential fit .. versionchanged:: 1.0.0 Deprecated :meth:`PersistenceLength.perform_fit` has now been removed. + .. versionchanged:: 2.0.0 + Former ``results`` are now stored as ``results.bond_autocorrelation`` """ def __init__(self, atomgroups, **kwargs): super(PersistenceLength, self).__init__( @@ -224,7 +225,7 @@ def _conclude(self): norm = np.linspace(n - 1, 1, n - 1) norm *= len(self._atomgroups) * self.n_frames - self.results = self._results / norm + self.results.bond_autocorrelation = self._results / norm self._calc_bond_length() self._perform_fit() @@ -241,12 +242,13 @@ def _calc_bond_length(self): def _perform_fit(self): """Fit the results to an exponential decay""" try: - self.results + self.results.bond_autocorrelation except AttributeError: raise NoDataError("Use the run method first") from None - self.x = np.arange(len(self.results)) * self.lb + self.x = np.arange(len(self.results.bond_autocorrelation)) * self.lb - self.lp = fit_exponential_decay(self.x, self.results) + self.lp = fit_exponential_decay(self.x, + self.results.bond_autocorrelation) self.fit = np.exp(-self.x/self.lp) @@ -265,8 +267,13 @@ def plot(self, ax=None): import matplotlib.pyplot as plt if ax is None: fig, ax = plt.subplots() - ax.plot(self.x, self.results, 'ro', label='Result') - ax.plot(self.x, self.fit, label='Fit') + ax.plot(self.x, + self.results.bond_autocorrelation, + 'ro', + label='Result') + ax.plot(self.x, + self.fit, + label='Fit') ax.set_xlabel(r'x') ax.set_ylabel(r'$C(x)$') ax.set_xlim(0.0, 40 * self.lb) diff --git a/package/doc/sphinx/source/documentation_pages/analysis_modules.rst b/package/doc/sphinx/source/documentation_pages/analysis_modules.rst index 95605bace5b..274c1f56779 100644 --- a/package/doc/sphinx/source/documentation_pages/analysis_modules.rst +++ b/package/doc/sphinx/source/documentation_pages/analysis_modules.rst @@ -5,28 +5,43 @@ Analysis modules **************** The :mod:`MDAnalysis.analysis` module contains code to carry out specific -analysis functionality. It is based on the core functionality (i.e. trajectory +analysis functionality for MD trajectories. +It is based on the core functionality (i.e. trajectory I/O, selections etc). The analysis modules can be used as examples for how to use MDAnalysis but also as working code for research projects; typically all contributed code has been used by the authors in their own work. +An analysis using the available modules +usually follows the same structure -Please see the individual module documentation for additional references and -citation information. +#. Import the desired module, since analysis modules are not imported + by default. +#. Initialize the module previously imported. +#. Run the analysis, optionally for specific trajectory slices +#. Access the analysis from the `results` attribute -These modules are not imported by default; in order to use them one has to -import them from :mod:`MDAnalysis.analysis`, for instance :: +.. code-block:: python - import MDAnalysis.analysis.align + from MDAnalysis.analysis import ExampleAnalysisModule # (e.g. RMSD) + + analysis_obj = ExampleAnalysisModule(universe, ...) + analysis_obj.run(start_frame, stop_frame, step) + print(analysis_obj.results) + + +Please see the individual module documentation for any specific caveats +and also read and cite the reference papers associated with these algorithms. .. rubric:: Additional dependencies Some of the modules in :mod:`MDAnalysis.analysis` require additional Python packages to enable full functionality. For example, :mod:`MDAnalysis.analysis.encore` provides more options if `scikit-learn`_ is -installed. These package are *not automatically installed* with -:program:`pip`(although one can add the ``[analysis]`` requirement to the -:program:`pip` command line to force their installation). If you install -MDAnalysis with :program:`conda` (see :ref:`installation-instructions`) then a +installed. If you installed MDAnalysis with +:program:`pip` (see :ref:`installation-instructions`) +these packages are *not automatically installed*. +Although, one can add the ``[analysis]`` tag to the +:program:`pip` command to force their installation. If you installed +MDAnalysis with :program:`conda` then a *full set of dependencies* is automatically installed. Other modules require external programs. For instance, the @@ -38,10 +53,13 @@ corresponding MDAnalysis module. .. _scikit-learn: http://scikit-learn.org/ .. _HOLE: http://www.holeprogram.org/ - Building blocks for Analysis ============================ +The building block for the analysis modules is +:class:`MDAnalysis.analysis.base.AnalysisBase`. +To build your own analysis class start by reading the documentation. + .. toctree:: :maxdepth: 1 diff --git a/testsuite/MDAnalysisTests/analysis/test_base.py b/testsuite/MDAnalysisTests/analysis/test_base.py index 95a7c8a19bf..71ea1aea3e3 100644 --- a/testsuite/MDAnalysisTests/analysis/test_base.py +++ b/testsuite/MDAnalysisTests/analysis/test_base.py @@ -20,6 +20,8 @@ # MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations. # J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787 # +from collections import UserDict + import pytest import numpy as np @@ -33,6 +35,48 @@ from MDAnalysisTests.util import no_deprecated_call +class Test_Results: + + @pytest.fixture + def results(self): + return base.Results(a=1, b=2) + + def test_get(self, results): + assert results.a == results["a"] == 1 + + def test_no_attr(self, results): + with pytest.raises(AttributeError): + results.c + + def test_set_attr(self, results): + value = [1, 2, 3, 4] + results.c = value + assert results.c == results["c"] == value + + def test_set_key(self, results): + value = [1, 2, 3, 4] + results["c"] = value + assert results.c == results["c"] == value + + @pytest.mark.parametrize('key', dir(UserDict) + ["data"]) + def test_existing_dict_attr(self, results, key): + msg = f"'{key}' is a protected dictionary attribute" + with pytest.raises(AttributeError, match=key): + results[key] = None + + @pytest.mark.parametrize('key', dir(UserDict) + ["data"]) + def test_wrong_init_type(self, key): + msg = f"'{key}' is a protected dictionary attribute" + with pytest.raises(AttributeError, match=msg): + base.Results(**{key: None}) + + @pytest.mark.parametrize('key', ("0123", "0j", "1.1", "{}", "a b")) + def test_weird_key(self, results, key): + msg = f"'{key}' is not a valid attribute" + with pytest.raises(ValueError, match=msg): + results[key] = None + + class FrameAnalysis(base.AnalysisBase): """Just grabs frame numbers of frames it goes over""" @@ -152,6 +196,11 @@ def simple_function(mobile): return mobile.center_of_geometry() +def test_results_type(u): + an = FrameAnalysis(u.trajectory) + assert type(an.results) == base.Results + + @pytest.mark.parametrize('start, stop, step, nframes', [ (None, None, 2, 49), (None, 50, 2, 25), @@ -168,17 +217,25 @@ def test_AnalysisFromFunction(u, start, stop, step, nframes): ana3 = base.AnalysisFromFunction(simple_function, u.trajectory, u.atoms) ana3.run(start=start, stop=stop, step=step) - results = [] + frames = [] + times = [] + timeseries = [] for ts in u.trajectory[start:stop:step]: - results.append(simple_function(u.atoms)) + frames.append(ts.frame) + times.append(ts.time) + timeseries.append(simple_function(u.atoms)) - results = np.asarray(results) + frames = np.asarray(frames) + times = np.asarray(times) + timeseries = np.asarray(timeseries) - assert np.size(results, 0) == nframes + assert np.size(timeseries, 0) == nframes for ana in (ana1, ana2, ana3): - assert_equal(results, ana.results) + assert_equal(frames, ana.results.frames) + assert_equal(times, ana.results.times) + assert_equal(timeseries, ana.results.timeseries) def mass_xyz(atomgroup1, atomgroup2, masses): @@ -191,7 +248,7 @@ def test_AnalysisFromFunction_args_content(u): another = mda.Universe(TPR, XTC).select_atoms("protein") ans = base.AnalysisFromFunction(mass_xyz, protein, another, masses) assert len(ans.args) == 3 - result = np.sum(ans.run().results) + result = np.sum(ans.run().results.timeseries) assert_almost_equal(result, -317054.67757345125, decimal=6) assert (ans.args[0] is protein) and (ans.args[1] is another) assert ans._trajectory is protein.universe.trajectory @@ -211,7 +268,7 @@ def test_analysis_class(): results.append(simple_function(u.atoms)) results = np.asarray(results) - assert_equal(results, ana.results) + assert_equal(results, ana.results.timeseries) with pytest.raises(ValueError): ana_class(2) diff --git a/testsuite/MDAnalysisTests/analysis/test_gnm.py b/testsuite/MDAnalysisTests/analysis/test_gnm.py index 0bb25e1b4c7..6521c08eb86 100644 --- a/testsuite/MDAnalysisTests/analysis/test_gnm.py +++ b/testsuite/MDAnalysisTests/analysis/test_gnm.py @@ -43,10 +43,9 @@ def test_gnm(universe, tmpdir): gnm = mda.analysis.gnm.GNMAnalysis(universe, ReportVector=output) gnm.run() result = gnm.results - assert len(result) == 10 - time, eigenvalues, eigenvectors = zip(*result) - assert_almost_equal(time, np.arange(0, 1000, 100), decimal=4) - assert_almost_equal(eigenvalues, + assert len(result.times) == 10 + assert_almost_equal(gnm.results.times, np.arange(0, 1000, 100), decimal=4) + assert_almost_equal(gnm.results.eigenvalues, [2.0287113e-15, 4.1471575e-15, 1.8539533e-15, 4.3810359e-15, 3.9607304e-15, 4.1289113e-15, 2.5501084e-15, 4.0498182e-15, 4.2058769e-15, 3.9839431e-15]) @@ -56,10 +55,9 @@ def test_gnm_run_step(universe): gnm = mda.analysis.gnm.GNMAnalysis(universe) gnm.run(step=3) result = gnm.results - assert len(result) == 4 - time, eigenvalues, eigenvectors = zip(*result) - assert_almost_equal(time, np.arange(0, 1200, 300), decimal=4) - assert_almost_equal(eigenvalues, + assert len(result.times) == 4 + assert_almost_equal(gnm.results.times, np.arange(0, 1200, 300), decimal=4) + assert_almost_equal(gnm.results.eigenvalues, [2.0287113e-15, 4.3810359e-15, 2.5501084e-15, 3.9839431e-15]) @@ -94,10 +92,9 @@ def test_closeContactGNMAnalysis(universe): gnm = mda.analysis.gnm.closeContactGNMAnalysis(universe, weights="size") gnm.run(stop=2) result = gnm.results - assert len(result) == 2 - time, eigenvalues, eigenvectors = zip(*result) - assert_almost_equal(time, (0, 100), decimal=4) - assert_almost_equal(eigenvalues, [0.1502614, 0.1426407]) + assert len(result.times) == 2 + assert_almost_equal(gnm.results.times, (0, 100), decimal=4) + assert_almost_equal(gnm.results.eigenvalues, [0.1502614, 0.1426407]) gen = gnm.generate_kirchoff() assert_almost_equal(gen[0], [16.326744128018923, -2.716098853586913, -1.94736842105263, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, @@ -121,10 +118,9 @@ def test_closeContactGNMAnalysis_weights_None(universe): gnm = mda.analysis.gnm.closeContactGNMAnalysis(universe, weights=None) gnm.run(stop=2) result = gnm.results - assert len(result) == 2 - time, eigenvalues, eigenvectors = zip(*result) - assert_almost_equal(time, (0, 100), decimal=4) - assert_almost_equal(eigenvalues, [2.4328739, 2.2967251]) + assert len(result.times) == 2 + assert_almost_equal(gnm.results.times, (0, 100), decimal=4) + assert_almost_equal(gnm.results.eigenvalues, [2.4328739, 2.2967251]) gen = gnm.generate_kirchoff() assert_almost_equal(gen[0], [303.0, -58.0, -37.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -1.0, diff --git a/testsuite/MDAnalysisTests/analysis/test_persistencelength.py b/testsuite/MDAnalysisTests/analysis/test_persistencelength.py index 864a2aed06c..ebc941fc729 100644 --- a/testsuite/MDAnalysisTests/analysis/test_persistencelength.py +++ b/testsuite/MDAnalysisTests/analysis/test_persistencelength.py @@ -62,14 +62,14 @@ def test_ag_ValueError(self, u): polymer.PersistenceLength(ags) def test_run(self, p_run): - assert len(p_run.results) == 280 + assert len(p_run.results.bond_autocorrelation) == 280 def test_lb(self, p_run): assert_almost_equal(p_run.lb, 1.485, 3) def test_fit(self, p_run): assert_almost_equal(p_run.lp, 6.504, 3) - assert len(p_run.fit) == len(p_run.results) + assert len(p_run.fit) == len(p_run.results.bond_autocorrelation) def test_raise_NoDataError(self, p): #Ensure that a NoDataError is raised if perform_fit() @@ -90,7 +90,7 @@ def test_plot_with_ax(self, p_run): ax2 = p_run.plot(ax=ax) assert ax2 is ax - + def test_current_axes(self, p_run): fig, ax = plt.subplots() ax2 = p_run.plot(ax=None)