Set operations for Groups

[khuston]

I apologize if this is already a feature or if has been discussed before, but I've often thought AtomGroup subtraction would be a convenient feature, in the same way __add__ is already overridden for AtomGroup. It would be defined in the same way set subtraction is defined for Python:

For AtomGroup instances:

return AtomGroup(list(set(self._atoms) - set(other._atoms)))

For Atom instances:

return AtomGroup(list(set(self._atoms) - set([other])))

It's a simple addition, but it seems to work. If it's alright, I can submit a pull request with the feature and tests.
Example usage:

In [2]: u.select_atoms('resname 58')
Out[2]: <AtomGroup with 300 atoms>

In [3]: u.select_atoms('resname 58 and prop z < 220')
Out[3]: <AtomGroup with 201 atoms>

In [4]: u.select_atoms('resname 58') - u.select_atoms('resname 58 and prop z < 220')
Out[4]: <AtomGroup with 99 atoms>

[richardjgowers]

Hey Kyle,

Yeah I think this would be great. I think @jbarnoud ? recently proposed this somewhere else, but I can't find it for the life of me.

We're changing how AtomGroups work quite soon, but I think most of the hard work here is coming up with solid tests (which won't change) rather than the actual AtomGroup methods (which will change). So yeah, put together something

[jbarnoud]

I have a work in progress based on the new topology scheme. I implemented all the set operators and I am working at testing them.

The current status is visible at https://github.com/jbarnoud/mdanalysis/commits/issue-363-set-operators.

[dotsdl]

Since AtomGroup is not technically a set, how are repeats handled? Are repeats just dropped and the result just sorted?

[mnmelo]

As a workaround while we decide and develop you can use:

# to emulate difference = grp1 - grp2
difference = grp1.select_atoms("not group other", other=grp2)

[jbarnoud]

@dotsdl Yes. When using set operations, the repeats are dropped in the result, and the index are sorted.

[khuston]

Thank you for that @mnmelo - I didn't know about the other keyword.

So since __eq__ for Atoms is based on the index, atoms with the same index will be dropped. I forgot about the case where a user adds AtomGroups together to get an AtomGroup with duplicate atoms.

[mnmelo]

As a clarification on the workaround: group is the keyword. other is an arbitrary name that must match a named argument to select_atoms (trhough which you pass an existing AtomGroup). You might as well write:

difference = grp1.select_atoms("not group khuston", khuston=grp2)

[richardjgowers]

One option here is to not use the - operator, but only allow it via a difference method. This way it's more explicit that a set difference operation has been done, rather than another interpretation of subtraction.

[jbarnoud]

I really like the idea of having all set operators available, and not only the corresponding functions. But I have to agree that, in the particular case of - may be confusing because of the asymmetry with + that do keep the duplicates.

[orbeckst]

My expectation for "-" would be that it

• keeps the original order
• removes any atom that is found in the subtracted group

i.e.,

def subtract(a, b):
    return [idx for idx in a.indices if idx in b.indices]

Admittedly, these are not set operations but AtomGroups are not sets either so I feel it's ok to be idiosyncratic. I would reserve explicit methods such as g.difference() or g.union() for actual set operations.

[richardjgowers]

@orbeckst So for example?

ag1 = [1, 2, 2, 5, 4, 3]

ag2 = [2, 4]

ag1 - ag2 = [1, 5, 3]

That makes sense to me (and is the opposite to +), and then having the set operations explicit via methods makes that clear too.

[orbeckst]

Yes, that's what makes sense to me.

Oliver Beckstein
email: [email protected]

Am Feb 23, 2016 um 3:24 schrieb Richard Gowers [email protected]:

@orbeckst So for example?

ag1 = [1, 2, 2, 5, 4, 3]

ag2 = [2, 4]

ag1 - ag2 = [1, 5, 3]
That makes sense to me (and is the opposite to +), and then having the set operations explicit via methods makes that clear too.

—
Reply to this email directly or view it on GitHub.

[richardjgowers]

@jbarnoud did some cool code here:

jbarnoud@63c1f2b

Should we add these directly as Group methods? Or have a setoperations submodule in MDAnalysis?