From db2355057abb96076ba9991fe59c41dfb364ab4a Mon Sep 17 00:00:00 2001 From: g2707 <55765804+shaivimalik@users.noreply.github.com> Date: Fri, 24 Mar 2023 02:02:28 +0530 Subject: [PATCH 1/9] performance improvements --- package/MDAnalysis/core/groups.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/package/MDAnalysis/core/groups.py b/package/MDAnalysis/core/groups.py index 67fee225638..06156d97c87 100644 --- a/package/MDAnalysis/core/groups.py +++ b/package/MDAnalysis/core/groups.py @@ -1071,7 +1071,7 @@ def center(self, weights, wrap=False, unwrap=False, compound='group'): return coords.mean(axis=0) # promote weights to dtype if required: weights = weights.astype(dtype, copy=False) - return (coords * weights[:, None]).sum(axis=0) / weights.sum() + return np.einsum('ij,ij->j',coords,weights[:, None]) / weights.sum() # When compound split caching gets implemented it will be clever to # preempt at this point whether or not stable sorting will be needed @@ -1100,7 +1100,7 @@ def center(self, weights, wrap=False, unwrap=False, compound='group'): _centers = _coords.mean(axis=1) else: _weights = weights[atom_mask] - _centers = (_coords * _weights[:, :, None]).sum(axis=1) + _centers = np.einsum('ijk,ijk->ik',_coords,_weights[:, :, None]) _centers /= _weights.sum(axis=1)[:, None] centers[compound_mask] = _centers if wrap: @@ -1864,7 +1864,7 @@ def unwrap(self, compound='fragments', reference='com', inplace=True): raise ValueError("Cannot perform unwrap with " "reference='com' because the total " "mass of the group is zero.") - refpos = np.sum(positions * masses[:, None], axis=0) + refpos = np.einsum('ij,ij->j',positions,masses[:, None]) refpos /= total_mass else: # reference == 'cog' refpos = positions.mean(axis=0) @@ -1893,8 +1893,8 @@ def unwrap(self, compound='fragments', reference='com', inplace=True): "reference='com' because the " "total mass of at least one of " "the {} is zero.".format(comp)) - refpos = np.sum(positions[atom_mask] - * masses[:, :, None], axis=1) + refpos = np.einsum('ijk,ijk->ik',positions[atom_mask], + masses[:, :, None]) refpos /= total_mass[:, None] else: # reference == 'cog' refpos = positions[atom_mask].mean(axis=1) From 8986ffe4d232e33d4ad185abb16b7200aedc26d0 Mon Sep 17 00:00:00 2001 From: g2707 <55765804+shaivimalik@users.noreply.github.com> Date: Tue, 28 Mar 2023 03:37:22 +0530 Subject: [PATCH 2/9] numpy changes --- package/MDAnalysis/analysis/bat.py | 21 +++++++++++---------- package/MDAnalysis/core/topologyattrs.py | 21 ++++++++++----------- package/MDAnalysis/lib/transformations.py | 12 +++++------- 3 files changed, 26 insertions(+), 28 deletions(-) diff --git a/package/MDAnalysis/analysis/bat.py b/package/MDAnalysis/analysis/bat.py index 0e1e19f5b50..9918162f386 100644 --- a/package/MDAnalysis/analysis/bat.py +++ b/package/MDAnalysis/analysis/bat.py @@ -372,13 +372,14 @@ def _single_frame(self): v01 = p1 - p0 v21 = p1 - p2 # Internal coordinates - r01 = np.sqrt(np.sum(v01 * - v01)) # Distance between first two root atoms - r12 = np.sqrt(np.sum(v21 * - v21)) # Distance between second two root atoms + r01 = np.sqrt(np.einsum('i,i->',v01,v01)) + # Distance between first two root atoms + r12 = np.sqrt(np.einsum('i,i->',v21,v21)) + # Distance between second two root atoms # Angle between root atoms - a012 = np.arccos(max(-1.,min(1.,np.sum(v01*v21)/\ - np.sqrt(np.sum(v01*v01)*np.sum(v21*v21))))) + a012 = np.arccos(max(-1.,min(1.,np.einsum('i,i->',v01,v21)/\ + np.sqrt(np.einsum('i,i->',v01,v01)* + np.einsum('i,i->',v21,v21))))) # External coordinates e = v01 / r01 phi = np.arctan2(e[1], e[0]) # Polar angle @@ -545,12 +546,12 @@ def Cartesian(self, bat_frame): cs_tor = np.cos(torsion) v21 = p1 - p2 - v21 /= np.sqrt(np.sum(v21 * v21)) + v21 /= np.sqrt(np.einsum('i,i->',v21,v21)) v32 = p2 - p3 - v32 /= np.sqrt(np.sum(v32 * v32)) + v32 /= np.sqrt(np.einsum('i,i->',v32,v32)) vp = np.cross(v32, v21) - cs = np.sum(v21 * v32) + cs = np.einsum('i,i->',v21,v32) sn = max(np.sqrt(1.0 - cs * cs), 0.0000000001) vp = vp / sn @@ -563,4 +564,4 @@ def Cartesian(self, bat_frame): @property def atoms(self): """The atomgroup for which BAT are computed (read-only property)""" - return self._ag + return self._ag \ No newline at end of file diff --git a/package/MDAnalysis/core/topologyattrs.py b/package/MDAnalysis/core/topologyattrs.py index e06cc63c3d1..5ff9a18fccc 100644 --- a/package/MDAnalysis/core/topologyattrs.py +++ b/package/MDAnalysis/core/topologyattrs.py @@ -1718,8 +1718,8 @@ def radius_of_gyration(group, wrap=False, **kwargs): else: recenteredpos = atomgroup.positions - com - rog_sq = np.sum(masses * np.sum(recenteredpos**2, - axis=1)) / atomgroup.total_mass() + rog_sq = np.einsum('ij,ij->',masses,np.einsum('ijk,ijk->ik', + recenteredpos,recenteredpos))/atomgroup.total_mass() return np.sqrt(rog_sq) @@ -2232,8 +2232,8 @@ def dipole_vector(group, ) - ref) else: recenteredpos = (atomgroup.positions - ref) - dipole_vector = np.sum(recenteredpos * charges[:, np.newaxis], - axis=0) + dipole_vector = np.einsum('ijk,ijk->jk',recenteredpos, + charges[:, np.newaxis]) else: (atom_masks, compound_masks, n_compounds) = atomgroup._split_by_compound_indices(compound) @@ -2248,10 +2248,9 @@ def dipole_vector(group, dipole_vector = np.empty((n_compounds, 3), dtype=np.float64) for compound_mask, atom_mask in zip(compound_masks, atom_masks): - dipole_vector[compound_mask] = np.sum( - (coords[atom_mask] - ref[compound_mask][:, None, :]) * - chgs[atom_mask][:, :, None], - axis=1) + dipole_vector[compound_mask] = np.einsum('ijk,ijk->ik', + (coords[atom_mask] - ref[compound_mask][:, None, :]), + chgs[atom_mask][:, :, None]) return dipole_vector @@ -2315,9 +2314,9 @@ def dipole_moment(group, **kwargs): dipole_vector = atomgroup.dipole_vector(**kwargs) if len(dipole_vector.shape) > 1: - dipole_moment = np.sqrt(np.sum(dipole_vector**2, axis=1)) + dipole_moment = np.sqrt(np.einsum('ij,ij->i',dipole_vector,dipole_vector)) else: - dipole_moment = np.sqrt(np.sum(dipole_vector**2)) + dipole_moment = np.sqrt(np.einsum('i,i->',dipole_vector,dipole_vector)) return dipole_moment @@ -3305,4 +3304,4 @@ class CMaps(_Connection): attrname = 'cmaps' singular = 'cmaps' transplants = defaultdict(list) - _n_atoms = 5 + _n_atoms = 5 \ No newline at end of file diff --git a/package/MDAnalysis/lib/transformations.py b/package/MDAnalysis/lib/transformations.py index 5f328e918f0..8509ee17627 100644 --- a/package/MDAnalysis/lib/transformations.py +++ b/package/MDAnalysis/lib/transformations.py @@ -941,9 +941,9 @@ def superimposition_matrix(v0, v1, scaling=False, usesvd=True): M[:3, :3] = R else: # compute symmetric matrix N - xx, yy, zz = np.sum(v0 * v1, axis=1) - xy, yz, zx = np.sum(v0 * np.roll(v1, -1, axis=0), axis=1) - xz, yx, zy = np.sum(v0 * np.roll(v1, -2, axis=0), axis=1) + xx, yy, zz = np.einsum('ij,ij->i', v0 , v1) + xy, yz, zx = np.einsum('ij,ij->i', v0, np.roll(v1, -1, axis=0)) + xz, yx, zy = np.einsum('ij,ij->i', v0, np.roll(v1, -2, axis=0)) N = ( (xx + yy + zz, 0.0, 0.0, 0.0), (yz - zy, xx - yy - zz, 0.0, 0.0), @@ -958,9 +958,7 @@ def superimposition_matrix(v0, v1, scaling=False, usesvd=True): # scale: ratio of rms deviations from centroid if scaling: - v0 *= v0 - v1 *= v1 - M[:3, :3] *= math.sqrt(np.sum(v1) / np.sum(v0)) + M[:3, :3] *= math.sqrt(np.einsum('ij,ij->',v1,v1) / np.einsum('ij,ij->',v0,v0)) # translation M[:3, 3] = t1 @@ -1827,4 +1825,4 @@ def rotaxis(a, b): _import_module('_transformations') # Documentation in HTML format can be generated with Epydoc -__docformat__ = "restructuredtext en" +__docformat__ = "restructuredtext en" \ No newline at end of file From 1a7d5fe018d13f8e4701e81c47c6e5767e144abc Mon Sep 17 00:00:00 2001 From: g2707 <55765804+shaivimalik@users.noreply.github.com> Date: Tue, 28 Mar 2023 04:04:23 +0530 Subject: [PATCH 3/9] fixes --- package/MDAnalysis/core/topologyattrs.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/package/MDAnalysis/core/topologyattrs.py b/package/MDAnalysis/core/topologyattrs.py index 5ff9a18fccc..f23d2315990 100644 --- a/package/MDAnalysis/core/topologyattrs.py +++ b/package/MDAnalysis/core/topologyattrs.py @@ -1718,7 +1718,7 @@ def radius_of_gyration(group, wrap=False, **kwargs): else: recenteredpos = atomgroup.positions - com - rog_sq = np.einsum('ij,ij->',masses,np.einsum('ijk,ijk->ik', + rog_sq = np.einsum('ij,ij->',masses,np.einsum('ij,ij->i', recenteredpos,recenteredpos))/atomgroup.total_mass() return np.sqrt(rog_sq) @@ -2232,7 +2232,7 @@ def dipole_vector(group, ) - ref) else: recenteredpos = (atomgroup.positions - ref) - dipole_vector = np.einsum('ijk,ijk->jk',recenteredpos, + dipole_vector = np.einsum('ij,ij->j',recenteredpos, charges[:, np.newaxis]) else: (atom_masks, compound_masks, From c50024bd552c50c7c657bbc3d710291d2bd47259 Mon Sep 17 00:00:00 2001 From: g2707 <55765804+shaivimalik@users.noreply.github.com> Date: Tue, 28 Mar 2023 04:25:53 +0530 Subject: [PATCH 4/9] minor fixes --- package/MDAnalysis/core/topologyattrs.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/package/MDAnalysis/core/topologyattrs.py b/package/MDAnalysis/core/topologyattrs.py index f23d2315990..a41ae176326 100644 --- a/package/MDAnalysis/core/topologyattrs.py +++ b/package/MDAnalysis/core/topologyattrs.py @@ -1718,7 +1718,7 @@ def radius_of_gyration(group, wrap=False, **kwargs): else: recenteredpos = atomgroup.positions - com - rog_sq = np.einsum('ij,ij->',masses,np.einsum('ij,ij->i', + rog_sq = np.einsum('i,i->',masses,np.einsum('ij,ij->i', recenteredpos,recenteredpos))/atomgroup.total_mass() return np.sqrt(rog_sq) From 3e5dd02e17b344cc4d09c51f67ea7967469e8874 Mon Sep 17 00:00:00 2001 From: g2707 <55765804+shaivimalik@users.noreply.github.com> Date: Thu, 30 Mar 2023 07:36:25 +0530 Subject: [PATCH 5/9] final changes --- package/AUTHORS | 1 + package/CHANGELOG | 5 ++++- package/MDAnalysis/analysis/bat.py | 3 +-- package/MDAnalysis/core/topologyattrs.py | 5 +++-- package/MDAnalysis/lib/transformations.py | 4 +++- 5 files changed, 12 insertions(+), 6 deletions(-) diff --git a/package/AUTHORS b/package/AUTHORS index e4d426d89c0..efefebd9555 100644 --- a/package/AUTHORS +++ b/package/AUTHORS @@ -203,6 +203,7 @@ Chronological list of authors - Patricio Barletta - Mikhail Glagolev 2023 + - Shaivi Malik - Christian Pfaendner - Pratham Chauhan - Meet Brijwani diff --git a/package/CHANGELOG b/package/CHANGELOG index 0e3f883848b..5f60d595896 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -16,7 +16,8 @@ The rules for this file: ??/??/?? IAlibay, pgbarletta, mglagolev, hmacdope, manuel.nuno.melo, chrispfae, ooprathamm, MeetB7, BFedder, v-parmar, MoSchaeffler, jbarnoud, jandom, - xhgchen, jaclark5, DrDomenicoMarson, AHMED-salah00, schlaicha, SophiaRuan + xhgchen, jaclark5, DrDomenicoMarson, AHMED-salah00, schlaicha, SophiaRuan, + g2707 * 2.5.0 Fixes @@ -50,6 +51,8 @@ Enhancements and SegmentGroup. (PR #3953) Changes + * einsum method for Einstein summation convention introduced to increase + performance (Issue #2063, PR #4089) * Add progress bars to track the progress of _conclude() functions (_conclude_simple() and _conclude_fft()) in msd.py (Issue #4070, PR #4072) * As per NEP29 the minimum supported NumPy version has been raised to 1.21 diff --git a/package/MDAnalysis/analysis/bat.py b/package/MDAnalysis/analysis/bat.py index 9918162f386..6d7e3cd9e57 100644 --- a/package/MDAnalysis/analysis/bat.py +++ b/package/MDAnalysis/analysis/bat.py @@ -378,8 +378,7 @@ def _single_frame(self): # Distance between second two root atoms # Angle between root atoms a012 = np.arccos(max(-1.,min(1.,np.einsum('i,i->',v01,v21)/\ - np.sqrt(np.einsum('i,i->',v01,v01)* - np.einsum('i,i->',v21,v21))))) + np.sqrt(np.einsum('i,i->',v01,v01)*np.einsum('i,i->',v21,v21))))) # External coordinates e = v01 / r01 phi = np.arctan2(e[1], e[0]) # Polar angle diff --git a/package/MDAnalysis/core/topologyattrs.py b/package/MDAnalysis/core/topologyattrs.py index a41ae176326..6796c56dd43 100644 --- a/package/MDAnalysis/core/topologyattrs.py +++ b/package/MDAnalysis/core/topologyattrs.py @@ -2249,8 +2249,9 @@ def dipole_vector(group, dipole_vector = np.empty((n_compounds, 3), dtype=np.float64) for compound_mask, atom_mask in zip(compound_masks, atom_masks): dipole_vector[compound_mask] = np.einsum('ijk,ijk->ik', - (coords[atom_mask] - ref[compound_mask][:, None, :]), - chgs[atom_mask][:, :, None]) + (coords[atom_mask]- + ref[compound_mask][:, None, :]), + chgs[atom_mask][:, :, None]) return dipole_vector diff --git a/package/MDAnalysis/lib/transformations.py b/package/MDAnalysis/lib/transformations.py index 8509ee17627..b6e88e18362 100644 --- a/package/MDAnalysis/lib/transformations.py +++ b/package/MDAnalysis/lib/transformations.py @@ -958,7 +958,9 @@ def superimposition_matrix(v0, v1, scaling=False, usesvd=True): # scale: ratio of rms deviations from centroid if scaling: - M[:3, :3] *= math.sqrt(np.einsum('ij,ij->',v1,v1) / np.einsum('ij,ij->',v0,v0)) + M[:3, :3] *= math.sqrt(np.einsum('ij,ij->',v1,v1) / + np.einsum('ij,ij->',v0,v0)) + # translation M[:3, 3] = t1 From 4b3b7b61e2bf083d314ccd1886d8d4ec7d7f9ece Mon Sep 17 00:00:00 2001 From: g2707 <128649094+g2707@users.noreply.github.com> Date: Thu, 30 Mar 2023 08:42:57 +0530 Subject: [PATCH 6/9] Update AUTHORS --- package/AUTHORS | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/package/AUTHORS b/package/AUTHORS index efefebd9555..5e8122153db 100644 --- a/package/AUTHORS +++ b/package/AUTHORS @@ -203,7 +203,6 @@ Chronological list of authors - Patricio Barletta - Mikhail Glagolev 2023 - - Shaivi Malik - Christian Pfaendner - Pratham Chauhan - Meet Brijwani @@ -216,6 +215,7 @@ Chronological list of authors - Josh Vermaas - Xiaoxu Ruan - Egor Marin + - Shaivi Malik External code ------------- From dcf731d8c6ed5f6b4dc833e050d0ac9107b6987b Mon Sep 17 00:00:00 2001 From: g2707 <128649094+g2707@users.noreply.github.com> Date: Thu, 30 Mar 2023 11:38:37 +0530 Subject: [PATCH 7/9] Update transformations.py --- package/MDAnalysis/lib/transformations.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/package/MDAnalysis/lib/transformations.py b/package/MDAnalysis/lib/transformations.py index b6e88e18362..f66d1c6ea2b 100644 --- a/package/MDAnalysis/lib/transformations.py +++ b/package/MDAnalysis/lib/transformations.py @@ -959,7 +959,7 @@ def superimposition_matrix(v0, v1, scaling=False, usesvd=True): # scale: ratio of rms deviations from centroid if scaling: M[:3, :3] *= math.sqrt(np.einsum('ij,ij->',v1,v1) / - np.einsum('ij,ij->',v0,v0)) + np.einsum('ij,ij->',v0,v0)) # translation @@ -1827,4 +1827,4 @@ def rotaxis(a, b): _import_module('_transformations') # Documentation in HTML format can be generated with Epydoc -__docformat__ = "restructuredtext en" \ No newline at end of file +__docformat__ = "restructuredtext en" From 0513994481333e9121c41939246a884ac7be02d1 Mon Sep 17 00:00:00 2001 From: g2707 <128649094+g2707@users.noreply.github.com> Date: Thu, 30 Mar 2023 11:45:45 +0530 Subject: [PATCH 8/9] Update topologyattrs.py --- package/MDAnalysis/core/topologyattrs.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/package/MDAnalysis/core/topologyattrs.py b/package/MDAnalysis/core/topologyattrs.py index 6796c56dd43..197e4fe96dc 100644 --- a/package/MDAnalysis/core/topologyattrs.py +++ b/package/MDAnalysis/core/topologyattrs.py @@ -1719,7 +1719,7 @@ def radius_of_gyration(group, wrap=False, **kwargs): recenteredpos = atomgroup.positions - com rog_sq = np.einsum('i,i->',masses,np.einsum('ij,ij->i', - recenteredpos,recenteredpos))/atomgroup.total_mass() + recenteredpos,recenteredpos))/atomgroup.total_mass() return np.sqrt(rog_sq) @@ -2233,7 +2233,7 @@ def dipole_vector(group, else: recenteredpos = (atomgroup.positions - ref) dipole_vector = np.einsum('ij,ij->j',recenteredpos, - charges[:, np.newaxis]) + charges[:, np.newaxis]) else: (atom_masks, compound_masks, n_compounds) = atomgroup._split_by_compound_indices(compound) @@ -2249,9 +2249,9 @@ def dipole_vector(group, dipole_vector = np.empty((n_compounds, 3), dtype=np.float64) for compound_mask, atom_mask in zip(compound_masks, atom_masks): dipole_vector[compound_mask] = np.einsum('ijk,ijk->ik', - (coords[atom_mask]- - ref[compound_mask][:, None, :]), - chgs[atom_mask][:, :, None]) + (coords[atom_mask]- + ref[compound_mask][:, None, :]), + chgs[atom_mask][:, :, None]) return dipole_vector @@ -3305,4 +3305,4 @@ class CMaps(_Connection): attrname = 'cmaps' singular = 'cmaps' transplants = defaultdict(list) - _n_atoms = 5 \ No newline at end of file + _n_atoms = 5 From 174a3cc4ebd07a444d7d4b10375d4c7840480ed0 Mon Sep 17 00:00:00 2001 From: g2707 <128649094+g2707@users.noreply.github.com> Date: Thu, 30 Mar 2023 11:50:49 +0530 Subject: [PATCH 9/9] Update bat.py --- package/MDAnalysis/analysis/bat.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/package/MDAnalysis/analysis/bat.py b/package/MDAnalysis/analysis/bat.py index 6d7e3cd9e57..80bc7ae4f44 100644 --- a/package/MDAnalysis/analysis/bat.py +++ b/package/MDAnalysis/analysis/bat.py @@ -378,7 +378,7 @@ def _single_frame(self): # Distance between second two root atoms # Angle between root atoms a012 = np.arccos(max(-1.,min(1.,np.einsum('i,i->',v01,v21)/\ - np.sqrt(np.einsum('i,i->',v01,v01)*np.einsum('i,i->',v21,v21))))) + np.sqrt(np.einsum('i,i->',v01,v01)*np.einsum('i,i->',v21,v21))))) # External coordinates e = v01 / r01 phi = np.arctan2(e[1], e[0]) # Polar angle @@ -563,4 +563,4 @@ def Cartesian(self, bat_frame): @property def atoms(self): """The atomgroup for which BAT are computed (read-only property)""" - return self._ag \ No newline at end of file + return self._ag