From dc84819d0a0667a5d4b7753cfca1b9b43901484b Mon Sep 17 00:00:00 2001 From: Richard Gowers Date: Fri, 30 Oct 2015 14:50:13 -0700 Subject: [PATCH 01/17] Removed unused/orphaned code from core/Selections Start of Issue #362 Removed: CASelection (wasn't accessible anyway?) BondedSelection (wasn't fully implemented) NUCLEICXSTAL (only the name existed?) --- package/MDAnalysis/core/Selection.py | 99 +++++++++++++--------------- 1 file changed, 44 insertions(+), 55 deletions(-) diff --git a/package/MDAnalysis/core/Selection.py b/package/MDAnalysis/core/Selection.py index f021698d803..6a5ced09c74 100644 --- a/package/MDAnalysis/core/Selection.py +++ b/package/MDAnalysis/core/Selection.py @@ -695,7 +695,7 @@ class NucleicSelection(Selection): * recognized (CHARMM in Gromacs): 'DA', 'DU', 'DC', 'DG', 'DT' .. versionchanged:: 0.8 - additional Gromacs selections (see also :class:`NucleicXstalSelection`) + additional Gromacs selections """ nucl_res = dict([(x, None) for x in [ 'ADE', 'URA', 'CYT', 'GUA', 'THY', 'DA', 'DC', 'DG', 'DT', 'RA', @@ -769,38 +769,6 @@ def __repr__(self): return "<'NucleicSugarSelection' >" -class CASelection(BackboneSelection): - """Select atoms named CA in protein residues (supposed to be the C-alphas) - """ - - def _apply(self, group): - return set([a for a in group.atoms if (a.name == "CA" and a.resname in self.prot_res)]) - - def __repr__(self): - return "<'CASelection' >" - - -class BondedSelection(Selection): - def __init__(self, sel): - Selection.__init__(self) - self.sel = sel - - def _apply(self, group): - res = self.sel._apply(group) - # Find all the atoms bonded to each - sel = [] - for atom in res: - for b1, b2 in group._bonds: - if atom.index == b1: - sel.append(group.atoms[b2]) - elif atom.index == b2: - sel.append(group.atoms[b1]) - return set(sel) - - def __repr__(self): - return "<'BondedSelection' to " + repr(self.sel) + " >" - - class PropertySelection(Selection): """Some of the possible properties: x, y, z, radius, mass, @@ -910,7 +878,6 @@ class SelectionParser: CYZONE = 'cyzone' POINT = 'point' BYRES = 'byres' - BONDED = 'bonded' BYNUM = 'bynum' PROP = 'prop' ATOM = 'atom' @@ -925,7 +892,6 @@ class SelectionParser: TYPE = 'type' PROTEIN = 'protein' NUCLEIC = 'nucleic' - NUCLEICXSTAL = 'nucleicxstal' BB = 'backbone' NBB = 'nucleicbackbone' BASE = 'nucleicbase' @@ -943,23 +909,51 @@ class SelectionParser: FULLSELGROUP = 'fullgroup' classdict = dict([ - (ALL, AllSelection), (NOT, NotSelection), (AND, AndSelection), (OR, OrSelection), - (SEGID, SegmentNameSelection), (RESID, ResidueIDSelection), (RESNUM, ResnumSelection), - (RESNAME, ResidueNameSelection), (NAME, AtomNameSelection), (ALTLOC, AltlocSelection), - (TYPE, AtomTypeSelection), (BYRES, ByResSelection), - (BYNUM, ByNumSelection), (PROP, PropertySelection), - (AROUND, AroundSelection), (SPHLAYER, SphericalLayerSelection), (SPHZONE, SphericalZoneSelection), - (CYLAYER, CylindricalLayerSelection), (CYZONE, CylindricalZoneSelection), - (POINT, PointSelection), (NUCLEIC, NucleicSelection), (PROTEIN, ProteinSelection), - (BB, BackboneSelection), (NBB, NucleicBackboneSelection), - (BASE, BaseSelection), (SUGAR, NucleicSugarSelection), - #(BONDED, BondedSelection), not supported yet, need a better way to walk the bond lists - (ATOM, AtomSelection), (SELGROUP, SelgroupSelection), (FULLSELGROUP, FullSelgroupSelection), - (SAME, SameSelection), (GLOBAL, GlobalSelection)]) + (ALL, AllSelection), + (NOT, NotSelection), + (AND, AndSelection), + (OR, OrSelection), + (SEGID, SegmentNameSelection), + (RESID, ResidueIDSelection), + (RESNUM, ResnumSelection), + (RESNAME, ResidueNameSelection), + (NAME, AtomNameSelection), + (ALTLOC, AltlocSelection), + (TYPE, AtomTypeSelection), + (BYRES, ByResSelection), + (BYNUM, ByNumSelection), + (PROP, PropertySelection), + (AROUND, AroundSelection), + (SPHLAYER, SphericalLayerSelection), + (SPHZONE, SphericalZoneSelection), + (CYLAYER, CylindricalLayerSelection), + (CYZONE, CylindricalZoneSelection), + (POINT, PointSelection), + (NUCLEIC, NucleicSelection), + (PROTEIN, ProteinSelection), + (BB, BackboneSelection), + (NBB, NucleicBackboneSelection), + (BASE, BaseSelection), + (SUGAR, NucleicSugarSelection), + (ATOM, AtomSelection), + (SELGROUP, SelgroupSelection), + (FULLSELGROUP, FullSelgroupSelection), + (SAME, SameSelection), + (GLOBAL, GlobalSelection)]) associativity = dict([(AND, "left"), (OR, "left")]) precedence = dict( - [(AROUND, 1), (SPHLAYER, 1), (SPHZONE, 1), (CYLAYER, 1), (CYZONE, 1), (POINT, 1), (BYRES, 1), (BONDED, 1), - (SAME, 1), (AND, 3), (OR, 3), (NOT, 5), (GLOBAL, 5)]) + [(AROUND, 1), + (SPHLAYER, 1), + (SPHZONE, 1), + (CYLAYER, 1), + (CYZONE, 1), + (POINT, 1), + (BYRES, 1), + (SAME, 1), + (AND, 3), + (OR, 3), + (NOT, 5), + (GLOBAL, 5)]) # Borg pattern: http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/66531 _shared_state = {} @@ -1047,9 +1041,6 @@ def __parse_subexp(self): y = self.__consume_token() z = self.__consume_token() return self.classdict[op](float(dist), float(x), float(y), float(z)) - elif op == self.BONDED: - exp = self.__parse_expression(self.precedence[op]) - return self.classdict[op](exp) elif op == self.LPAREN: exp = self.__parse_expression(0) self.__expect(self.RPAREN) @@ -1069,8 +1060,6 @@ def __parse_subexp(self): return self.classdict[op]() elif op == self.NUCLEIC: return self.classdict[op]() - elif op == self.NUCLEICXSTAL: - return self.classdict[op]() elif op == self.ALL: return self.classdict[op]() elif op == self.BB: From 7ff8710ddadeb635bc4bdbb902407690c52dc9c2 Mon Sep 17 00:00:00 2001 From: Richard Gowers Date: Fri, 30 Oct 2015 15:14:46 -0700 Subject: [PATCH 02/17] Added tests for point selection --- .../MDAnalysisTests/test_atomselections.py | 37 +++++++++++++++++++ 1 file changed, 37 insertions(+) diff --git a/testsuite/MDAnalysisTests/test_atomselections.py b/testsuite/MDAnalysisTests/test_atomselections.py index 06ecd939e6a..47cfbfeac1c 100644 --- a/testsuite/MDAnalysisTests/test_atomselections.py +++ b/testsuite/MDAnalysisTests/test_atomselections.py @@ -14,6 +14,7 @@ # J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787 # import MDAnalysis +import MDAnalysis as mda import MDAnalysis.core.Selection from MDAnalysis.tests.datafiles import PSF, DCD, PRMpbc, TRJpbc_bz2, PSF_NAMD, PDB_NAMD, GRO, NUCL, TPR, XTC @@ -21,6 +22,8 @@ from numpy.testing import * from nose.plugins.attrib import attr +from MDAnalysis.lib.distances import distance_array + from MDAnalysisTests.plugins.knownfailure import knownfailure class TestSelectionsCHARMM(TestCase): @@ -388,3 +391,37 @@ def test_nucleicsugar(self): rna = self.universe.select_atoms("nucleicsugar") assert_equal(rna.n_residues, 23) assert_equal(rna.n_atoms, rna.n_residues * 5) + + +class TestPointSelection(object): + def setUp(self): + self.u = mda.Universe(GRO) + + def tearDown(self): + del self.u + + def test_point(self): + # The example selection + ag = self.u.select_atoms('point 5.0 5.0 5.0 3.5') + + d = distance_array(np.array([[5.0, 5.0, 5.0]], dtype=np.float32), + self.u.atoms.positions, + box=self.u.dimensions) + + idx = np.where(d < 3.5)[1] + + assert_equal(set(ag.indices), set(idx)) + + def test_point_2(self): + ag1 = self.u.atoms[:10000] + + ag2 = ag1.select_atoms('point 5.0 5.0 5.0 3.5') + + d = distance_array(np.array([[5.0, 5.0, 5.0]], dtype=np.float32), + ag1.positions, + box=self.u.dimensions) + + idx = np.where(d < 3.5)[1] + + assert_equal(set(ag2.indices), set(idx)) + From 6d389eb6e6f0edc400a9ce6cafe5491465bb9c92 Mon Sep 17 00:00:00 2001 From: Richard Gowers Date: Fri, 30 Oct 2015 15:18:14 -0700 Subject: [PATCH 03/17] Removed unused CompositeSelection --- package/MDAnalysis/core/Selection.py | 28 ---------------------------- 1 file changed, 28 deletions(-) diff --git a/package/MDAnalysis/core/Selection.py b/package/MDAnalysis/core/Selection.py index 6a5ced09c74..c895e243ce7 100644 --- a/package/MDAnalysis/core/Selection.py +++ b/package/MDAnalysis/core/Selection.py @@ -468,34 +468,6 @@ def __repr__(self): return "<'PointSelection' " + repr(self.cutoff) + " Ang around " + repr(self.ref) + ">" -class CompositeSelection(Selection): - def __init__(self, name=None, type=None, resname=None, resid=None, segid=None): - Selection.__init__(self) - self.name = name - self.type = type - self.resname = resname - self.resid = resid - self.segid = segid - - def _apply(self, group): - res = [] - for a in group.atoms: - add = True - if self.name is not None and a.name != self.name: - add = False - if self.type is not None and a.type != self.type: - add = False - if self.resname is not None and a.resname != self.resname: - add = False - if self.resid is not None and a.resid != self.resid: - add = False - if self.segid is not None and a.segid != self.segid: - add = False - if add: - res.append(a) - return set(res) - - class AtomSelection(Selection): def __init__(self, name, resid, segid): Selection.__init__(self) From e815630314399e5d9acc82dff74cb721ae754e48 Mon Sep 17 00:00:00 2001 From: Richard Gowers Date: Fri, 30 Oct 2015 15:52:10 -0700 Subject: [PATCH 04/17] Added tests for prop selection --- .../MDAnalysisTests/test_atomselections.py | 64 +++++++++++++++++++ 1 file changed, 64 insertions(+) diff --git a/testsuite/MDAnalysisTests/test_atomselections.py b/testsuite/MDAnalysisTests/test_atomselections.py index 47cfbfeac1c..83f0ca7f8de 100644 --- a/testsuite/MDAnalysisTests/test_atomselections.py +++ b/testsuite/MDAnalysisTests/test_atomselections.py @@ -425,3 +425,67 @@ def test_point_2(self): assert_equal(set(ag2.indices), set(idx)) + +class TestPropSelection(object): + plurals = {'mass': 'masses', + 'charge': 'charges'} + + def _check_lt(self, prop, ag): + setattr(ag[::2], self.plurals[prop], 500.0) + + sel = ag.select_atoms('prop {} < 500.0'.format(prop)) + + assert_equal(set(sel.indices), + set(ag[getattr(ag, self.plurals[prop]) < 500.0].indices)) + + def _check_le(self, prop, ag): + setattr(ag[::2], self.plurals[prop], 500.0) + + sel = ag.select_atoms('prop {} <= 500.0'.format(prop)) + + assert_equal(set(sel.indices), + set(ag[getattr(ag, self.plurals[prop]) <= 500.0].indices)) + + def _check_gt(self, prop, ag): + setattr(ag[::2], self.plurals[prop], 500.0) + + sel = ag.select_atoms('prop {} > 500.0'.format(prop)) + + assert_equal(set(sel.indices), + set(ag[getattr(ag, self.plurals[prop]) > 500.0].indices)) + + def _check_ge(self, prop, ag): + setattr(ag[::2], self.plurals[prop], 500.0) + + sel = ag.select_atoms('prop {} >= 500.0'.format(prop)) + + assert_equal(set(sel.indices), + set(ag[getattr(ag, self.plurals[prop]) >= 500.0].indices)) + + def _check_eq(self, prop, ag): + setattr(ag[::2], self.plurals[prop], 500.0) + + sel = ag.select_atoms('prop {} == 500.0'.format(prop)) + + assert_equal(set(sel.indices), + set(ag[getattr(ag, self.plurals[prop]) == 500.0].indices)) + + def _check_ne(self, prop, ag): + setattr(ag[::2], self.plurals[prop], 500.0) + + sel = ag.select_atoms('prop {} != 500.0'.format(prop)) + + assert_equal(set(sel.indices), + set(ag[getattr(ag, self.plurals[prop]) != 500.0].indices)) + + def test_props(self): + u = mda.Universe(GRO) + + for prop in ['mass', 'charge']: + for ag in [u.atoms, u.atoms[:100]]: + yield self._check_lt, prop, ag + yield self._check_le, prop, ag + yield self._check_gt, prop, ag + yield self._check_ge, prop, ag + yield self._check_eq, prop, ag + yield self._check_ne, prop, ag From 0ad9c10f1cac75fa102f805974dd491551d1f8c8 Mon Sep 17 00:00:00 2001 From: Richard Gowers Date: Sun, 1 Nov 2015 13:24:04 -0800 Subject: [PATCH 05/17] Added BONDED back into to core.Selections Not working yet though Removed dunderscores throughout core.Selections --- package/MDAnalysis/core/Selection.py | 151 ++++++++++++++++----------- 1 file changed, 89 insertions(+), 62 deletions(-) diff --git a/package/MDAnalysis/core/Selection.py b/package/MDAnalysis/core/Selection.py index c895e243ce7..18fcecfbf72 100644 --- a/package/MDAnalysis/core/Selection.py +++ b/package/MDAnalysis/core/Selection.py @@ -485,6 +485,24 @@ def __repr__(self): return "<'AtomSelection' " + repr(self.segid) + " " + repr(self.resid) + " " + repr(self.name) + " >" +class BondedSelection(Selection): + def __init__(self, sel): + Selection.__init__(self) + self.sel = sel + + def _apply(self, group): + res = self.sel._apply(group) + sel = [] + for atom in res: + print group.bonds + for b1, b2 in group.bonds: + if atom == b1: + sel.append(b1) + elif atom == b2: + sel.append(b2) + return set(sel) + + class SelgroupSelection(Selection): def __init__(self, selgroup): Selection.__init__(self) @@ -808,7 +826,7 @@ def _apply(self, group): sel_indices = np.array([a.index for a in Selection._group_atoms]) result_set = group.atoms[np.where(np.in1d(p[sel_indices], p[res_indices]))[0]]._atoms else: - self.__error(self.prop, expected=False) + self._error(self.prop, expected=False) return set(result_set) def __repr__(self): return "<'SameSelection' of "+ repr(self.prop)+" >" @@ -849,6 +867,7 @@ class SelectionParser: CYLAYER = 'cylayer' CYZONE = 'cyzone' POINT = 'point' + BONDED = 'bonded' BYRES = 'byres' BYNUM = 'bynum' PROP = 'prop' @@ -884,6 +903,7 @@ class SelectionParser: (ALL, AllSelection), (NOT, NotSelection), (AND, AndSelection), + (BONDED, BondedSelection), (OR, OrSelection), (SEGID, SegmentNameSelection), (RESID, ResidueIDSelection), @@ -921,6 +941,7 @@ class SelectionParser: (CYZONE, 1), (POINT, 1), (BYRES, 1), + (BONDED, 1), (SAME, 1), (AND, 3), (OR, 3), @@ -935,98 +956,102 @@ def __new__(cls, *p, **k): self.__dict__ = cls._shared_state return self - def __peek_token(self): + def _peek_token(self): """Looks at the next token in our token stream.""" return self.tokens[0] - def __consume_token(self): + def _consume_token(self): """Pops off the next token in our token stream.""" return self.tokens.pop(0) - def __error(self, token, expected=True): + def _error(self, token, expected=True): """Stops parsing and reports an error.""" if expected: raise ParseError("Parsing error- '" + self.selectstr + "'\n" + repr(token) + " expected") else: raise ParseError("Parsing error- '" + self.selectstr + "'\n" + repr(token) + " unexpected") - def __expect(self, token): - if self.__peek_token() == token: - self.__consume_token() + def _expect(self, token): + if self._peek_token() == token: + self._consume_token() else: - self.__error(token) + self._error(token) def parse(self, selectstr, selgroups): self.selectstr = selectstr self.selgroups = selgroups self.tokens = selectstr.replace('(', ' ( ').replace(')', ' ) ').split() + [self.EOF] - parsetree = self.__parse_expression(0) - self.__expect(self.EOF) + parsetree = self._parse_expression(0) + self._expect(self.EOF) return parsetree - def __parse_expression(self, p): - exp1 = self.__parse_subexp() - while self.__peek_token() in (self.AND, self.OR) and self.precedence[self.__peek_token()] >= p: # bin ops - op = self.__consume_token() + def _parse_expression(self, p): + exp1 = self._parse_subexp() + while self._peek_token() in (self.AND, self.OR) and self.precedence[self._peek_token()] >= p: # bin ops + op = self._consume_token() if self.associativity[op] == "right": q = self.precedence[op] else: q = 1 + self.precedence[op] - exp2 = self.__parse_expression(q) + exp2 = self._parse_expression(q) exp1 = self.classdict[op](exp1, exp2) return exp1 - def __parse_subexp(self): - op = self.__consume_token() + def _parse_subexp(self): + op = self._consume_token() if op in (self.NOT, self.BYRES, self.GLOBAL): # unary operators - exp = self.__parse_expression(self.precedence[op]) + exp = self._parse_expression(self.precedence[op]) return self.classdict[op](exp) elif op in (self.AROUND): - dist = self.__consume_token() - exp = self.__parse_expression(self.precedence[op]) + dist = self._consume_token() + exp = self._parse_expression(self.precedence[op]) return self.classdict[op](exp, float(dist)) elif op in (self.SPHLAYER): - inRadius = self.__consume_token() - exRadius = self.__consume_token() - exp = self.__parse_expression(self.precedence[op]) + inRadius = self._consume_token() + exRadius = self._consume_token() + exp = self._parse_expression(self.precedence[op]) return self.classdict[op](exp, float(inRadius), float(exRadius)) elif op in (self.SPHZONE): - dist = self.__consume_token() - exp = self.__parse_expression(self.precedence[op]) + dist = self._consume_token() + exp = self._parse_expression(self.precedence[op]) return self.classdict[op](exp, float(dist)) elif op in (self.CYLAYER): - inRadius = self.__consume_token() - exRadius = self.__consume_token() - zmax = self.__consume_token() - zmin = self.__consume_token() - exp = self.__parse_expression(self.precedence[op]) + inRadius = self._consume_token() + exRadius = self._consume_token() + zmax = self._consume_token() + zmin = self._consume_token() + exp = self._parse_expression(self.precedence[op]) return self.classdict[op](exp, float(inRadius), float(exRadius), float(zmax), float(zmin)) elif op in (self.CYZONE): - exRadius = self.__consume_token() - zmax = self.__consume_token() - zmin = self.__consume_token() - exp = self.__parse_expression(self.precedence[op]) + exRadius = self._consume_token() + zmax = self._consume_token() + zmin = self._consume_token() + exp = self._parse_expression(self.precedence[op]) return self.classdict[op](exp, float(exRadius), float(zmax), float(zmin)) elif op in (self.POINT): - dist = self.__consume_token() - x = self.__consume_token() - y = self.__consume_token() - z = self.__consume_token() + dist = self._consume_token() + x = self._consume_token() + y = self._consume_token() + z = self._consume_token() return self.classdict[op](float(dist), float(x), float(y), float(z)) + elif op == self.BONDED: + exp = self._parse_expression(self.precedence[op]) + return self.classdict[op](exp) elif op == self.LPAREN: - exp = self.__parse_expression(0) - self.__expect(self.RPAREN) + exp = self._parse_expression(0) + self._expect(self.RPAREN) return exp elif op in (self.SEGID, self.RESNAME, self.NAME, self.TYPE, self.ALTLOC): - data = self.__consume_token() + data = self._consume_token() if data in ( - self.LPAREN, - self.RPAREN, self.AND, self.OR, self.NOT, self.SEGID, self.RESID, self.RESNAME, self.NAME, - self.TYPE): - self.__error("Identifier") + self.LPAREN, self.RPAREN, + self.AND, self.OR, self.NOT, + self.SEGID, self.RESID, self.RESNAME, + self.NAME, self.TYPE): + self._error("Identifier") return self.classdict[op](data) elif op in (self.SELGROUP, self.FULLSELGROUP): - grpname = self.__consume_token() + grpname = self._consume_token() return self.classdict[op](self.selgroups[grpname]) elif op == self.PROTEIN: return self.classdict[op]() @@ -1043,7 +1068,7 @@ def __parse_subexp(self): elif op == self.SUGAR: return self.classdict[op]() elif op in (self.RESID, self.RESNUM, self.BYNUM): # can operate on ranges X:Y or X-Y - data = self.__consume_token() + data = self._consume_token() try: lower = int(data) upper = None @@ -1051,18 +1076,18 @@ def __parse_subexp(self): selrange = re.match("(\d+)[:-](\d+)", data) # check if in appropriate format 'lower:upper' or 'lower-upper' if not selrange: - self.__error(op) + self._error(op) lower, upper = map(int, selrange.groups()) return self.classdict[op](lower, upper) elif op == self.PROP: - prop = self.__consume_token() + prop = self._consume_token() if prop == "abs": abs = True - prop = self.__consume_token() + prop = self._consume_token() else: abs = False - oper = self.__consume_token() - value = float(self.__consume_token()) + oper = self._consume_token() + value = float(self._consume_token()) ops = dict([ (self.GT, np.greater), (self.LT, np.less), (self.GE, np.greater_equal), (self.LE, np.less_equal), @@ -1070,21 +1095,23 @@ def __parse_subexp(self): if oper in ops.keys(): return self.classdict[op](prop, ops[oper], value, abs) elif op == self.ATOM: - segid = self.__consume_token() - resid = int(self.__consume_token()) - name = self.__consume_token() + segid = self._consume_token() + resid = int(self._consume_token()) + name = self._consume_token() return self.classdict[op](name, resid, segid) elif op == self.SAME: - prop = self.__consume_token() - self.__expect("as") - if prop in ("name", "type", "resname", "resid", "segid", "mass", "charge", "radius", "bfactor", - "resnum", "residue", "segment", "fragment", "x", "y", "z"): - exp = self.__parse_expression(self.precedence[op]) + prop = self._consume_token() + self._expect("as") + if prop in ("name", "type", "resname", "resid", "segid", + "mass", "charge", "radius", "bfactor", + "resnum", "residue", "segment", "fragment", + "x", "y", "z"): + exp = self._parse_expression(self.precedence[op]) return self.classdict[op](exp, prop) else: - self.__error(prop, expected=False) + self._error(prop, expected=False) else: - self.__error(op, expected=False) + self._error(op, expected=False) # The module level instance Parser = SelectionParser() From 516c8c46285f1022a66c037fa259830d96bac3db Mon Sep 17 00:00:00 2001 From: Richard Gowers Date: Sun, 1 Nov 2015 13:29:47 -0800 Subject: [PATCH 06/17] Removed unused Selection.__init__ --- package/MDAnalysis/core/Selection.py | 30 +--------------------------- 1 file changed, 1 insertion(+), 29 deletions(-) diff --git a/package/MDAnalysis/core/Selection.py b/package/MDAnalysis/core/Selection.py index 18fcecfbf72..80cb1cddaae 100644 --- a/package/MDAnalysis/core/Selection.py +++ b/package/MDAnalysis/core/Selection.py @@ -36,12 +36,7 @@ from ..lib.mdamath import triclinic_vectors -class Selection: - def __init__(self): - # This allows you to build a Selection without tying it to a particular group yet - # Updatable means every timestep - self.update = False # not used at the moment - +class Selection(object): def __repr__(self): return "<" + self.__class__.__name__ + ">" @@ -89,16 +84,12 @@ def apply(self, group): class AllSelection(Selection): - def __init__(self): - Selection.__init__(self) - def _apply(self, group): return set(group.atoms[:]) class NotSelection(Selection): def __init__(self, sel): - Selection.__init__(self) self.sel = sel def _apply(self, group): @@ -111,7 +102,6 @@ def __repr__(self): class AndSelection(Selection): def __init__(self, lsel, rsel): - Selection.__init__(self) self.rsel = rsel self.lsel = lsel @@ -124,7 +114,6 @@ def __repr__(self): class OrSelection(Selection): def __init__(self, lsel, rsel): - Selection.__init__(self) self.rsel = rsel self.lsel = lsel @@ -136,7 +125,6 @@ def __repr__(self): class GlobalSelection(Selection): def __init__(self, sel): - Selection.__init__(self) self.sel = sel def _apply(self, group): @@ -145,7 +133,6 @@ def _apply(self, group): class AroundSelection(Selection): def __init__(self, sel, cutoff, periodic=None): - Selection.__init__(self) self.sel = sel self.cutoff = cutoff self.sqdist = cutoff * cutoff @@ -212,7 +199,6 @@ def __repr__(self): class SphericalLayerSelection(Selection): def __init__(self, sel, inRadius, exRadius, periodic=None): - Selection.__init__(self) self.sel = sel self.inRadius = inRadius self.exRadius = exRadius @@ -273,7 +259,6 @@ def __repr__(self): class SphericalZoneSelection(Selection): def __init__(self, sel, cutoff, periodic=None): - Selection.__init__(self) self.sel = sel self.cutoff = cutoff self.sqdist = cutoff * cutoff @@ -327,7 +312,6 @@ def __repr__(self): class _CylindricalSelection(Selection): def __init__(self, sel, exRadius, zmax, zmin, periodic=None): - Selection.__init__(self) self.sel = sel self.exRadius = exRadius self.exRadiusSq = exRadius * exRadius @@ -397,7 +381,6 @@ def __repr__(self): class CylindricalZoneSelection(_CylindricalSelection): def __init__(self, sel, exRadius, zmax, zmin, periodic=None): - Selection.__init__(self) _CylindricalSelection.__init__(self, sel, exRadius, zmax, zmin, periodic) def __repr__(self): @@ -406,7 +389,6 @@ def __repr__(self): class CylindricalLayerSelection(_CylindricalSelection): def __init__(self, sel, inRadius, exRadius, zmax, zmin, periodic=None): - Selection.__init__(self) _CylindricalSelection.__init__(self, sel, exRadius, zmax, zmin, periodic) self.inRadius = inRadius self.inRadiusSq = inRadius * inRadius @@ -418,7 +400,6 @@ def __repr__(self): class PointSelection(Selection): def __init__(self, x, y, z, cutoff, periodic=None): - Selection.__init__(self) self.ref = np.array((float(x), float(y), float(z))) self.cutoff = float(cutoff) self.cutoffsq = float(cutoff) * float(cutoff) @@ -470,7 +451,6 @@ def __repr__(self): class AtomSelection(Selection): def __init__(self, name, resid, segid): - Selection.__init__(self) self.name = name self.resid = resid self.segid = segid @@ -487,7 +467,6 @@ def __repr__(self): class BondedSelection(Selection): def __init__(self, sel): - Selection.__init__(self) self.sel = sel def _apply(self, group): @@ -505,7 +484,6 @@ def _apply(self, group): class SelgroupSelection(Selection): def __init__(self, selgroup): - Selection.__init__(self) self._grp = selgroup def _apply(self, group): @@ -520,7 +498,6 @@ def __repr__(self): @deprecate(old_name='fullgroup', new_name='global group') class FullSelgroupSelection(Selection): def __init__(self, selgroup): - Selection.__init__(self) self._grp = selgroup def _apply(self, group): @@ -532,7 +509,6 @@ def __repr__(self): class StringSelection(Selection): def __init__(self, field): - Selection.__init__(self) self._field = field def _apply(self, group): @@ -579,7 +555,6 @@ def __init__(self, altLoc): class ByResSelection(Selection): def __init__(self, sel): - Selection.__init__(self) self.sel = sel def _apply(self, group): @@ -597,7 +572,6 @@ def __repr__(self): class _RangeSelection(Selection): def __init__(self, lower, upper): - Selection.__init__(self) self.lower = lower self.upper = upper @@ -765,7 +739,6 @@ class PropertySelection(Selection): """ def __init__(self, prop, operator, value, abs=False): - Selection.__init__(self) self.prop = prop self.operator = operator self.value = value @@ -802,7 +775,6 @@ class SameSelection(Selection): # When adding new keywords here don't forget to also add them to the # case statement under the SAME op, where they are first checked. def __init__(self, sel, prop): - Selection.__init__(self) self.sel = sel self.prop = prop def _apply(self, group): From b4ffb7b0179f634315827a9758b3f594b87143a6 Mon Sep 17 00:00:00 2001 From: Richard Gowers Date: Sun, 1 Nov 2015 13:31:11 -0800 Subject: [PATCH 07/17] removed unused associativity feature from Selection --- package/MDAnalysis/core/Selection.py | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/package/MDAnalysis/core/Selection.py b/package/MDAnalysis/core/Selection.py index 80cb1cddaae..c24caa010f1 100644 --- a/package/MDAnalysis/core/Selection.py +++ b/package/MDAnalysis/core/Selection.py @@ -904,7 +904,6 @@ class SelectionParser: (FULLSELGROUP, FullSelgroupSelection), (SAME, SameSelection), (GLOBAL, GlobalSelection)]) - associativity = dict([(AND, "left"), (OR, "left")]) precedence = dict( [(AROUND, 1), (SPHLAYER, 1), @@ -961,10 +960,7 @@ def _parse_expression(self, p): exp1 = self._parse_subexp() while self._peek_token() in (self.AND, self.OR) and self.precedence[self._peek_token()] >= p: # bin ops op = self._consume_token() - if self.associativity[op] == "right": - q = self.precedence[op] - else: - q = 1 + self.precedence[op] + q = 1 + self.precedence[op] exp2 = self._parse_expression(q) exp1 = self.classdict[op](exp1, exp2) return exp1 From 6a5e447702923ca1c0f9fafee70f72f18430e284 Mon Sep 17 00:00:00 2001 From: Richard Gowers Date: Sun, 1 Nov 2015 14:19:20 -0800 Subject: [PATCH 08/17] Fixed bonded selection --- package/MDAnalysis/core/Selection.py | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/package/MDAnalysis/core/Selection.py b/package/MDAnalysis/core/Selection.py index c24caa010f1..1b344ff29d0 100644 --- a/package/MDAnalysis/core/Selection.py +++ b/package/MDAnalysis/core/Selection.py @@ -473,12 +473,9 @@ def _apply(self, group): res = self.sel._apply(group) sel = [] for atom in res: - print group.bonds - for b1, b2 in group.bonds: - if atom == b1: - sel.append(b1) - elif atom == b2: - sel.append(b2) + for b in group.bonds: + if atom in b: + sel.append(b.partner(atom)) return set(sel) From 5320fb7a5ca71f6983e5992d05916026da2603e1 Mon Sep 17 00:00:00 2001 From: Richard Gowers Date: Sun, 1 Nov 2015 14:54:17 -0800 Subject: [PATCH 09/17] Added test for bonded selections --- .../MDAnalysisTests/test_atomselections.py | 38 ++++++++++++++++--- 1 file changed, 33 insertions(+), 5 deletions(-) diff --git a/testsuite/MDAnalysisTests/test_atomselections.py b/testsuite/MDAnalysisTests/test_atomselections.py index 83f0ca7f8de..d8909b86382 100644 --- a/testsuite/MDAnalysisTests/test_atomselections.py +++ b/testsuite/MDAnalysisTests/test_atomselections.py @@ -13,19 +13,32 @@ # MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations. # J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787 # -import MDAnalysis -import MDAnalysis as mda -import MDAnalysis.core.Selection -from MDAnalysis.tests.datafiles import PSF, DCD, PRMpbc, TRJpbc_bz2, PSF_NAMD, PDB_NAMD, GRO, NUCL, TPR, XTC import numpy as np -from numpy.testing import * +from numpy.testing import( + TestCase, + dec, + assert_equal, + assert_array_almost_equal, + assert_, + assert_array_equal, +) from nose.plugins.attrib import attr +import MDAnalysis +import MDAnalysis as mda +import MDAnalysis.core.Selection from MDAnalysis.lib.distances import distance_array +from MDAnalysis.tests.datafiles import ( + PSF, DCD, + PRMpbc, TRJpbc_bz2, + PSF_NAMD, PDB_NAMD, + GRO, NUCL, TPR, XTC +) from MDAnalysisTests.plugins.knownfailure import knownfailure + class TestSelectionsCHARMM(TestCase): def setUp(self): """Set up the standard AdK system in implicit solvent.""" @@ -489,3 +502,18 @@ def test_props(self): yield self._check_ge, prop, ag yield self._check_eq, prop, ag yield self._check_ne, prop, ag + + +class TestBondedSelection(object): + def setUp(self): + self.u = mda.Universe(PSF, DCD) + + def tearDown(self): + del self.u + + def test_bonded_1(self): + ag = self.u.select_atoms('type 2 and bonded name N') + + assert_(len(ag) == 3) + + From 4ceb1071386262743a6ed14f2aad8997fc1984c7 Mon Sep 17 00:00:00 2001 From: Richard Gowers Date: Sun, 1 Nov 2015 15:09:23 -0800 Subject: [PATCH 10/17] Added warning for BondedSelection without bonds --- package/MDAnalysis/core/Selection.py | 5 +++++ testsuite/MDAnalysisTests/test_atomselections.py | 9 +++++++++ 2 files changed, 14 insertions(+) diff --git a/package/MDAnalysis/core/Selection.py b/package/MDAnalysis/core/Selection.py index 1b344ff29d0..aaf1e7f53a6 100644 --- a/package/MDAnalysis/core/Selection.py +++ b/package/MDAnalysis/core/Selection.py @@ -29,6 +29,7 @@ import numpy as np from numpy.lib.utils import deprecate from Bio.KDTree import KDTree +import warnings from .AtomGroup import AtomGroup, Universe from MDAnalysis.core import flags @@ -471,6 +472,10 @@ def __init__(self, sel): def _apply(self, group): res = self.sel._apply(group) + # Check if we have bonds + if not group.bonds: + warnings.warn("Bonded selection has 0 bonds") + sel = [] for atom in res: for b in group.bonds: diff --git a/testsuite/MDAnalysisTests/test_atomselections.py b/testsuite/MDAnalysisTests/test_atomselections.py index d8909b86382..b2859e37ecc 100644 --- a/testsuite/MDAnalysisTests/test_atomselections.py +++ b/testsuite/MDAnalysisTests/test_atomselections.py @@ -22,13 +22,16 @@ assert_array_almost_equal, assert_, assert_array_equal, + assert_warns, ) from nose.plugins.attrib import attr +import warnings import MDAnalysis import MDAnalysis as mda import MDAnalysis.core.Selection from MDAnalysis.lib.distances import distance_array +from MDAnalysis.core.topologyobjects import TopologyGroup from MDAnalysis.tests.datafiles import ( PSF, DCD, @@ -516,4 +519,10 @@ def test_bonded_1(self): assert_(len(ag) == 3) + def test_nobonds_warns(self): + self.u.bonds = TopologyGroup([]) + + assert_warns(UserWarning, self.u.select_atoms, 'type 2 and bonded name N') + + From 05ca00e65f15c4759eb8f995582d7d0345de3fd9 Mon Sep 17 00:00:00 2001 From: Richard Gowers Date: Sun, 1 Nov 2015 15:18:15 -0800 Subject: [PATCH 11/17] Removed unused bitwise operators from Selections --- package/MDAnalysis/core/Selection.py | 12 ------------ 1 file changed, 12 deletions(-) diff --git a/package/MDAnalysis/core/Selection.py b/package/MDAnalysis/core/Selection.py index aaf1e7f53a6..790d3e8591d 100644 --- a/package/MDAnalysis/core/Selection.py +++ b/package/MDAnalysis/core/Selection.py @@ -41,18 +41,6 @@ class Selection(object): def __repr__(self): return "<" + self.__class__.__name__ + ">" - def __and__(self, other): - return AndSelection(self, other) - - def __or__(self, other): - return OrSelection(self, other) - - def __invert__(self): - return NotSelection(self) - - def __hash__(self): - return hash(repr(self)) - def _apply(self, group): # This is an error raise NotImplementedError("No _apply function defined for " + repr(self.__class__.__name__)) From 3de57fcfa61b35cb9961271ad546dc70ebd84206 Mon Sep 17 00:00:00 2001 From: Richard Gowers Date: Sun, 1 Nov 2015 15:55:29 -0800 Subject: [PATCH 12/17] Added more tests for Selections --- .../MDAnalysisTests/test_atomselections.py | 68 +++++++++++++++++-- 1 file changed, 61 insertions(+), 7 deletions(-) diff --git a/testsuite/MDAnalysisTests/test_atomselections.py b/testsuite/MDAnalysisTests/test_atomselections.py index b2859e37ecc..0d354356e47 100644 --- a/testsuite/MDAnalysisTests/test_atomselections.py +++ b/testsuite/MDAnalysisTests/test_atomselections.py @@ -44,7 +44,10 @@ class TestSelectionsCHARMM(TestCase): def setUp(self): - """Set up the standard AdK system in implicit solvent.""" + """Set up the standard AdK system in implicit solvent. + + Geometry here is orthogonal + """ self.universe = MDAnalysis.Universe(PSF, DCD) def tearDown(self): @@ -158,6 +161,17 @@ def test_cyzone(self): sel = self.universe.select_atoms('cyzone 6.0 10 -10 bynum 1281') assert_equal(len(sel), 166) + def test_point(self): + ag = self.universe.select_atoms('point 5.0 5.0 5.0 3.5') + + d = distance_array(np.array([[5.0, 5.0, 5.0]], dtype=np.float32), + self.universe.atoms.positions, + box=self.universe.dimensions) + + idx = np.where(d < 3.5)[1] + + assert_equal(set(ag.indices), set(idx)) + def test_prop(self): sel = self.universe.select_atoms('prop y <= 16') sel2 = self.universe.select_atoms('prop abs z < 8') @@ -180,7 +194,7 @@ def test_bynum(self): assert_equal(subsel[-1].index, 4) # TODO: - # add more test cases for byres, bynum, point + # add more test cases for byres, bynum # and also for selection keywords such as 'nucleic' def test_same_resname(self): @@ -409,20 +423,59 @@ def test_nucleicsugar(self): assert_equal(rna.n_atoms, rna.n_residues * 5) -class TestPointSelection(object): +class TestTriclinicSelections(object): + """Non-KDTree based selections + + This system has triclinic geometry so won't use KDTree based selections + """ def setUp(self): self.u = mda.Universe(GRO) + self.box = self.u.dimensions def tearDown(self): del self.u - def test_point(self): + def test_around(self): + r1 = self.u.select_atoms('resid 1') + + ag = self.u.select_atoms('around 5.0 resid 1') + + d = distance_array(self.u.atoms.positions, r1.positions, box=self.box) + idx = set(np.where(d < 5.0)[0]) + + # Around doesn't include atoms from the reference group + idx.difference_update(set(r1.indices)) + assert_(idx == set(ag.indices)) + + def test_sphlayer(self): + r1 = self.u.select_atoms('resid 1') + cog = r1.center_of_geometry().reshape(1, 3) + + ag = self.u.select_atoms('sphlayer 2.4 6.0 resid 1') + + d = distance_array(self.u.atoms.positions, cog, box=self.box) + idx = set(np.where((d > 2.4) & (d < 6.0))[0]) + + assert_(idx == set(ag.indices)) + + def test_sphzone(self): + r1 = self.u.select_atoms('resid 1') + cog = r1.center_of_geometry().reshape(1, 3) + + ag = self.u.select_atoms('sphzone 5.0 resid 1') + + d = distance_array(self.u.atoms.positions, cog, box=self.box) + idx = set(np.where(d < 5.0)[0]) + + assert_(idx == set(ag.indices)) + + def test_point_1(self): # The example selection ag = self.u.select_atoms('point 5.0 5.0 5.0 3.5') d = distance_array(np.array([[5.0, 5.0, 5.0]], dtype=np.float32), self.u.atoms.positions, - box=self.u.dimensions) + box=self.box) idx = np.where(d < 3.5)[1] @@ -435,7 +488,7 @@ def test_point_2(self): d = distance_array(np.array([[5.0, 5.0, 5.0]], dtype=np.float32), ag1.positions, - box=self.u.dimensions) + box=self.box) idx = np.where(d < 3.5)[1] @@ -522,7 +575,8 @@ def test_bonded_1(self): def test_nobonds_warns(self): self.u.bonds = TopologyGroup([]) - assert_warns(UserWarning, self.u.select_atoms, 'type 2 and bonded name N') + assert_warns(UserWarning, + self.u.select_atoms, 'type 2 and bonded name N') From 5df50ffd0bf06cde7c99d1d8c1382e91d519269c Mon Sep 17 00:00:00 2001 From: Richard Gowers Date: Sun, 1 Nov 2015 16:01:03 -0800 Subject: [PATCH 13/17] Removed unused expected error from Selections And general code tidying --- package/MDAnalysis/core/Selection.py | 46 ++++++++++++++++------------ 1 file changed, 26 insertions(+), 20 deletions(-) diff --git a/package/MDAnalysis/core/Selection.py b/package/MDAnalysis/core/Selection.py index 790d3e8591d..2e05bfcfba5 100644 --- a/package/MDAnalysis/core/Selection.py +++ b/package/MDAnalysis/core/Selection.py @@ -788,7 +788,7 @@ def _apply(self, group): sel_indices = np.array([a.index for a in Selection._group_atoms]) result_set = group.atoms[np.where(np.in1d(p[sel_indices], p[res_indices]))[0]]._atoms else: - self._error(self.prop, expected=False) + self._error(self.prop) return set(result_set) def __repr__(self): return "<'SameSelection' of "+ repr(self.prop)+" >" @@ -925,12 +925,10 @@ def _consume_token(self): """Pops off the next token in our token stream.""" return self.tokens.pop(0) - def _error(self, token, expected=True): + def _error(self, token): """Stops parsing and reports an error.""" - if expected: - raise ParseError("Parsing error- '" + self.selectstr + "'\n" + repr(token) + " expected") - else: - raise ParseError("Parsing error- '" + self.selectstr + "'\n" + repr(token) + " unexpected") + raise ParseError("Parsing error- '{}'\n {} expected" + "".format(self.selectstr, token)) def _expect(self, token): if self._peek_token() == token: @@ -941,14 +939,16 @@ def _expect(self, token): def parse(self, selectstr, selgroups): self.selectstr = selectstr self.selgroups = selgroups - self.tokens = selectstr.replace('(', ' ( ').replace(')', ' ) ').split() + [self.EOF] + self.tokens = selectstr.replace('(', ' ( ').replace(')', ' ) ') + self.tokens = self.tokens.split() + [self.EOF] parsetree = self._parse_expression(0) self._expect(self.EOF) return parsetree def _parse_expression(self, p): exp1 = self._parse_subexp() - while self._peek_token() in (self.AND, self.OR) and self.precedence[self._peek_token()] >= p: # bin ops + while (self._peek_token() in (self.AND, self.OR) and + self.precedence[self._peek_token()] >= p): # bin ops op = self._consume_token() q = 1 + self.precedence[op] exp2 = self._parse_expression(q) @@ -960,33 +960,37 @@ def _parse_subexp(self): if op in (self.NOT, self.BYRES, self.GLOBAL): # unary operators exp = self._parse_expression(self.precedence[op]) return self.classdict[op](exp) - elif op in (self.AROUND): + elif op == self.AROUND: dist = self._consume_token() exp = self._parse_expression(self.precedence[op]) return self.classdict[op](exp, float(dist)) - elif op in (self.SPHLAYER): + elif op == self.SPHLAYER: inRadius = self._consume_token() exRadius = self._consume_token() exp = self._parse_expression(self.precedence[op]) return self.classdict[op](exp, float(inRadius), float(exRadius)) - elif op in (self.SPHZONE): + elif op == self.SPHZONE: dist = self._consume_token() exp = self._parse_expression(self.precedence[op]) return self.classdict[op](exp, float(dist)) - elif op in (self.CYLAYER): + elif op == self.CYLAYER: inRadius = self._consume_token() exRadius = self._consume_token() zmax = self._consume_token() zmin = self._consume_token() exp = self._parse_expression(self.precedence[op]) - return self.classdict[op](exp, float(inRadius), float(exRadius), float(zmax), float(zmin)) - elif op in (self.CYZONE): + return self.classdict[op](exp, + float(inRadius), float(exRadius), + float(zmax), float(zmin)) + elif op == self.CYZONE: exRadius = self._consume_token() zmax = self._consume_token() zmin = self._consume_token() exp = self._parse_expression(self.precedence[op]) - return self.classdict[op](exp, float(exRadius), float(zmax), float(zmin)) - elif op in (self.POINT): + return self.classdict[op](exp, + float(exRadius), + float(zmax), float(zmin)) + elif op == self.POINT: dist = self._consume_token() x = self._consume_token() y = self._consume_token() @@ -1025,14 +1029,16 @@ def _parse_subexp(self): return self.classdict[op]() elif op == self.SUGAR: return self.classdict[op]() - elif op in (self.RESID, self.RESNUM, self.BYNUM): # can operate on ranges X:Y or X-Y + elif op in (self.RESID, self.RESNUM, self.BYNUM): + # can operate on ranges X:Y or X-Y data = self._consume_token() try: lower = int(data) upper = None except ValueError: + # check if in appropriate format 'lower:upper' or 'lower-upper' selrange = re.match("(\d+)[:-](\d+)", - data) # check if in appropriate format 'lower:upper' or 'lower-upper' + data) if not selrange: self._error(op) lower, upper = map(int, selrange.groups()) @@ -1067,9 +1073,9 @@ def _parse_subexp(self): exp = self._parse_expression(self.precedence[op]) return self.classdict[op](exp, prop) else: - self._error(prop, expected=False) + self._error(prop) else: - self._error(op, expected=False) + self._error(op) # The module level instance Parser = SelectionParser() From 9f5b98afab3ace73f0075ceca744cc16e3746578 Mon Sep 17 00:00:00 2001 From: Richard Gowers Date: Sun, 1 Nov 2015 16:15:46 -0800 Subject: [PATCH 14/17] Removed reprs from core/Selection --- package/MDAnalysis/core/Selection.py | 114 ++++++--------------------- 1 file changed, 22 insertions(+), 92 deletions(-) diff --git a/package/MDAnalysis/core/Selection.py b/package/MDAnalysis/core/Selection.py index 2e05bfcfba5..13c8f449fd4 100644 --- a/package/MDAnalysis/core/Selection.py +++ b/package/MDAnalysis/core/Selection.py @@ -38,19 +38,15 @@ class Selection(object): - def __repr__(self): - return "<" + self.__class__.__name__ + ">" - def _apply(self, group): # This is an error - raise NotImplementedError("No _apply function defined for " + repr(self.__class__.__name__)) + raise NotImplementedError("No _apply function defined") def apply(self, group): # Cache the result for future use # atoms is from Universe # returns AtomGroup - if not (isinstance(group, Universe) or isinstance(group, AtomGroup)): - raise Exception("Must pass in an AtomGroup or Universe to the Selection") + # make a set of all the atoms in the group # XXX this should be static to all the class members Selection._group_atoms = set(group.atoms) @@ -85,9 +81,6 @@ def _apply(self, group): notsel = self.sel._apply(group) return (set(group.atoms[:]) - notsel) - def __repr__(self): - return "<'NotSelection' " + repr(self.sel) + ">" - class AndSelection(Selection): def __init__(self, lsel, rsel): @@ -97,9 +90,6 @@ def __init__(self, lsel, rsel): def _apply(self, group): return self.lsel._apply(group) & self.rsel._apply(group) - def __repr__(self): - return "<'AndSelection' " + repr(self.lsel) + "," + repr(self.rsel) + ">" - class OrSelection(Selection): def __init__(self, lsel, rsel): @@ -109,8 +99,6 @@ def __init__(self, lsel, rsel): def _apply(self, group): return self.lsel._apply(group) | self.rsel._apply(group) - def __repr__(self): - return "<'OrSelection' " + repr(self.lsel) + "," + repr(self.rsel) + ">" class GlobalSelection(Selection): def __init__(self, sel): @@ -182,9 +170,6 @@ def _apply_distmat(self, group): np.any(dist <= self.cutoff, axis=1).nonzero()[0]] # make list np array and use fancy indexing? return set(res_atoms) - def __repr__(self): - return "<'AroundSelection' " + repr(self.cutoff) + " around " + repr(self.sel) + ">" - class SphericalLayerSelection(Selection): def __init__(self, sel, inRadius, exRadius, periodic=None): @@ -241,10 +226,6 @@ def _apply_distmat(self, group): box = None return set(res_atoms) - def __repr__(self): - return "<'SphericalLayerSelection' inner radius " + repr(self.inRadius) + ", external radius " + repr( - self.exRadius) + " centered in " + repr(self.sel) + ">" - class SphericalZoneSelection(Selection): def __init__(self, sel, cutoff, periodic=None): @@ -295,9 +276,6 @@ def _apply_distmat(self, group): box = None return set(res_atoms) - def __repr__(self): - return "<'SphericalZoneSelection' radius " + repr(self.cutoff) + " centered in " + repr(self.sel) + ">" - class _CylindricalSelection(Selection): def __init__(self, sel, exRadius, zmax, zmin, periodic=None): @@ -364,17 +342,11 @@ def _apply_distmat(self, group): res_atoms = set(Selection._group_atoms_list[ndx] for ndx in ndxs) return res_atoms - def __repr__(self): - return "<'CylindricalSelection' radius " + repr(self.exRadius) + ", zmax " + repr( - self.zmax) + ", zmin " + repr(self.zmin) + ">" class CylindricalZoneSelection(_CylindricalSelection): def __init__(self, sel, exRadius, zmax, zmin, periodic=None): _CylindricalSelection.__init__(self, sel, exRadius, zmax, zmin, periodic) - def __repr__(self): - return "<'CylindricalZoneSelection' radius " + repr(self.exRadius) + ", zmax " + repr( - self.zmax) + ", zmin " + repr(self.zmin) + ">" class CylindricalLayerSelection(_CylindricalSelection): def __init__(self, sel, inRadius, exRadius, zmax, zmin, periodic=None): @@ -382,10 +354,6 @@ def __init__(self, sel, inRadius, exRadius, zmax, zmin, periodic=None): self.inRadius = inRadius self.inRadiusSq = inRadius * inRadius - def __repr__(self): - return "<'CylindricalLayerSelection' inner radius " + repr(self.inRadius) + ", external radius " + repr( - self.exRadius) + ", zmax " + repr(self.zmax) + ", zmin " + repr(self.zmin) + ">" - class PointSelection(Selection): def __init__(self, x, y, z, cutoff, periodic=None): @@ -434,9 +402,6 @@ def _apply_distmat(self, group): # make list np array and use fancy indexing? return set(res_atoms) - def __repr__(self): - return "<'PointSelection' " + repr(self.cutoff) + " Ang around " + repr(self.ref) + ">" - class AtomSelection(Selection): def __init__(self, name, resid, segid): @@ -450,9 +415,6 @@ def _apply(self, group): return set([a]) return set([]) - def __repr__(self): - return "<'AtomSelection' " + repr(self.segid) + " " + repr(self.resid) + " " + repr(self.name) + " >" - class BondedSelection(Selection): def __init__(self, sel): @@ -481,9 +443,6 @@ def _apply(self, group): res_atoms = [i for i in self._grp if i.index in common] return set(res_atoms) - def __repr__(self): - return "<" + repr(self.__class__.__name__) + ">" - @deprecate(old_name='fullgroup', new_name='global group') class FullSelgroupSelection(Selection): @@ -493,9 +452,6 @@ def __init__(self, selgroup): def _apply(self, group): return set(self._grp._atoms) - def __repr__(self): - return "<" + repr(self.__class__.__name__) + ">" - class StringSelection(Selection): def __init__(self, field): @@ -507,10 +463,8 @@ def _apply(self, group): wc_pos = value.find('*') # This returns -1, so if it's not in value then use the whole of value if wc_pos == -1: wc_pos = None - return set([a for a in group.atoms if getattr(a, self._field)[:wc_pos] == value[:wc_pos]]) - - def __repr__(self): - return "<" + repr(self.__class__.__name__) + ": " + repr(getattr(self, self._field)) + ">" + return set([a for a in group.atoms + if getattr(a, self._field)[:wc_pos] == value[:wc_pos]]) class AtomNameSelection(StringSelection): @@ -556,18 +510,12 @@ def _apply(self, group): sel.append(atom) return set(sel) - def __repr__(self): - return "<'ByResSelection'>" - class _RangeSelection(Selection): def __init__(self, lower, upper): self.lower = lower self.upper = upper - def __repr__(self): - return "<'" + self.__class__.__name__ + "' " + repr(self.lower) + ":" + repr(self.upper) + " >" - class ResidueIDSelection(_RangeSelection): def _apply(self, group): @@ -634,9 +582,6 @@ class ProteinSelection(Selection): def _apply(self, group): return set([a for a in group.atoms if a.resname in self.prot_res]) - def __repr__(self): - return "<'ProteinSelection' >" - class NucleicSelection(Selection): """A nucleic selection consists of all atoms in nucleic acid residues with recognized residue names. @@ -658,26 +603,21 @@ class NucleicSelection(Selection): def _apply(self, group): return set([a for a in group.atoms if a.resname in self.nucl_res]) - def __repr__(self): - return "<'NucleicSelection' >" - class BackboneSelection(ProteinSelection): """A BackboneSelection contains all atoms with name 'N', 'CA', 'C', 'O'. - This excludes OT* on C-termini (which are included by, eg VMD's backbone selection). + This excludes OT* on C-termini + (which are included by, eg VMD's backbone selection). """ bb_atoms = dict([(x, None) for x in ['N', 'CA', 'C', 'O']]) def _apply(self, group): return set([a for a in group.atoms if (a.name in self.bb_atoms and a.resname in self.prot_res)]) - def __repr__(self): - return "<'BackboneSelection' >" - class NucleicBackboneSelection(NucleicSelection): - """A NucleicBackboneSelection contains all atoms with name "P", "C5'", C3'", "O3'", "O5'". + """Contains all atoms with name "P", "C5'", C3'", "O3'", "O5'". These atoms are only recognized if they are in a residue matched by the :class:`NucleicSelection`. @@ -685,10 +625,9 @@ class NucleicBackboneSelection(NucleicSelection): bb_atoms = dict([(x, None) for x in ["P", "C5'", "C3'", "O3'", "O5'"]]) def _apply(self, group): - return set([a for a in group.atoms if (a.name in self.bb_atoms and a.resname in self.nucl_res)]) - - def __repr__(self): - return "<'NucleicBackboneSelection' >" + return set([a for a in group.atoms + if (a.name in self.bb_atoms and + a.resname in self.nucl_res)]) class BaseSelection(NucleicSelection): @@ -705,29 +644,26 @@ class BaseSelection(NucleicSelection): 'O2', 'N4', 'O4', 'C5M']]) def _apply(self, group): - return set([a for a in group.atoms if (a.name in self.base_atoms and a.resname in self.nucl_res)]) - - def __repr__(self): - return "<'BaseSelection' >" + return set([a for a in group.atoms + if (a.name in self.base_atoms and + a.resname in self.nucl_res)]) class NucleicSugarSelection(NucleicSelection): - """A NucleicSugarSelection contains all atoms with name C1', C2', C3', C4', O2', O4', O3'. + """Contains all atoms with name C1', C2', C3', C4', O2', O4', O3'. """ sug_atoms = dict([(x, None) for x in ['C1\'', 'C2\'', 'C3\'', 'C4\'', 'O4\'']]) def _apply(self, group): - return set([a for a in group.atoms if (a.name in self.sug_atoms and a.resname in self.nucl_res)]) - - def __repr__(self): - return "<'NucleicSugarSelection' >" + return set([a for a in group.atoms + if (a.name in self.sug_atoms and + a.resname in self.nucl_res)]) class PropertySelection(Selection): """Some of the possible properties: x, y, z, radius, mass, """ - def __init__(self, prop, operator, value, abs=False): self.prop = prop self.operator = operator @@ -748,18 +684,13 @@ def _apply(self, group): res = res[0] result_set = [group.atoms[i] for i in res] elif self.prop == "mass": - result_set = [a for a in group.atoms if self.operator(a.mass, self.value)] + result_set = [a for a in group.atoms + if self.operator(a.mass, self.value)] elif self.prop == "charge": - result_set = [a for a in group.atoms if self.operator(a.charge, self.value)] + result_set = [a for a in group.atoms + if self.operator(a.charge, self.value)] return set(result_set) - def __repr__(self): - if self.abs: - abs_str = " abs " - else: - abs_str = "" - return "<'PropertySelection' " + abs_str + repr(self.prop) + " " + repr(self.operator.__name__) + " " + repr( - self.value) + ">" class SameSelection(Selection): # When adding new keywords here don't forget to also add them to the @@ -767,6 +698,7 @@ class SameSelection(Selection): def __init__(self, sel, prop): self.sel = sel self.prop = prop + def _apply(self, group): res = self.sel._apply(group) if not res: @@ -790,8 +722,6 @@ def _apply(self, group): else: self._error(self.prop) return set(result_set) - def __repr__(self): - return "<'SameSelection' of "+ repr(self.prop)+" >" class ParseError(Exception): From d16c6bd3f82cfd4b37b1bce951d0ed2aa3f49583 Mon Sep 17 00:00:00 2001 From: Richard Gowers Date: Mon, 2 Nov 2015 08:17:24 -0800 Subject: [PATCH 15/17] Fixed PointSelection KDTree (CNS was removed) Added DistanceSelection base class for distance based selections Appropriate backend is now chosen on init, useful so we can hack in and choose a different one for testing More tests for distance based selections --- package/MDAnalysis/core/Selection.py | 129 +++++++---------- .../MDAnalysisTests/test_atomselections.py | 133 +++++++++++++++++- 2 files changed, 185 insertions(+), 77 deletions(-) diff --git a/package/MDAnalysis/core/Selection.py b/package/MDAnalysis/core/Selection.py index 13c8f449fd4..09940bcc2f8 100644 --- a/package/MDAnalysis/core/Selection.py +++ b/package/MDAnalysis/core/Selection.py @@ -108,20 +108,34 @@ def _apply(self, group): sel = self.sel._apply(group.universe) return sel -class AroundSelection(Selection): - def __init__(self, sel, cutoff, periodic=None): - self.sel = sel - self.cutoff = cutoff - self.sqdist = cutoff * cutoff - if periodic is None: - self.periodic = flags['use_periodic_selections'] +class DistanceSelection(Selection): + """Base class for distance search based selections - def _apply(self, group): - # make choosing _fast/_slow configurable (while testing) + Grabs the flags for this selection + - 'use_KDTree_routines' + - 'use_periodic_selections' + + Populates the `_apply` method with either + - _apply_KDTree + - _apply_distmat + """ + def __init__(self): if flags['use_KDTree_routines'] in (True, 'fast', 'always'): - return self._apply_KDTree(group) + self._apply = self._apply_KDTree else: - return self._apply_distmat(group) + self._apply = self._apply_distmat + + self.periodic = flags['use_periodic_selections'] + # KDTree doesn't support periodic + if self.periodic: + self._apply = self._apply_distmat + + +class AroundSelection(DistanceSelection): + def __init__(self, sel, cutoff): + super(AroundSelection, self).__init__() + self.sel = sel + self.cutoff = cutoff def _apply_KDTree(self, group): """KDTree based selection is about 7x faster than distmat for typical problems. @@ -160,7 +174,7 @@ def _apply_distmat(self, group): sel_coor = Selection.coord[sel_indices] sys_coor = Selection.coord[sys_indices] if self.periodic: - box = group.dimensions[:3] # ignored with KDTree + box = group.dimensions[:3] else: box = None @@ -171,20 +185,12 @@ def _apply_distmat(self, group): return set(res_atoms) -class SphericalLayerSelection(Selection): - def __init__(self, sel, inRadius, exRadius, periodic=None): +class SphericalLayerSelection(DistanceSelection): + def __init__(self, sel, inRadius, exRadius): + super(SphericalLayerSelection, self).__init__() self.sel = sel self.inRadius = inRadius self.exRadius = exRadius - if periodic is None: - self.periodic = flags['use_periodic_selections'] - - def _apply(self, group): - # make choosing _fast/_slow configurable (while testing) - if flags['use_KDTree_routines'] in (True, 'fast', 'always'): - return self._apply_KDTree(group) - else: - return self._apply_distmat(group) def _apply_KDTree(self, group): """Selection using KDTree but periodic = True not supported. @@ -195,8 +201,7 @@ def _apply_KDTree(self, group): sel_atoms = self.sel._apply(group) sel_CoG = AtomGroup(sel_atoms).center_of_geometry() self.ref = np.array((sel_CoG[0], sel_CoG[1], sel_CoG[2])) - if self.periodic: - pass # or warn? -- no periodic functionality with KDTree search + kdtree = KDTree(dim=3, bucket_size=10) kdtree.set_coords(sys_coor) @@ -206,6 +211,7 @@ def _apply_KDTree(self, group): found_IntIndices = kdtree.get_indices() found_indices = list(set(found_ExtIndices) - set(found_IntIndices)) res_atoms = [self._group_atoms_list[i] for i in found_indices] + return set(res_atoms) def _apply_distmat(self, group): @@ -220,27 +226,15 @@ def _apply_distmat(self, group): str(sel_CoG[2]) + " " + str(self.inRadius) sel = sys_ag.select_atoms(sel_CoG_str) res_atoms = AtomGroup(set(sel)) - if self.periodic: - box = group.dimensions[:3] # ignored with KDTree - else: - box = None + return set(res_atoms) -class SphericalZoneSelection(Selection): - def __init__(self, sel, cutoff, periodic=None): +class SphericalZoneSelection(DistanceSelection): + def __init__(self, sel, cutoff): + super(SphericalZoneSelection, self).__init__() self.sel = sel self.cutoff = cutoff - self.sqdist = cutoff * cutoff - if periodic is None: - self.periodic = flags['use_periodic_selections'] - - def _apply(self, group): - # make choosing _fast/_slow configurable (while testing) - if flags['use_KDTree_routines'] in (True, 'fast', 'always'): - return self._apply_KDTree(group) - else: - return self._apply_distmat(group) def _apply_KDTree(self, group): """Selection using KDTree but periodic = True not supported. @@ -251,8 +245,6 @@ def _apply_KDTree(self, group): sel_atoms = self.sel._apply(group) # group is wrong, should be universe (?!) sel_CoG = AtomGroup(sel_atoms).center_of_geometry() self.ref = np.array((sel_CoG[0], sel_CoG[1], sel_CoG[2])) - if self.periodic: - pass # or warn? -- no periodic functionality with KDTree search kdtree = KDTree(dim=3, bucket_size=10) kdtree.set_coords(sys_coor) @@ -270,15 +262,12 @@ def _apply_distmat(self, group): self.cutoff) sel = sys_ag.select_atoms(sel_CoG_str) res_atoms = AtomGroup(set(sel)) - if self.periodic: - box = group.dimensions[:3] # ignored with KDTree - else: - box = None + return set(res_atoms) class _CylindricalSelection(Selection): - def __init__(self, sel, exRadius, zmax, zmin, periodic=None): + def __init__(self, sel, exRadius, zmax, zmin): self.sel = sel self.exRadius = exRadius self.exRadiusSq = exRadius * exRadius @@ -287,10 +276,6 @@ def __init__(self, sel, exRadius, zmax, zmin, periodic=None): self.periodic = flags['use_periodic_selections'] def _apply(self, group): - #KDTree function not implementable - return self._apply_distmat(group) - - def _apply_distmat(self, group): sel_atoms = self.sel._apply(group) sel_CoG = AtomGroup(sel_atoms).center_of_geometry() coords = AtomGroup(Selection._group_atoms_list).positions @@ -344,44 +329,36 @@ def _apply_distmat(self, group): class CylindricalZoneSelection(_CylindricalSelection): - def __init__(self, sel, exRadius, zmax, zmin, periodic=None): - _CylindricalSelection.__init__(self, sel, exRadius, zmax, zmin, periodic) + def __init__(self, sel, exRadius, zmax, zmin): + _CylindricalSelection.__init__(self, sel, exRadius, zmax, zmin) class CylindricalLayerSelection(_CylindricalSelection): - def __init__(self, sel, inRadius, exRadius, zmax, zmin, periodic=None): - _CylindricalSelection.__init__(self, sel, exRadius, zmax, zmin, periodic) + def __init__(self, sel, inRadius, exRadius, zmax, zmin): + _CylindricalSelection.__init__(self, sel, exRadius, zmax, zmin) self.inRadius = inRadius self.inRadiusSq = inRadius * inRadius -class PointSelection(Selection): - def __init__(self, x, y, z, cutoff, periodic=None): - self.ref = np.array((float(x), float(y), float(z))) - self.cutoff = float(cutoff) - self.cutoffsq = float(cutoff) * float(cutoff) - if periodic is None: - self.periodic = flags['use_periodic_selections'] - - def _apply(self, group): - # make choosing _fast/_slow configurable (while testing) - if flags['use_KDTree_routines'] in ('always',): - return self._apply_KDTree(group) - else: - return self._apply_distmat(group) +class PointSelection(DistanceSelection): + def __init__(self, x, y, z, cutoff): + super(PointSelection, self).__init__() + self.ref = np.array([x, y, z]) + self.cutoff = cutoff def _apply_KDTree(self, group): """Selection using KDTree but periodic = True not supported. (KDTree routine is ca 15% slower than the distance matrix one) """ sys_indices = np.array([a.index for a in self._group_atoms_list]) - sys_coor = Selection.coord[sys_indices] - if self.periodic: - pass # or warn? -- no periodic functionality with KDTree search - CNS = CoordinateNeighborSearch(sys_coor) # cache the KDTree for this selection/frame? - found_indices = CNS.search(self.ref, self.cutoff) - res_atoms = [self._group_atoms_list[i] for i in found_indices] # make list np array and use fancy indexing? + kdtree = KDTree(dim=3, bucket_size=10) + kdtree.set_coords(Selection.coord[sys_indices]) + kdtree.search(self.ref, self.cutoff) + found_indices = kdtree.get_indices() + + # make list np array and use fancy indexing? + res_atoms = [self._group_atoms_list[i] for i in found_indices] return set(res_atoms) def _apply_distmat(self, group): diff --git a/testsuite/MDAnalysisTests/test_atomselections.py b/testsuite/MDAnalysisTests/test_atomselections.py index 0d354356e47..17cdcec112e 100644 --- a/testsuite/MDAnalysisTests/test_atomselections.py +++ b/testsuite/MDAnalysisTests/test_atomselections.py @@ -32,12 +32,14 @@ import MDAnalysis.core.Selection from MDAnalysis.lib.distances import distance_array from MDAnalysis.core.topologyobjects import TopologyGroup +from MDAnalysis.core.Selection import Parser from MDAnalysis.tests.datafiles import ( PSF, DCD, PRMpbc, TRJpbc_bz2, PSF_NAMD, PDB_NAMD, - GRO, NUCL, TPR, XTC + GRO, NUCL, TPR, XTC, + TRZ_psf, TRZ, ) from MDAnalysisTests.plugins.knownfailure import knownfailure @@ -423,6 +425,135 @@ def test_nucleicsugar(self): assert_equal(rna.n_atoms, rna.n_residues * 5) +class TestDistanceSelections(object): + """Both KDTree and distmat selections on orthogonal system + + Selections to check: + - Around + - SphericalLayer + - SphericalZone + - Point + + Cylindrical methods don't use KDTree + """ + def setUp(self): + self.u = mda.Universe(TRZ_psf, TRZ) + + def tearDown(self): + del self.u + + def choosemeth(self, sel, meth): + """hack in the desired apply method""" + if meth == 'kdtree': + sel._apply = sel._apply_KDTree + elif meth == 'distmat': + sel._apply = sel._apply_distmat + return sel + + def _check_around(self, meth): + sel = Parser.parse('around 5.0 resid 1', self.u.atoms) + sel = self.choosemeth(sel, meth) + result = sel.apply(self.u.atoms) + + r1 = self.u.select_atoms('resid 1') + cog = r1.center_of_geometry().reshape(1, 3) + + # Ugly hack + # KDTree doesn't do periodicity + # distmat does + if meth == 'distmat': + box = self.u.dimensions + else: + box = None + + d = distance_array(self.u.atoms.positions, r1.positions, box=box) + ref = set(np.where(d < 5.0)[0]) + + # Around doesn't include atoms from the reference group + ref.difference_update(set(r1.indices)) + assert_(ref == set(result.indices)) + + + def test_around(self): + for meth in ['kdtree', 'distmat']: + yield self._check_around, meth + + def _check_spherical_layer(self, meth): + sel = Parser.parse('sphlayer 2.4 6.0 resid 1' , self.u.atoms) + sel = self.choosemeth(sel, meth) + result = sel.apply(self.u.atoms) + + r1 = self.u.select_atoms('resid 1') + cog = r1.center_of_geometry().reshape(1, 3) + + # Ugly hack + # KDTree doesn't do periodicity + # distmat does + if meth == 'distmat': + box = self.u.dimensions + else: + box = None + + d = distance_array(self.u.atoms.positions, cog, box=box) + ref = set(np.where((d > 2.4) & (d < 6.0))[0]) + + assert_(ref == set(result.indices)) + + def test_spherical_layer(self): + for meth in ['kdtree', 'distmat']: + yield self._check_spherical_layer, meth + + def _check_spherical_zone(self, meth): + sel = Parser.parse('sphzone 5.0 resid 1', self.u.atoms) + sel = self.choosemeth(sel, meth) + result = sel.apply(self.u.atoms) + + r1 = self.u.select_atoms('resid 1') + cog = r1.center_of_geometry().reshape(1, 3) + + # Ugly hack + # KDTree doesn't do periodicity + # distmat does + if meth == 'distmat': + box = self.u.dimensions + else: + box = None + + d = distance_array(self.u.atoms.positions, cog, box=box) + ref = set(np.where(d < 5.0)[0]) + + assert_(ref == set(result.indices)) + + def test_spherical_zone(self): + for meth in ['kdtree', 'distmat']: + yield self._check_spherical_zone, meth + + def _check_point(self, meth): + sel = Parser.parse('point 5.0 5.0 5.0 3.0', self.u.atoms) + sel = self.choosemeth(sel, meth) + result = sel.apply(self.u.atoms) + + # Ugly hack + # KDTree doesn't do periodicity + # distmat does + if meth == 'distmat': + box = self.u.dimensions + else: + box = None + + d = distance_array(np.array([[5.0, 5.0, 5.0]], dtype=np.float32), + self.u.atoms.positions, + box=box) + ref = set(np.where(d < 3.5)[1]) + + assert_(ref == set(result.indices)) + + def test_point(self): + # Doesn't work with KDTree currently + for meth in ['kdtree', 'distmat']: + yield self._check_point, meth + + class TestTriclinicSelections(object): """Non-KDTree based selections From c1ba9d20c08475aa3594c5d1dd1dfb82785ad645 Mon Sep 17 00:00:00 2001 From: Richard Gowers Date: Mon, 2 Nov 2015 08:22:54 -0800 Subject: [PATCH 16/17] Added test for byres --- testsuite/MDAnalysisTests/test_atomselections.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/testsuite/MDAnalysisTests/test_atomselections.py b/testsuite/MDAnalysisTests/test_atomselections.py index 17cdcec112e..e6862ebbf91 100644 --- a/testsuite/MDAnalysisTests/test_atomselections.py +++ b/testsuite/MDAnalysisTests/test_atomselections.py @@ -196,8 +196,11 @@ def test_bynum(self): assert_equal(subsel[-1].index, 4) # TODO: - # add more test cases for byres, bynum # and also for selection keywords such as 'nucleic' + def test_byres(self): + sel = self.universe.select_atoms('byres bynum 0:5') + + assert_equal(len(sel), len(self.universe.residues[0])) def test_same_resname(self): """Test the 'same ... as' construct (Issue 217)""" @@ -232,8 +235,6 @@ def test_same_segment(self): #cleanup self.universe.residues.set_segids("4AKE") - - def test_empty_selection(self): """Test that empty selection can be processed (see Issue 12)""" assert_equal(len(self.universe.select_atoms('resname TRP')), 0) # no Trp in AdK From 88f4798c51c9633b167c69d1c7ae7580fe73c529 Mon Sep 17 00:00:00 2001 From: Richard Gowers Date: Mon, 2 Nov 2015 09:59:49 -0800 Subject: [PATCH 17/17] Reworked spherical selections, more efficient Added distance selection tests for with/without PBC Finishes Issue #362 --- package/CHANGELOG | 4 +- package/MDAnalysis/core/AtomGroup.py | 114 ++++++++++-------- package/MDAnalysis/core/Selection.py | 39 +++--- .../MDAnalysisTests/test_atomselections.py | 85 +++++-------- 4 files changed, 126 insertions(+), 116 deletions(-) diff --git a/package/CHANGELOG b/package/CHANGELOG index e875064fe21..691a583d81d 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -30,6 +30,8 @@ Enhancement * GRO file reading approximately 50% faster for large files (Issue #212) * GRO file writing now will write velocities where possible (Issue #494) + * Added bonded selection (Issue #362) + * Spherical layer and spherical zone selections now much faster (Issue #362) Changes * built html doc files are no longer version controlled (Issue #491) @@ -45,7 +47,7 @@ Fixes * Fixed non-compliant Amber NCDFWriter (Issue #488) * Fixed many Timestep methods failing when positions weren't present (Issue #512) - + * Fixed PointSelection using KDTree (Issue #362) 10/08/15 diff --git a/package/MDAnalysis/core/AtomGroup.py b/package/MDAnalysis/core/AtomGroup.py index 74cb5542655..c40365b02ce 100644 --- a/package/MDAnalysis/core/AtomGroup.py +++ b/package/MDAnalysis/core/AtomGroup.py @@ -4383,38 +4383,40 @@ def select_atoms(self, sel, *othersel, **selgroups): **Simple selections** protein, backbone, nucleic, nucleicbackbone - selects all atoms that belong to a standard set of residues; a protein - is identfied by a hard-coded set of residue names so it may not - work for esoteric residues. + selects all atoms that belong to a standard set of residues; + a protein is identfied by a hard-coded set of residue names so + it may not work for esoteric residues. segid *seg-name* - select by segid (as given in the topology), e.g. ``segid 4AKE`` or ``segid DMPC`` + select by segid (as given in the topology), e.g. ``segid 4AKE`` + or ``segid DMPC`` resid *residue-number-range* - resid can take a single residue number or a range of numbers. A range - consists of two numbers separated by a colon (inclusive) such - as ``resid 1:5``. A residue number ("resid") is taken directly from the - topology. + resid can take a single residue number or a range of numbers. A + range consists of two numbers separated by a colon (inclusive) + such as ``resid 1:5``. A residue number ("resid") is taken + directly from the topology. resnum *resnum-number-range* - resnum is the canonical residue number; typically it is set to the residue id - in the original PDB structure. + resnum is the canonical residue number; typically it is set to + the residue id in the original PDB structure. resname *residue-name* select by residue name, e.g. ``resname LYS`` name *atom-name* - select by atom name (as given in the topology). Often, this is force - field dependent. Example: ``name CA`` (for Cα atoms) or ``name OW`` (for SPC water oxygen) + select by atom name (as given in the topology). Often, this is + force field dependent. Example: ``name CA`` (for Cα atoms) + or ``name OW`` (for SPC water oxygen) type *atom-type* - select by atom type; this is either a string or a number and depends on - the force field; it is read from the topology file (e.g. the CHARMM PSF - file contains numeric atom types). It has non-sensical values when a - PDB or GRO file is used as a topology. + select by atom type; this is either a string or a number and + depends on the force field; it is read from the topology file + (e.g. the CHARMM PSF file contains numeric atom types). It has + non-sensical values when a PDB or GRO file is used as a topology atom *seg-name* *residue-number* *atom-name* a selector for a single atom consisting of segid resid atomname, - e.g. ``DMPC 1 C2`` selects the C2 carbon of the first residue of the DMPC - segment + e.g. ``DMPC 1 C2`` selects the C2 carbon of the first residue of + the DMPC segment altloc *alternative-location* a selection for atoms where alternative locations are available, which is often the case with high-resolution crystal structures - e.g. `resid 4 and resname ALA and altloc B` selects only the atoms - of ALA-4 that have an altloc B record. + e.g. `resid 4 and resname ALA and altloc B` selects only the + atoms of ALA-4 that have an altloc B record. **Boolean** @@ -4423,65 +4425,81 @@ def select_atoms(self, sel, *othersel, **selgroups): all atoms that aren't part of a protein and, or - combine two selections according to the rules of boolean algebra, - e.g. ``protein and not (resname ALA or resname LYS)`` selects all atoms - that belong to a protein, but are not in a lysine or alanine residue + combine two selections according to the rules of boolean + algebra, e.g. ``protein and not (resname ALA or resname LYS)`` + selects all atoms that belong to a protein, but are not in a + lysine or alanine residue **Geometric** around *distance* *selection* selects all atoms a certain cutoff away from another selection, - e.g. ``around 3.5 protein`` selects all atoms not belonging to protein - that are within 3.5 Angstroms from the protein + e.g. ``around 3.5 protein`` selects all atoms not belonging to + protein that are within 3.5 Angstroms from the protein point *x* *y* *z* *distance* selects all atoms within a cutoff of a point in space, make sure - coordinate is separated by spaces, e.g. ``point 5.0 5.0 5.0 3.5`` selects - all atoms within 3.5 Angstroms of the coordinate (5.0, 5.0, 5.0) + coordinate is separated by spaces, + e.g. ``point 5.0 5.0 5.0 3.5`` selects all atoms within 3.5 + Angstroms of the coordinate (5.0, 5.0, 5.0) prop [abs] *property* *operator* *value* - selects atoms based on position, using *property* **x**, **y**, or - **z** coordinate. Supports the **abs** keyword (for absolute value) and - the following *operators*: **<, >, <=, >=, ==, !=**. For example, ``prop z >= 5.0`` - selects all atoms with z coordinate greater than 5.0; ``prop abs z <= 5.0`` - selects all atoms within -5.0 <= z <= 5.0. + selects atoms based on position, using *property* **x**, **y**, + or **z** coordinate. Supports the **abs** keyword (for absolute + value) and the following *operators*: **<, >, <=, >=, ==, !=**. + For example, ``prop z >= 5.0`` selects all atoms with z + coordinate greater than 5.0; ``prop abs z <= 5.0`` selects all + atoms within -5.0 <= z <= 5.0. + sphzone *radius* *selection* + Selects all atoms that are within *radius* of the center of + geometry of *selection* + sphlayer *inner radius* *outer radius* *selection* + Similar to sphzone, but also excludes atoms that are within + *inner radius* of the selection COG **Connectivity** byres *selection* selects all atoms that are in the same segment and residue as selection, e.g. specify the subselection after the byres keyword + bonded *selection* + selects all atoms that are bonded to selection + eg: ``select name H bonded name O`` selects only hydrogens + bonded to oxygens **Index** bynum *index-range* selects all atoms within a range of (1-based) inclusive indices, - e.g. ``bynum 1`` selects the first atom in the universe; ``bynum 5:10`` - selects atoms 5 through 10 inclusive. All atoms in the - :class:`MDAnalysis.Universe` are consecutively numbered, and the index - runs from 1 up to the total number of atoms. + e.g. ``bynum 1`` selects the first atom in the universe; + ``bynum 5:10`` selects atoms 5 through 10 inclusive. All atoms + in the :class:`MDAnalysis.Universe` are consecutively numbered, + and the index runs from 1 up to the total number of atoms. **Preexisting selections** group *group-name* - selects the atoms in the :class:`AtomGroup` passed to the function as an - argument named *group-name*. Only the atoms common to *group-name* and the - instance :meth:`~select_atoms` was called from will be considered. - *group-name* will be included in the parsing just by comparison of atom indices. - This means that it is up to the user to make sure they were defined in an - appropriate :class:`Universe`. + selects the atoms in the :class:`AtomGroup` passed to the + function as an argument named *group-name*. Only the atoms + common to *group-name* and the instance :meth:`~select_atoms` + was called from will be considered. *group-name* will be + included in the parsing just by comparison of atom indices. + This means that it is up to the user to make sure they were + defined in an appropriate :class:`Universe`. fullgroup *group-name* - just like the ``group`` keyword with the difference that all the atoms of - *group-name* are included. The resulting selection may therefore have atoms - that were initially absent from the instance :meth:`~select_atoms` was - called from. - + just like the ``group`` keyword with the difference that all the + atoms of *group-name* are included. The resulting selection may + therefore have atoms that were initially absent from the + instance :meth:`~select_atoms` was called from. .. versionchanged:: 0.7.4 Added *resnum* selection. .. versionchanged:: 0.8.1 Added *group* and *fullgroup* selections. + .. versionchanged:: 0.13.0 + Added *bonded* selection """ - from . import Selection # can ONLY import in method, otherwise cyclical import! + # can ONLY import in method, otherwise cyclical import! + from . import Selection atomgrp = Selection.Parser.parse(sel, selgroups).apply(self) if len(othersel) == 0: diff --git a/package/MDAnalysis/core/Selection.py b/package/MDAnalysis/core/Selection.py index 09940bcc2f8..646ecd1c6e2 100644 --- a/package/MDAnalysis/core/Selection.py +++ b/package/MDAnalysis/core/Selection.py @@ -215,17 +215,20 @@ def _apply_KDTree(self, group): return set(res_atoms) def _apply_distmat(self, group): - sel_atoms = self.sel._apply(group) # group is wrong, should be universe (?!) + # group is wrong, should be universe (?!) + sel_atoms = self.sel._apply(group) sel_CoG = AtomGroup(sel_atoms).center_of_geometry() - sys_atoms_list = [a for a in (self._group_atoms)] # list needed for back-indexing + sel_CoG = sel_CoG.reshape(1, 3).astype(np.float32) + # list needed for back-indexing + sys_atoms_list = [a for a in (self._group_atoms)] sys_ag = AtomGroup(sys_atoms_list) - sel_CoG_str = \ - str("point ") +\ - str(sel_CoG[0]) + " " + str(sel_CoG[1]) + " " + str(sel_CoG[2]) + " " +\ - str(self.exRadius) + " and not point " + str(sel_CoG[0]) + " " + str(sel_CoG[1]) + " " + \ - str(sel_CoG[2]) + " " + str(self.inRadius) - sel = sys_ag.select_atoms(sel_CoG_str) - res_atoms = AtomGroup(set(sel)) + + box = sys_ag.dimensions if self.periodic else None + d = distances.distance_array(sel_CoG, + sys_ag.positions, + box=box) + idx = np.where((d < self.exRadius) & (d > self.inRadius))[1] + res_atoms = sys_ag[idx] return set(res_atoms) @@ -254,14 +257,20 @@ def _apply_KDTree(self, group): return set(res_atoms) def _apply_distmat(self, group): - sel_atoms = self.sel._apply(group) # group is wrong, should be universe (?!) + # group is wrong, should be universe (?!) + sel_atoms = self.sel._apply(group) sel_CoG = AtomGroup(sel_atoms).center_of_geometry() - sys_atoms_list = [a for a in (self._group_atoms)] # list needed for back-indexing + sel_CoG = sel_CoG.reshape(1, 3).astype(np.float32) + # list needed for back-indexing + sys_atoms_list = [a for a in (self._group_atoms)] sys_ag = AtomGroup(sys_atoms_list) - sel_CoG_str = str("point ") + str(sel_CoG[0]) + " " + str(sel_CoG[1]) + " " + str(sel_CoG[2]) + " " + str( - self.cutoff) - sel = sys_ag.select_atoms(sel_CoG_str) - res_atoms = AtomGroup(set(sel)) + + box = sys_ag.dimensions if self.periodic else None + d = distances.distance_array(sel_CoG, + sys_ag.positions, + box=box) + idx = np.where(d < self.cutoff)[1] + res_atoms = sys_ag[idx] return set(res_atoms) diff --git a/testsuite/MDAnalysisTests/test_atomselections.py b/testsuite/MDAnalysisTests/test_atomselections.py index e6862ebbf91..5bea7974a3b 100644 --- a/testsuite/MDAnalysisTests/test_atomselections.py +++ b/testsuite/MDAnalysisTests/test_atomselections.py @@ -437,111 +437,93 @@ class TestDistanceSelections(object): Cylindrical methods don't use KDTree """ + methods = [('kdtree', False), + ('distmat', True), + ('distmat', False)] + def setUp(self): self.u = mda.Universe(TRZ_psf, TRZ) def tearDown(self): del self.u - def choosemeth(self, sel, meth): + def choosemeth(self, sel, meth, periodic): """hack in the desired apply method""" if meth == 'kdtree': sel._apply = sel._apply_KDTree elif meth == 'distmat': sel._apply = sel._apply_distmat + + if periodic: + sel.periodic = True + else: + sel.periodic = False + return sel - def _check_around(self, meth): + def _check_around(self, meth, periodic): sel = Parser.parse('around 5.0 resid 1', self.u.atoms) - sel = self.choosemeth(sel, meth) + sel = self.choosemeth(sel, meth, periodic) result = sel.apply(self.u.atoms) r1 = self.u.select_atoms('resid 1') cog = r1.center_of_geometry().reshape(1, 3) - # Ugly hack - # KDTree doesn't do periodicity - # distmat does - if meth == 'distmat': - box = self.u.dimensions - else: - box = None - - d = distance_array(self.u.atoms.positions, r1.positions, box=box) + box = self.u.dimensions if periodic else None + d = distance_array(self.u.atoms.positions, r1.positions, + box=box) ref = set(np.where(d < 5.0)[0]) # Around doesn't include atoms from the reference group ref.difference_update(set(r1.indices)) assert_(ref == set(result.indices)) - def test_around(self): - for meth in ['kdtree', 'distmat']: - yield self._check_around, meth + for meth, periodic in self.methods: + yield self._check_around, meth, periodic - def _check_spherical_layer(self, meth): + def _check_spherical_layer(self, meth, periodic): sel = Parser.parse('sphlayer 2.4 6.0 resid 1' , self.u.atoms) - sel = self.choosemeth(sel, meth) + sel = self.choosemeth(sel, meth, periodic) result = sel.apply(self.u.atoms) r1 = self.u.select_atoms('resid 1') cog = r1.center_of_geometry().reshape(1, 3) - # Ugly hack - # KDTree doesn't do periodicity - # distmat does - if meth == 'distmat': - box = self.u.dimensions - else: - box = None - + box = self.u.dimensions if periodic else None d = distance_array(self.u.atoms.positions, cog, box=box) ref = set(np.where((d > 2.4) & (d < 6.0))[0]) assert_(ref == set(result.indices)) def test_spherical_layer(self): - for meth in ['kdtree', 'distmat']: - yield self._check_spherical_layer, meth + for meth, periodic in self.methods: + yield self._check_spherical_layer, meth, periodic - def _check_spherical_zone(self, meth): + def _check_spherical_zone(self, meth, periodic): sel = Parser.parse('sphzone 5.0 resid 1', self.u.atoms) - sel = self.choosemeth(sel, meth) + sel = self.choosemeth(sel, meth, periodic) result = sel.apply(self.u.atoms) r1 = self.u.select_atoms('resid 1') cog = r1.center_of_geometry().reshape(1, 3) - # Ugly hack - # KDTree doesn't do periodicity - # distmat does - if meth == 'distmat': - box = self.u.dimensions - else: - box = None - + box = self.u.dimensions if periodic else None d = distance_array(self.u.atoms.positions, cog, box=box) ref = set(np.where(d < 5.0)[0]) assert_(ref == set(result.indices)) def test_spherical_zone(self): - for meth in ['kdtree', 'distmat']: - yield self._check_spherical_zone, meth + for meth, periodic in self.methods: + yield self._check_spherical_zone, meth, periodic - def _check_point(self, meth): + def _check_point(self, meth, periodic): sel = Parser.parse('point 5.0 5.0 5.0 3.0', self.u.atoms) - sel = self.choosemeth(sel, meth) + sel = self.choosemeth(sel, meth, periodic) result = sel.apply(self.u.atoms) - # Ugly hack - # KDTree doesn't do periodicity - # distmat does - if meth == 'distmat': - box = self.u.dimensions - else: - box = None - + box = self.u.dimensions if periodic else None d = distance_array(np.array([[5.0, 5.0, 5.0]], dtype=np.float32), self.u.atoms.positions, box=box) @@ -550,9 +532,8 @@ def _check_point(self, meth): assert_(ref == set(result.indices)) def test_point(self): - # Doesn't work with KDTree currently - for meth in ['kdtree', 'distmat']: - yield self._check_point, meth + for meth, periodic in self.methods: + yield self._check_point, meth, periodic class TestTriclinicSelections(object):