From 4e3243dd3724e9455649f9fcb69f59ef1481ab0f Mon Sep 17 00:00:00 2001 From: Lily Wang Date: Thu, 6 May 2021 16:42:36 -0700 Subject: [PATCH 1/2] add results --- package/CHANGELOG | 2 + package/MDAnalysis/analysis/helix_analysis.py | 58 ++++++++++--------- .../analysis/test_helix_analysis.py | 36 ++++++------ 3 files changed, 50 insertions(+), 46 deletions(-) diff --git a/package/CHANGELOG b/package/CHANGELOG index abaf73f9f2a..80b613196df 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -175,6 +175,8 @@ Enhancements checking if it can be used in parallel analysis. (Issue #2996, PR #2950) Changes + * `helix_analysis.HELANAL` now uses the `analysis.base.Results` class to + store results attributes (Issue #3261, #3267) * `analysis.align.AlignTraj` and `analysis.align.AverageStructure` now store their result attributes using the `analysis.base.Results` class (Issues #3278 and #3261) diff --git a/package/MDAnalysis/analysis/helix_analysis.py b/package/MDAnalysis/analysis/helix_analysis.py index 5e8c516ab7e..f7eca993d8b 100644 --- a/package/MDAnalysis/analysis/helix_analysis.py +++ b/package/MDAnalysis/analysis/helix_analysis.py @@ -29,7 +29,7 @@ :Year: 2020 :Copyright: GNU Public License v3 -.. versionadded:: 1.0.0 +.. versionadded:: 2.0.0 This module contains code to analyse protein helices using the HELANAL_ algorithm @@ -63,7 +63,9 @@ helanal = hel.HELANAL(u, select='name CA and resnum 161-187') helanal.run() - print(helanal.summary) +All computed properties are available in ``.results``:: + + print(helanal.results.summary) Alternatively, you can analyse several helices at once by passing in multiple selection strings:: @@ -309,42 +311,42 @@ class HELANAL(AnalysisBase): Attributes ---------- - local_twists : array or list of arrays + results.local_twists : array or list of arrays The local twist angle from atom i+1 to i+2. Each array has shape (n_frames, n_residues-3) - local_nres_per_turn : array or list of arrays + results.local_nres_per_turn : array or list of arrays Number of residues per turn, based on local_twist. Each array has shape (n_frames, n_residues-3) - local_axes : array or list of arrays + results.local_axes : array or list of arrays The length-wise helix axis of the local window. Each array has shape (n_frames, n_residues-3, 3) - local_heights : array or list of arrays + results.local_heights : array or list of arrays The rise of each local helix. Each array has shape (n_frames, n_residues-3) - local_helix_directions : array or list of arrays + results.local_helix_directions : array or list of arrays The unit vector from each local origin to atom i+1. Each array has shape (n_frames, n_residues-2, 3) - local_origins :array or list of arrays + results.local_origins :array or list of arrays The projected origin for each helix. Each array has shape (n_frames, n_residues-2, 3) - local_screw_angles : array or list of arrays + results.local_screw_angles : array or list of arrays The local screw angle for each helix. Each array has shape (n_frames, n_residues-2) - local_bends : array or list of arrays + results.local_bends : array or list of arrays The angles between local helix axes, 3 windows apart. Each array has shape (n_frames, n_residues-6) - all_bends : array or list of arrays + results.all_bends : array or list of arrays The angles between local helix axes. Each array has shape (n_frames, n_residues-3, n_residues-3) - global_axis : array or list of arrays + results.global_axis : array or list of arrays The length-wise axis for the overall helix. This points at the first helix window in the helix, so it runs opposite to the direction of the residue numbers. Each array has shape (n_frames, 3) - global_tilts : array or list of arrays + results.global_tilts : array or list of arrays The angle between the global axis and the reference axis. Each array has shape (n_frames,) - summary : dict or list of dicts + results.summary : dict or list of dicts Summary of stats for each property: the mean, the sample standard deviation, and the mean absolute deviation. """ @@ -418,49 +420,49 @@ def _prepare(self): for key, dims in self.attr_shapes.items(): empty = [self._zeros_per_frame( dims, n_positions=n) for n in n_res] - setattr(self, key, empty) + setattr(self.results, key, empty) - self.global_axis = [self._zeros_per_frame((3,)) for n in n_res] - self.all_bends = [self._zeros_per_frame((n-3, n-3)) for n in n_res] + self.results.global_axis = [self._zeros_per_frame((3,)) for n in n_res] + self.results.all_bends = [self._zeros_per_frame((n-3, n-3)) for n in n_res] def _single_frame(self): _f = self._frame_index for i, ag in enumerate(self.atomgroups): results = helix_analysis(ag.positions, ref_axis=self.ref_axis) for key, value in results.items(): - attr = getattr(self, key) + attr = getattr(self.results, key) attr[i][_f] = value def _conclude(self): # compute tilt of global axes - self.global_tilts = [] + self.results.global_tilts = tilts = [] norm_ref = (self.ref_axis**2).sum() ** 0.5 - for axes in self.global_axis: + for axes in self.results.global_axis: cos = np.matmul(self.ref_axis, axes.T) / \ (mdamath.pnorm(axes)*norm_ref) cos = np.clip(cos, -1.0, 1.0) - self.global_tilts.append(np.rad2deg(np.arccos(cos))) + tilts.append(np.rad2deg(np.arccos(cos))) global_attrs = ['global_axis', 'global_tilts', 'all_bends'] attrnames = list(self.attr_shapes.keys()) + global_attrs # summarise - self.summary = [] + self.results.summary = [] for i in range(len(self.atomgroups)): stats = {} for name in attrnames: - attr = getattr(self, name) + attr = getattr(self.results, name) mean = attr[i].mean(axis=0) dev = np.abs(attr[i]-mean) stats[name] = {'mean': mean, 'sample_sd': attr[i].std(axis=0, ddof=1), 'abs_dev': dev.mean(axis=0)} - self.summary.append(stats) + self.results.summary.append(stats) # flatten? if len(self.atomgroups) == 1 and self._flatten: for name in attrnames + ['summary']: - attr = getattr(self, name) - setattr(self, name, attr[0]) + attr = getattr(self.results, name) + setattr(self.results, name, attr[0]) def universe_from_origins(self): """ @@ -471,7 +473,7 @@ def universe_from_origins(self): Universe or list of Universes """ try: - origins = self.local_origins + origins = self.results.local_origins except AttributeError: raise ValueError('Call run() before universe_from_origins') @@ -485,6 +487,6 @@ def universe_from_origins(self): atom_resindex=np.arange(n_res), trajectory=True).load_new(xyz) universe.append(u) - if not isinstance(self.local_origins, list): + if not isinstance(self.results.local_origins, list): universe = universe[0] return universe diff --git a/testsuite/MDAnalysisTests/analysis/test_helix_analysis.py b/testsuite/MDAnalysisTests/analysis/test_helix_analysis.py index f96691d5292..57bd538fa23 100644 --- a/testsuite/MDAnalysisTests/analysis/test_helix_analysis.py +++ b/testsuite/MDAnalysisTests/analysis/test_helix_analysis.py @@ -351,7 +351,7 @@ def helanal(self, psf_ca): return ha.run(start=10, stop=80) def test_regression_summary(self, helanal): - bends = helanal.summary['all_bends'] + bends = helanal.results.summary['all_bends'] old_helanal = read_bending_matrix(HELANAL_BENDING_MATRIX_SUBSET) assert_almost_equal(np.triu(bends['mean'], 1), old_helanal['Mean'], decimal=1) @@ -368,11 +368,11 @@ def test_regression_values(self): for key, value in HELANAL_SINGLE_DATA.items(): if 'summary' in key: - data = getattr(ha, key.split()[0]) + data = getattr(ha.results, key.split()[0]) calculated = [data.mean(), data.std(ddof=1), np.fabs(data-data.mean()).mean()] else: - calculated = getattr(ha, key)[0] + calculated = getattr(ha.results, key)[0] msg = 'Mismatch between calculated and reference {}' assert_almost_equal(calculated, value, @@ -385,10 +385,10 @@ def test_multiple_selections(self, psf_ca): ha.run() n_frames = len(psf_ca.universe.trajectory) assert len(ha.atomgroups) == 2 - assert len(ha.summary) == 2 - assert len(ha.all_bends) == 2 - assert ha.all_bends[0].shape == (n_frames, 8, 8) - assert ha.all_bends[1].shape == (n_frames, 18, 18) + assert len(ha.results.summary) == 2 + assert len(ha.results.all_bends) == 2 + assert ha.results.all_bends[0].shape == (n_frames, 8, 8) + assert ha.results.all_bends[1].shape == (n_frames, 18, 18) def test_universe_from_origins(self, helanal): u = helanal.universe_from_origins() @@ -464,22 +464,22 @@ def test_len_groups_short(self, psf_ca): def test_helanal_zigzag(self, zigzag, ref_axis, screw_angles): ha = hel.HELANAL(zigzag, select="all", ref_axis=ref_axis, flatten_single_helix=True).run() - assert_almost_equal(ha.local_twists, 180, decimal=4) - assert_almost_equal(ha.local_nres_per_turn, 2, decimal=4) - assert_almost_equal(ha.global_axis, [[0, 0, -1]], decimal=4) + assert_almost_equal(ha.results.local_twists, 180, decimal=4) + assert_almost_equal(ha.results.local_nres_per_turn, 2, decimal=4) + assert_almost_equal(ha.results.global_axis, [[0, 0, -1]], decimal=4) # all 0 vectors - assert_almost_equal(ha.local_axes, 0, decimal=4) - assert_almost_equal(ha.local_bends, 0, decimal=4) - assert_almost_equal(ha.all_bends, 0, decimal=4) - assert_almost_equal(ha.local_heights, 0, decimal=4) - assert_almost_equal(ha.local_helix_directions[0][0::2], + assert_almost_equal(ha.results.local_axes, 0, decimal=4) + assert_almost_equal(ha.results.local_bends, 0, decimal=4) + assert_almost_equal(ha.results.all_bends, 0, decimal=4) + assert_almost_equal(ha.results.local_heights, 0, decimal=4) + assert_almost_equal(ha.results.local_helix_directions[0][0::2], [[-1, 0, 0]]*49, decimal=4) - assert_almost_equal(ha.local_helix_directions[0][1::2], + assert_almost_equal(ha.results.local_helix_directions[0][1::2], [[1, 0, 0]]*49, decimal=4) origins = zigzag.atoms.positions[1:-1].copy() origins[:, 0] = 0 - assert_almost_equal(ha.local_origins[0], origins, decimal=4) - assert_almost_equal(ha.local_screw_angles[0], + assert_almost_equal(ha.results.local_origins[0], origins, decimal=4) + assert_almost_equal(ha.results.local_screw_angles[0], screw_angles*49, decimal=4) From d63d3899ad8626a877ceae25549273f3db2bc9a0 Mon Sep 17 00:00:00 2001 From: Lily Wang <31115101+lilyminium@users.noreply.github.com> Date: Thu, 6 May 2021 17:23:08 -0700 Subject: [PATCH 2/2] Apply suggestions from code review Co-authored-by: Oliver Beckstein --- package/MDAnalysis/analysis/helix_analysis.py | 10 +++++----- .../MDAnalysisTests/analysis/test_helix_analysis.py | 4 ++-- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/package/MDAnalysis/analysis/helix_analysis.py b/package/MDAnalysis/analysis/helix_analysis.py index f7eca993d8b..7518ee6c376 100644 --- a/package/MDAnalysis/analysis/helix_analysis.py +++ b/package/MDAnalysis/analysis/helix_analysis.py @@ -420,7 +420,7 @@ def _prepare(self): for key, dims in self.attr_shapes.items(): empty = [self._zeros_per_frame( dims, n_positions=n) for n in n_res] - setattr(self.results, key, empty) + self.results[key] = empty self.results.global_axis = [self._zeros_per_frame((3,)) for n in n_res] self.results.all_bends = [self._zeros_per_frame((n-3, n-3)) for n in n_res] @@ -430,7 +430,7 @@ def _single_frame(self): for i, ag in enumerate(self.atomgroups): results = helix_analysis(ag.positions, ref_axis=self.ref_axis) for key, value in results.items(): - attr = getattr(self.results, key) + attr = self.results[key] attr[i][_f] = value def _conclude(self): @@ -450,7 +450,7 @@ def _conclude(self): for i in range(len(self.atomgroups)): stats = {} for name in attrnames: - attr = getattr(self.results, name) + attr = self.results[name] mean = attr[i].mean(axis=0) dev = np.abs(attr[i]-mean) stats[name] = {'mean': mean, @@ -461,8 +461,8 @@ def _conclude(self): # flatten? if len(self.atomgroups) == 1 and self._flatten: for name in attrnames + ['summary']: - attr = getattr(self.results, name) - setattr(self.results, name, attr[0]) + attr = self.results[name] + self.results[name] = attr[0] def universe_from_origins(self): """ diff --git a/testsuite/MDAnalysisTests/analysis/test_helix_analysis.py b/testsuite/MDAnalysisTests/analysis/test_helix_analysis.py index 57bd538fa23..1e39147be81 100644 --- a/testsuite/MDAnalysisTests/analysis/test_helix_analysis.py +++ b/testsuite/MDAnalysisTests/analysis/test_helix_analysis.py @@ -368,11 +368,11 @@ def test_regression_values(self): for key, value in HELANAL_SINGLE_DATA.items(): if 'summary' in key: - data = getattr(ha.results, key.split()[0]) + data = ha.results[key.split()[0]] calculated = [data.mean(), data.std(ddof=1), np.fabs(data-data.mean()).mean()] else: - calculated = getattr(ha.results, key)[0] + calculated = ha.results[key][0] msg = 'Mismatch between calculated and reference {}' assert_almost_equal(calculated, value,