diff --git a/.travis.yml b/.travis.yml index 768474f4..8aa0c95b 100644 --- a/.travis.yml +++ b/.travis.yml @@ -18,16 +18,15 @@ env: - SETUP_CMD="pmda --cov pmda" # mdanalysis develop from source (see below), which needs # minimal CONDA_MDANALYSIS_DEPENDENCIES - #- CONDA_DEPENDENCIES="mdanalysis mdanalysistests dask joblib mock codecov cython hypothesis sphinx" - #- CONDA_MDANALYSIS_DEPENDENCIES="cython mmtf-python six biopython networkx scipy griddataformats gsd hypothesis" - - CONDA_MDANALYSIS_DEPENDENCIES="mdanalysis mdanalysistests" - # pin to matplotlib 3.2 should be removed once mda > 1.0.0 is available. - - CONDA_DEPENDENCIES="${CONDA_MDANALYSIS_DEPENDENCIES} dask distributed joblib mock codecov matplotlib=3.2" + # - CONDA_DEPENDENCIES="mdanalysis mdanalysistests dask joblib pytest-pep8 mock codecov cython hypothesis sphinx" + # - CONDA_MDANALYSIS_DEPENDENCIES="cython mmtf-python six biopython networkx scipy griddataformats gsd hypothesis" + - CONDA_MDANALYSIS_DEPENDENCIES="mmtf-python biopython networkx cython matplotlib scipy griddataformats hypothesis gsd" + - CONDA_DEPENDENCIES="${CONDA_MDANALYSIS_DEPENDENCIES} dask distributed joblib pytest-pep8 mock codecov" - CONDA_CHANNELS='conda-forge' - CONDA_CHANNEL_PRIORITY=True # install development version of MDAnalysis (needed until the test # files for analysis.rdf are available in release 0.19.0) - #- PIP_DEPENDENCIES="git+https://github.com/MDAnalysis/mdanalysis#egg=mdanalysis&subdirectory=package git+https://github.com/MDAnalysis/mdanalysis#egg=mdanalysistests&subdirectory=testsuite" + - PIP_DEPENDENCIES="git+https://github.com/MDAnalysis/mdanalysis#egg=mdanalysis&subdirectory=package git+https://github.com/MDAnalysis/mdanalysis#egg=mdanalysistests&subdirectory=testsuite" - NUMPY_VERSION=stable - BUILD_CMD='python setup.py develop' - CODECOV=false diff --git a/CHANGELOG b/CHANGELOG index 58b1ab97..69cfae53 100644 --- a/CHANGELOG +++ b/CHANGELOG @@ -35,6 +35,14 @@ Fixes Changes * requires MDAnalysis >= 1.0.0 (#122) * dropped official support for Python 3.5 (2.7 and >= 3.6 are supported) + * requires MDAnalysis >= 2.0.0 (#132) + * dropped official support for Python<=3.5 (#132) + * with support of serialization of Universe, in ParallelAnalysisBase: (#132) + - we pickle/unpickle Universes, instead of creating new ones. + - __init__ only takes Universe as then argument. + - _single_frame takes no argument. + - timing for generating universe is ditched. + - timing for creating dask graph is added. 10/14/2019 VOD555, nawtrey diff --git a/pmda/contacts.py b/pmda/contacts.py index 510cc0e9..b4b08913 100644 --- a/pmda/contacts.py +++ b/pmda/contacts.py @@ -182,8 +182,6 @@ def is_any_closer(r, r0, dist=2.5): :inherited-members: """ -from __future__ import absolute_import, division - import MDAnalysis as mda from MDAnalysis.analysis.contacts import (contact_matrix, hard_cut_q, radius_cut_q, soft_cut_q) @@ -237,7 +235,8 @@ def __init__(self, respective functions for reasonable values. """ universe = mobiles[0].universe - super(Contacts, self).__init__(universe, mobiles) + super().__init__(universe) + self._mobiles = mobiles if method == 'hard_cut': self.fraction_contacts = hard_cut_q @@ -267,13 +266,13 @@ def __init__(self, def _prepare(self): self.timeseries = None - def _single_frame(self, ts, atomgroups): - grA, grB = atomgroups + def _single_frame(self): + grA, grB = self._mobiles # compute distance array for a frame d = distance_array(grA.positions, grB.positions) y = np.empty(len(self.r0) + 1) - y[0] = ts.frame + y[0] = self._ts.frame for i, (initial_contacts, r0) in enumerate(zip(self.initial_contacts, self.r0)): # select only the contacts that were formed in the reference state diff --git a/pmda/custom.py b/pmda/custom.py index 23056def..cab045a6 100644 --- a/pmda/custom.py +++ b/pmda/custom.py @@ -16,9 +16,6 @@ same universe and return a value. """ -from __future__ import absolute_import - -from MDAnalysis.core.groups import AtomGroup from MDAnalysis.core.universe import Universe from MDAnalysis.coordinates.base import ProtoReader import numpy as np @@ -44,7 +41,7 @@ class AnalysisFromFunction(ParallelAnalysisBase): >>> return mda.analysis.align.rotation_matrix(mobile, ref)[0] >>> # now run an analysis using the standalone function >>> rot = AnalysisFromFunction(rotation_matrix, - trajectory, mobile, ref).run() + universe, mobile, ref).run() >>> print(rot.results) Raises @@ -64,40 +61,24 @@ def __init__(self, function, universe, *args, **kwargs): to be 'mobile' Atomgroups if they belong to the same universe. All other Atomgroups are assumed to be reference. Here 'mobile' means they will be iterated over. - Universe : :class:`~MDAnalysis.core.groups.Universe` + universe : :class:`~MDAnalysis.core.groups.Universe` a :class:`MDAnalysis.core.groups.Universe` (the - `atomgroups` must belong to this Universe) + `atomgroups` in other args must belong to this Universe) *args : list arguments for ``function`` **kwargs : dict - keyword arguments for ``function``. keyword arguments with name - 'universe' or 'atomgroups' will be ignored! Mobile atomgroups to - analyze can not be passed as keyword arguments currently. - + keyword arguments for ``function``. """ - self.function = function - - # collect all atomgroups with the same trajectory object as universe - trajectory = universe.trajectory - arg_ags = [] - self.other_args = [] - for arg in args: - if isinstance(arg, - AtomGroup) and arg.universe.trajectory == trajectory: - arg_ags.append(arg) - else: - self.other_args.append(arg) - - super(AnalysisFromFunction, self).__init__(universe, arg_ags) + super().__init__(universe) + self.args = args self.kwargs = kwargs def _prepare(self): self.results = [] - def _single_frame(self, ts, atomgroups): - args = atomgroups + self.other_args - return self.function(*args, **self.kwargs) + def _single_frame(self): + return self.function(*self.args, **self.kwargs) def _conclude(self): self.results = np.concatenate(self._results) @@ -132,7 +113,7 @@ def analysis_class(function): >>> def RotationMatrix(mobile, ref): >>> return mda.analysis.align.rotation_matrix(mobile, ref)[0] - >>> rot = RotationMatrix(u.trajectory, mobile, ref, step=2).run() + >>> rot = RotationMatrix(u, mobile, ref, step=2).run() >>> print(rot.results) See Also @@ -144,13 +125,14 @@ def analysis_class(function): class WrapperClass(AnalysisFromFunction): """Custom Analysis Function""" - def __init__(self, trajectory=None, *args, **kwargs): - if not (isinstance(trajectory, ProtoReader) or isinstance( - trajectory, Universe)): - print(type(trajectory)) + def __init__(self, universe=None, *args, **kwargs): + if not isinstance(universe, Universe): + print(type(universe)) raise ValueError( - "First argument needs to be an MDAnalysis reader object.") - super(WrapperClass, self).__init__(function, trajectory, *args, - **kwargs) + "First argument needs to be an MDAnalysis Universe.") + super().__init__(function, + universe, + *args, + **kwargs) return WrapperClass diff --git a/pmda/density.py b/pmda/density.py index 461ae1f7..65050f5c 100644 --- a/pmda/density.py +++ b/pmda/density.py @@ -22,9 +22,6 @@ MDAnalysis.analysis.density.density_from_Universe """ - -from __future__ import absolute_import - import numpy as np from MDAnalysis.lib.util import fixedwidth_bins @@ -241,7 +238,7 @@ def __init__(self, atomgroup, delta=1.0, atomselection=None, parameters=None, gridcenter=None, xdim=None, ydim=None, zdim=None): u = atomgroup.universe - super(DensityAnalysis, self).__init__(u, (atomgroup, )) + super().__init__(u) self._atomgroup = atomgroup self._delta = delta self._atomselection = atomselection @@ -259,10 +256,14 @@ def __init__(self, atomgroup, delta=1.0, atomselection=None, elif not updating and atomselection is not None: raise ValueError("""With updating=False, the atomselection='{}' is not used and should be None""".format(atomselection)) + elif updating and atomselection is not None: + self._select_atomgroup = atomgroup.select_atoms(atomselection, + updating=True) + else: + self._select_atomgroup = atomgroup def _prepare(self): - coord = self.current_coordinates(self._atomgroup, self._atomselection, - self._updating) + coord = self._select_atomgroup.positions if self._gridcenter is not None: # Generate a copy of smin/smax from coords to later check if the # defined box might be too small for the selection @@ -294,10 +295,10 @@ def _prepare(self): self._arange = arange self._bins = bins - def _single_frame(self, ts, atomgroups): - coord = self.current_coordinates(atomgroups[0], self._atomselection, - self._updating) - result = np.histogramdd(coord, bins=self._bins, range=self._arange, + def _single_frame(self): + result = np.histogramdd(self._select_atomgroup.positions, + bins=self._bins, + range=self._arange, normed=False) return result[0] @@ -330,12 +331,3 @@ def _reduce(res, result_single_frame): else: res += result_single_frame return res - - @staticmethod - def current_coordinates(atomgroup, atomselection, updating): - """Retrieves the current coordinates of all atoms in the chosen atom - selection. - Note: currently required to allow for updating selections""" - ag = atomgroup if not updating else atomgroup.select_atoms( - atomselection) - return ag.positions diff --git a/pmda/hbond_analysis.py b/pmda/hbond_analysis.py index db2ab4d4..41a73c09 100644 --- a/pmda/hbond_analysis.py +++ b/pmda/hbond_analysis.py @@ -26,8 +26,6 @@ :members: """ -from __future__ import absolute_import, division - import numpy as np from MDAnalysis.lib.distances import capped_distance, calc_angles @@ -174,7 +172,8 @@ def __init__(self, universe, donors_sel=None, hydrogens_sel=None, """ ag = universe.atoms - super(HydrogenBondAnalysis, self).__init__(universe, (ag, )) + super().__init__(universe) + self._atomgroup = ag self.donors_sel = donors_sel self.hydrogens_sel = hydrogens_sel self.acceptors_sel = acceptors_sel @@ -220,7 +219,7 @@ def guess_hydrogens(self, selection='all', max_mass=1.1, min_charge=0.3): :attr:`hydrogens_sel`. """ - u = self._universe() + u = self._universe ag = u.select_atoms(selection) hydrogens_ag = ag[ np.logical_and( @@ -284,7 +283,7 @@ def guess_donors(self, selection='all', max_charge=-0.5): hydrogens_sel = self.guess_hydrogens() else: hydrogens_sel = self.hydrogens_sel - u = self._universe() + u = self._universe hydrogens_ag = u.select_atoms(hydrogens_sel) ag = hydrogens_ag.residues.atoms.select_atoms( @@ -337,7 +336,7 @@ def guess_acceptors(self, selection='all', max_charge=-0.5): attribute :attr:`acceptors_sel`. """ - u = self._universe() + u = self._universe ag = u.select_atoms(selection) acceptors_ag = ag[ag.charges < max_charge] acceptors_list = np.unique( @@ -395,7 +394,7 @@ def _get_dh_pairs(self, u): return donors, hydrogens def _prepare(self): - u = mda.Universe(self._top, self._traj) + u = self._universe self.hbonds = [] self.frames = np.arange(self.start, self.stop, self.step) self.timesteps = (self.frames*u.trajectory.dt) + u.trajectory[0].time @@ -412,10 +411,9 @@ def _prepare(self): self._donors_ids = donors.ids self._hydrogens_ids = hydrogens.ids - def _single_frame(self, ts, atomgroups): - u = atomgroups[0].universe - - box = ts.dimensions + def _single_frame(self): + u = self._universe + box = self._ts.dimensions # Update donor-hydrogen pairs if necessary if self.update_selections: @@ -461,7 +459,7 @@ def _single_frame(self, ts, atomgroups): # Store data on hydrogen bonds found at this frame hbonds = [[], [], [], [], [], []] - hbonds[0].extend(np.full_like(hbond_donors, ts.frame)) + hbonds[0].extend(np.full_like(hbond_donors, self._ts.frame)) hbonds[1].extend(hbond_donors.ids) hbonds[2].extend(hbond_hydrogens.ids) hbonds[3].extend(hbond_acceptors.ids) @@ -510,7 +508,7 @@ def count_by_type(self): resname and atom type of the donor and acceptor atoms in a hydrogen bond. """ - u = self._universe() + u = self._universe d = u.atoms[self.hbonds[:, 1].astype(np.int)] a = u.atoms[self.hbonds[:, 3].astype(np.int)] @@ -543,7 +541,7 @@ def count_by_ids(self): hydrogen atom id and acceptor atom id in a hydrogen bond. """ - u = self._universe() + u = self._universe d = u.atoms[self.hbonds[:, 1].astype(np.int)] h = u.atoms[self.hbonds[:, 2].astype(np.int)] a = u.atoms[self.hbonds[:, 3].astype(np.int)] @@ -560,14 +558,6 @@ def count_by_ids(self): return unique_hbonds - def _universe(self): - # A Universe containing position information is needed for guessing - # donors and acceptors. - u = mda.Universe(self._top) - if not hasattr(u.atoms, 'positions'): - u.load_new(self._positions) - return u - @staticmethod def _reduce(res, result_single_frame): """ Use numpy array append to combine results""" diff --git a/pmda/leaflet.py b/pmda/leaflet.py index b0e0bbed..eac91da7 100644 --- a/pmda/leaflet.py +++ b/pmda/leaflet.py @@ -71,11 +71,10 @@ class LeafletFinder(ParallelAnalysisBase): """ def __init__(self, universe, atomgroups): + super().__init__(universe) self._atomgroup = atomgroups self._results = list() - super(LeafletFinder, self).__init__(universe, (atomgroups,)) - def _find_connected_components(self, data, cutoff=15.0): """Perform the Connected Components discovery for the atoms in data. @@ -157,7 +156,7 @@ def _find_connected_components(self, data, cutoff=15.0): return comp # pylint: disable=arguments-differ - def _single_frame(self, ts, atomgroups, scheduler_kwargs, n_jobs, + def _single_frame(self, scheduler_kwargs, n_jobs, cutoff=15.0): """Perform computation on a single trajectory frame. @@ -187,7 +186,7 @@ def _single_frame(self, ts, atomgroups, scheduler_kwargs, n_jobs, """ # Get positions of the atoms in the atomgroup and find their number. - atoms = ts.positions[atomgroups.indices] + atoms = self._atomgroup.positions matrix_size = atoms.shape[0] arranged_coord = list() part_size = int(matrix_size / n_jobs) @@ -199,10 +198,13 @@ def _single_frame(self, ts, atomgroups, scheduler_kwargs, n_jobs, [i, j])) # Distribute the data over the available cores, apply the map function # and execute. - parAtoms = db.from_sequence(arranged_coord, - npartitions=len(arranged_coord)) - parAtomsMap = parAtoms.map_partitions(self._find_connected_components, - cutoff=cutoff) + with timeit() as prepare_dask: + parAtoms = db.from_sequence(arranged_coord, + npartitions=len(arranged_coord)) + parAtomsMap = parAtoms.map_partitions( + self._find_connected_components, + cutoff=cutoff) + self.prepare_dask_total += prepare_dask.elapsed Components = parAtomsMap.compute(**scheduler_kwargs) # Gather the results and start the reduction. TODO: think if it can go @@ -277,11 +279,11 @@ def run(self, if scheduler == 'multiprocessing': scheduler_kwargs['num_workers'] = n_jobs - with timeit() as b_universe: - universe = mda.Universe(self._top, self._traj) - start, stop, step = self._trajectory.check_slice_indices( start, stop, step) + + self.prepare_dask_total = 0 + with timeit() as total: with timeit() as prepare: self._prepare() @@ -291,15 +293,14 @@ def run(self, times_io = [] for frame in range(start, stop, step): with timeit() as b_io: - ts = universe.trajectory[frame] + self._ts = self._universe.trajectory[frame] times_io.append(b_io.elapsed) with timeit() as b_compute: components = self. \ - _single_frame(ts=ts, - atomgroups=self._atomgroup, - scheduler_kwargs=scheduler_kwargs, - n_jobs=n_jobs, - cutoff=cutoff) + _single_frame( + scheduler_kwargs=scheduler_kwargs, + n_jobs=n_jobs, + cutoff=cutoff) timings.append(b_compute.elapsed) leaflet1 = self._atomgroup[components[0]] leaflet2 = self._atomgroup[components[1]] @@ -308,7 +309,7 @@ def run(self, self._conclude() self.timing = Timing(times_io, np.hstack(timings), total.elapsed, - b_universe.elapsed, prepare.elapsed, + prepare.elapsed, self.prepare_dask_total, conclude.elapsed) return self diff --git a/pmda/parallel.py b/pmda/parallel.py index 361edc4d..6c13b436 100644 --- a/pmda/parallel.py +++ b/pmda/parallel.py @@ -14,12 +14,10 @@ classes. """ -from __future__ import absolute_import, division from contextlib import contextmanager import warnings import time -from six.moves import range import MDAnalysis as mda from dask.delayed import delayed import dask @@ -35,7 +33,7 @@ class Timing(object): store various timeing results of obtained during a parallel analysis run """ - def __init__(self, io, compute, total, universe, prepare, + def __init__(self, io, compute, total, prepare, prepare_dask, conclude, wait=None, io_block=None, compute_block=None): self._io = io @@ -44,8 +42,8 @@ def __init__(self, io, compute, total, universe, prepare, self._compute_block = compute_block self._total = total self._cumulate = np.sum(io) + np.sum(compute) - self._universe = universe self._prepare = prepare + self._prepare_dask = prepare_dask self._conclude = conclude self._wait = wait @@ -83,16 +81,16 @@ def cumulate_time(self): """ return self._cumulate - @property - def universe(self): - """time to create a universe for each block""" - return self._universe - @property def prepare(self): """time to prepare""" return self._prepare + @property + def prepare_dask(self): + """time to submit jobs to dask""" + return self._prepare_dask + @property def conclude(self): """time to conclude""" @@ -136,16 +134,16 @@ class ParallelAnalysisBase(object): class NewAnalysis(ParallelAnalysisBase): def __init__(self, atomgroup, parameter): self._ag = atomgroup - super(NewAnalysis, self).__init__(atomgroup.universe, - self._ag) + self.parameter = parameter + super().__init__(atomgroup.universe) - def _single_frame(self, ts, agroups): + def _single_frame(self): # REQUIRED - # called for every frame. ``ts`` contains the current time step - # and ``agroups`` a tuple of atomgroups that are updated to the - # current frame. Return result of `some_function` for a single - # frame - return some_function(agroups[0], self._parameter) + # called for every frame. It can read all the attributes of + # of itself e.g. current timestep (``self._ts``), atomgroups + # (``self._ag``) that are updated to the current frame and etc. + # Return result of `some_function` for a single frame + return some_function(self._ag, self._parameter) def _conclude(self): # REQUIRED @@ -155,7 +153,7 @@ def _conclude(self): # sensible new variable. self.results = np.concatenate(self._results) # Apply normalisation and averaging to results here if wanted. - self.results /= np.sum(self.results + self.results /= np.sum(self.results) @staticmethod def _reduce(res, result_single_frame): @@ -189,16 +187,13 @@ def _reduce(res, result_single_frame): """ - def __init__(self, universe, atomgroups): + def __init__(self, universe): """Parameters ---------- - Universe : :class:`~MDAnalysis.core.groups.Universe` + universe : :class:`~MDAnalysis.core.groups.Universe` a :class:`MDAnalysis.core.groups.Universe` (the `atomgroups` must belong to this Universe) - atomgroups : tuple of :class:`~MDAnalysis.core.groups.AtomGroup` - atomgroups that are iterated in parallel - Attributes ---------- _results : list @@ -207,10 +202,8 @@ def __init__(self, universe, atomgroups): :meth:`pmda.parallel.ParallelAnalysisBase._single_frame`. """ + self._universe = universe self._trajectory = universe.trajectory - self._top = universe.filename - self._traj = universe.trajectory.filename - self._indices = [ag.indices for ag in atomgroups] @contextmanager def readonly_attributes(self): @@ -232,11 +225,15 @@ def __setattr__(self, key, val): # guards to stop people assigning to self when they shouldn't # if locked, the only attribute you can modify is _attr_lock # if self._attr_lock isn't set, default to unlocked - if key == '_attr_lock' or not getattr(self, '_attr_lock', False): - super(ParallelAnalysisBase, self).__setattr__(key, val) + + # keys that can be changed + if key in ['_ts', 'prepare_dask_total'] or \ + key == '_attr_lock' or \ + not getattr(self, '_attr_lock', False): + super().__setattr__(key, val) else: # raise HalError("I'm sorry Dave, I'm afraid I can't do that") - raise AttributeError("Can't set attribute at this time") + raise AttributeError("Can't set '{}' at this time".format(key)) def _conclude(self): """Finalise the results you've gathered. @@ -253,24 +250,15 @@ def _prepare(self): """additional preparation to run""" pass # pylint: disable=unnecessary-pass - def _single_frame(self, ts, atomgroups): + def _single_frame(self): """Perform computation on a single trajectory frame. Must return computed values as a list. You can only **read** from member variables stored in ``self``. Changing them during - a run will result in undefined behavior. `ts` and any of the + a run will result in undefined behavior. `self._ts` and any of the atomgroups can be changed (but changes will be overwritten when the next time step is read). - Parameters - ---------- - ts : :class:`~MDAnalysis.coordinates.base.Timestep` - The current coordinate frame (time step) in the - trajectory. - atomgroups : tuple - Tuple of :class:`~MDAnalysis.core.groups.AtomGroup` - instances that are updated to the current frame. - Returns ------- values : anything @@ -374,18 +362,15 @@ def run(self, blocks = [] _blocks = [] with self.readonly_attributes(): - for bslice in slices: - task = delayed( - self._dask_helper, pure=False)( - bslice, - self._indices, - self._top, - self._traj, ) - blocks.append(task) - # save the frame numbers for each block - _blocks.append(range(bslice.start, - bslice.stop, bslice.step)) - blocks = delayed(blocks) + with timeit() as prepare_dask: + for bslice in slices: + task = delayed(self._dask_helper, pure=False)(bslice) + blocks.append(task) + # save the frame numbers for each block + _blocks.append(range(bslice.start, + bslice.stop, bslice.step)) + blocks = delayed(blocks) + time_prepare_dask = prepare_dask.elapsed # record the time when scheduler starts working wait_start = time.time() @@ -393,7 +378,7 @@ def run(self, # hack to handle n_frames == 0 in this framework if len(res) == 0: # everything else wants list of block tuples - res = [([], [], [], 0, wait_start, 0, 0)] + res = [([], [], [], wait_start, 0, 0)] # record conclude time with timeit() as conclude: self._results = np.asarray([el[0] for el in res]) @@ -404,23 +389,25 @@ def run(self, self.timing = Timing( np.hstack([el[1] for el in res]), np.hstack([el[2] for el in res]), total.elapsed, - np.array([el[3] for el in res]), time_prepare, + time_prepare, + time_prepare_dask, conclude.elapsed, # waiting time = wait_end - wait_start - np.array([el[4]-wait_start for el in res]), - np.array([el[5] for el in res]), - np.array([el[6] for el in res])) + np.array([el[3] - wait_start for el in res]), + np.array([el[4] for el in res]), + np.array([el[5] for el in res])) + + # To make sure the trajectory is reset to initial state, + # if we are not running the analysis through the whole trajectory. + # With this, we get the same result (state of the trajectory) from + # ParallelAnalysisBase and MDAnalysis.AnalaysisBase. + self._trajectory.rewind() return self - def _dask_helper(self, bslice, indices, top, traj): + def _dask_helper(self, bslice): """helper function to actually setup dask graph""" # wait_end needs to be first line for accurate timing wait_end = time.time() - # record time to generate universe and atom groups - with timeit() as b_universe: - u = mda.Universe(top, traj) - agroups = [u.atoms[idx] for idx in indices] - res = [] times_io = [] times_compute = [] @@ -431,16 +418,16 @@ def _dask_helper(self, bslice, indices, top, traj): with timeit() as b_io: # explicit instead of 'for ts in u.trajectory[bslice]' # so that we can get accurate timing. - ts = u.trajectory[i] + self._ts = self._trajectory[i] # record compute time per frame with timeit() as b_compute: - res = self._reduce(res, self._single_frame(ts, agroups)) + res = self._reduce(res, self._single_frame()) times_io.append(b_io.elapsed) times_compute.append(b_compute.elapsed) # calculate io and compute time per block return np.asarray(res), np.asarray(times_io), np.asarray( - times_compute), b_universe.elapsed, wait_end, np.sum( + times_compute), wait_end, np.sum( times_io), np.sum(times_compute) @staticmethod diff --git a/pmda/rdf.py b/pmda/rdf.py index a8aabde6..2ff887a7 100644 --- a/pmda/rdf.py +++ b/pmda/rdf.py @@ -25,9 +25,6 @@ :inherited-members: """ - -from __future__ import absolute_import, division - import numpy as np from MDAnalysis.lib import distances @@ -93,8 +90,9 @@ class InterRDF(ParallelAnalysisBase): def __init__(self, g1, g2, nbins=75, range=(0.0, 15.0), exclusion_block=None): u = g1.universe - super(InterRDF, self).__init__(u, (g1, g2)) + super().__init__(u) + self._atomgroups = (g1, g2) # collect all atomgroups with the same trajectory object as universe self.nA = len(g1) self.nB = len(g2) @@ -117,22 +115,21 @@ def _prepare(self): # Set the max range to filter the search radius self._maxrange = self.rdf_settings['range'][1] - def _single_frame(self, ts, atomgroups): - g1, g2 = atomgroups - u = g1.universe + def _single_frame(self): + g1, g2 = self._atomgroups pairs, dist = distances.capped_distance(g1.positions, g2.positions, self._maxrange, - box=u.dimensions) + box=self._universe.dimensions) # If provided exclusions, create a mask of _result which # lets us take these out. if self._exclusion_block is not None: - idxA = pairs[:, 0]//self._exclusion_block[0] - idxB = pairs[:, 1]//self._exclusion_block[1] + idxA = pairs[:, 0] // self._exclusion_block[0] + idxB = pairs[:, 1] // self._exclusion_block[1] mask = np.where(idxA != idxB)[0] dist = dist[mask] count = np.histogram(dist, **self.rdf_settings)[0] - volume = u.trajectory.ts.volume + volume = self._universe.trajectory.ts.volume return np.array([count, np.array(volume, dtype=np.float64)]) @@ -150,7 +147,7 @@ def _conclude(self, ): # Volume in each radial shell vol = np.power(self.edges[1:], 3) - np.power(self.edges[:-1], 3) - vol *= 4/3.0 * np.pi + vol *= 4 / 3.0 * np.pi # Average number density box_vol = self.volume / self.n_frames @@ -261,7 +258,9 @@ def __init__(self, u, ags, for pair in ags: atomgroups.append(pair[0]) atomgroups.append(pair[1]) - super(InterRDF_s, self).__init__(u, atomgroups) + super().__init__(u) + + self._atomgroups = atomgroups # List of pairs of AtomGroups self._density = density @@ -289,22 +288,25 @@ 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)] + def _single_frame(self): + ags = [ + [self._atomgroups[2 * i], + self._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) + box=self._universe + .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 + volume = self._universe.trajectory.ts.volume return np.array([np.array(count), np.array(volume, dtype=np.float64)]) @@ -313,7 +315,7 @@ def _conclude(self): 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 + vol *= 4 / 3.0 * np.pi # Empty lists to restore indices, RDF rdf = [] diff --git a/pmda/rms/rmsd.py b/pmda/rms/rmsd.py index c6358f58..340ab1a5 100644 --- a/pmda/rms/rmsd.py +++ b/pmda/rms/rmsd.py @@ -24,8 +24,6 @@ :inherited-members: """ -from __future__ import absolute_import - from MDAnalysis.analysis import rms import numpy as np @@ -127,7 +125,8 @@ class RMSD(ParallelAnalysisBase): """ def __init__(self, mobile, ref, superposition=True): universe = mobile.universe - super(RMSD, self).__init__(universe, (mobile, )) + super().__init__(universe) + self._mobile = mobile self._ref_pos = ref.positions.copy() self.superposition = superposition @@ -137,7 +136,7 @@ def _prepare(self): def _conclude(self): self.rmsd = np.vstack(self._results) - def _single_frame(self, ts, atomgroups): - return (ts.frame, ts.time, - rms.rmsd(atomgroups[0].positions, self._ref_pos, + def _single_frame(self): + return (self._ts.frame, self._ts.time, + rms.rmsd(self._mobile.positions, self._ref_pos, superposition=self.superposition)) diff --git a/pmda/rms/rmsf.py b/pmda/rms/rmsf.py index ad563b55..60ea819b 100644 --- a/pmda/rms/rmsf.py +++ b/pmda/rms/rmsf.py @@ -24,9 +24,6 @@ MDAnalysis.analysis.rms.RMSF """ - -from __future__ import absolute_import, division - import numpy as np from pmda.parallel import ParallelAnalysisBase @@ -172,14 +169,12 @@ class RMSF(ParallelAnalysisBase): def __init__(self, atomgroup): u = atomgroup.universe - super(RMSF, self).__init__(u, (atomgroup, )) + super().__init__(u) self._atomgroup = atomgroup - self._top = u.filename - self._traj = u.trajectory.filename - def _single_frame(self, ts, atomgroups): + def _single_frame(self): # mean and sum of squares calculations done in _reduce() - return atomgroups[0] + return self._atomgroup def _conclude(self): """ diff --git a/pmda/test/test_hydrogenbonds_analysis.py b/pmda/test/test_hydrogenbonds_analysis.py index 2917b523..e7851175 100644 --- a/pmda/test/test_hydrogenbonds_analysis.py +++ b/pmda/test/test_hydrogenbonds_analysis.py @@ -104,7 +104,7 @@ def test_count_by_ids(self, h): def test_universe(self, h, universe): ref = universe.atoms.positions h.run(n_jobs=4, n_blocks=4) - u = h._universe() + u = h._universe assert_array_almost_equal(u.atoms.positions, ref) @@ -269,5 +269,5 @@ def test_universe(self): universe = MDAnalysis.Universe(GRO) ref = universe.atoms.positions h = HydrogenBondAnalysis(universe, **self.kwargs) - u = h._universe() + u = h._universe assert_array_almost_equal(u.atoms.positions, ref) diff --git a/pmda/test/test_parallel.py b/pmda/test/test_parallel.py index 0cefe72b..789c0c6b 100644 --- a/pmda/test/test_parallel.py +++ b/pmda/test/test_parallel.py @@ -24,23 +24,23 @@ def test_timeing(): io = np.arange(5) compute = np.arange(5) + 1 total = 5 - universe = np.arange(2) prepare = 3 + prepare_dask = 4 conclude = 6 wait = 12 io_block = np.sum(io) compute_block = np.sum(compute) timing = parallel.Timing(io, compute, total, - universe, prepare, conclude, wait, + prepare, prepare_dask, conclude, wait, io_block, compute_block,) np.testing.assert_equal(timing.io, io) np.testing.assert_equal(timing.compute, compute) np.testing.assert_equal(timing.total, total) - np.testing.assert_equal(timing.universe, universe) np.testing.assert_equal(timing.cumulate_time, np.sum(io) + np.sum(compute)) np.testing.assert_equal(timing.prepare, prepare) + np.testing.assert_equal(timing.prepare_dask, prepare_dask) np.testing.assert_equal(timing.conclude, conclude) np.testing.assert_equal(timing.wait, wait) np.testing.assert_equal(timing.io_block, io_block) @@ -50,7 +50,8 @@ def test_timeing(): class NoneAnalysis(parallel.ParallelAnalysisBase): def __init__(self, atomgroup): universe = atomgroup.universe - super(NoneAnalysis, self).__init__(universe, (atomgroup, )) + super().__init__(universe) + self._atomgroup = atomgroup def _prepare(self): pass @@ -58,8 +59,8 @@ def _prepare(self): def _conclude(self): self.res = np.concatenate(self._results) - def _single_frame(self, ts, atomgroups): - return ts.frame + def _single_frame(self): + return self._ts.frame @pytest.fixture @@ -72,7 +73,7 @@ def analysis(): @pytest.mark.parametrize('n_jobs', (1, 2)) def test_all_frames(analysis, n_jobs): analysis.run(n_jobs=n_jobs) - u = mda.Universe(analysis._top, analysis._traj) + u = analysis._universe assert len(analysis.res) == u.trajectory.n_frames @@ -84,7 +85,7 @@ def test_sub_frames(analysis, n_jobs): @pytest.mark.parametrize('n_jobs', (1, 2, 3)) def test_no_frames(analysis, n_jobs): - u = mda.Universe(analysis._top, analysis._traj) + u = analysis._universe n_frames = u.trajectory.n_frames with pytest.warns(UserWarning): analysis.run(start=n_frames, stop=n_frames+1, n_jobs=n_jobs) @@ -95,7 +96,6 @@ def test_no_frames(analysis, n_jobs): np.testing.assert_equal(analysis.timing.io_block, [0]) np.testing.assert_equal(analysis.timing.compute_block, [0]) np.testing.assert_equal(analysis.timing.wait, [0]) - assert analysis.timing.universe == 0 def test_scheduler(analysis, scheduler): @@ -103,7 +103,7 @@ def test_scheduler(analysis, scheduler): def test_nframes_less_nblocks_warning(analysis): - u = mda.Universe(analysis._top, analysis._traj) + u = analysis._universe n_frames = u.trajectory.n_frames with pytest.raises(ValueError): analysis.run(stop=2, n_blocks=4, n_jobs=2) @@ -124,7 +124,7 @@ def test_guess_nblocks(analysis): @pytest.mark.parametrize('n_blocks', np.arange(1, 11)) def test_blocks(analysis, n_blocks): analysis.run(n_blocks=n_blocks) - u = mda.Universe(analysis._top, analysis._traj) + u = analysis._universe n_frames = u.trajectory.n_frames start, stop, step = u.trajectory.check_slice_indices( None, None, None) @@ -137,7 +137,7 @@ def test_blocks(analysis, n_blocks): def test_attrlock(): u = mda.Universe(PSF, DCD) - pab = parallel.ParallelAnalysisBase(u, (u.atoms,)) + pab = parallel.ParallelAnalysisBase(u) # Should initially be allowed to set attributes pab.thing1 = 24 diff --git a/setup.cfg b/setup.cfg index b602da63..a835515d 100644 --- a/setup.cfg +++ b/setup.cfg @@ -15,3 +15,4 @@ versionfile_source = pmda/_version.py versionfile_build = pmda/_version.py tag_prefix = parentdir_prefix = pmda + diff --git a/setup.py b/setup.py index 0d465d87..54cd9810 100644 --- a/setup.py +++ b/setup.py @@ -8,8 +8,16 @@ # Released under the GNU Public Licence, v2 or any higher version from setuptools import setup, find_packages +import sys import versioneer +# Make sure I have the right Python version. +if sys.version_info[:2] < (3, 6): + print('PMDA requires Python 3.6 or better. Python {0:d}.{1:d} detected'.format(* + sys.version_info[:2])) + print('Please upgrade your version of Python.') + sys.exit(-1) + with open('README.rst', 'r') as f: long_description = f.read() @@ -33,8 +41,6 @@ 'Programming Language :: Python', 'Topic :: Scientific/Engineering', 'Topic :: Software Development :: Libraries :: Python Modules', - 'Programming Language :: Python :: 2', - 'Programming Language :: Python :: 2.7', 'Programming Language :: Python :: 3', 'Programming Language :: Python :: 3.6', 'Programming Language :: Python :: 3.7', @@ -48,7 +54,8 @@ }, packages=find_packages(), install_requires=[ - 'MDAnalysis>=1.0.0', + # 'MDAnalysis>=2.0.0', + 'mdanalysis @ git+https://github.com/MDAnalysis/mdanalysis#egg=mdanalysis&subdirectory=package', 'dask>=0.18', 'distributed', 'six', @@ -58,5 +65,6 @@ ], tests_require=[ 'pytest', - 'MDAnalysisTests>=1.0.0', # keep + # 'MDAnalysisTests>=2.0.0', # keep + 'mdanalysistests @ git+https://github.com/MDAnalysis/mdanalysis#egg=mdanalysistests&subdirectory=testsuite' ], )