Accept lists for bond indices, raise warning/error when add_TopologyAttr gets something unusable

[lilyminium]

Expected behavior

It is possible to add bonds to a bond-less Universe (or replace bonds in a bonded Universe) with u.add_TopologyAttr in a very specific way. I would expect errors or warnings to pop up when u.add_TopologyAttris passed something it can't work with. I would also not expect lists of indices to act so differently from tuples of indices.

Actual behavior

No warnings or errors, until the u.bonds attribute is accessed and either returns an AttributeError or an empty TopologyGroup.

Code to reproduce the behavior

I was trying to make a Universe of water.

>>> import MDAnalysis as mda
>>> import numpy as np
>>> sol = mda.Universe.empty(3000,
                         n_residues=1000, 
                         atom_resindex=np.repeat(range(1000), 3), 
                         residue_segindex=[0]*1000,
                         trajectory=True)
>>> sol
<Universe with 3000 atoms>
>>> sol.add_TopologyAttr('name', ['O', 'H1', 'H2']*n_residues)
>>> sol.atoms.names
array(['O', 'H1', 'H2', ..., 'O', 'H1', 'H2'], dtype=object)

Attempts

1. Trying a list of lists of indices (bonds attribute is available but broken)

>>> bonds = []
>>> for o in range(0, n_atoms, 3): 
...        bonds.extend([[o, o+1], [o, o+2]])
...
>>> sol.add_TopologyAttr('bonds', bonds)
>>> sol.bonds
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
MDAnalysis/core/topologyattrs.py in get_atoms(self, ag)
   1692             unique_bonds = set(itertools.chain(
-> 1693                 *[self._bondDict[a] for a in ag.ix]))
   1694         except TypeError:

TypeError: unhashable type: 'list'

During handling of the above exception, another exception occurred:

TypeError                                 Traceback (most recent call last)
<ipython-input-239-076c3c5d8d4d> in <module>
----> 1 sol.bonds

MDAnalysis/core/universe.py in bonds(self)
    702     def bonds(self):
    703         """Bonds between atoms"""
--> 704         return self.atoms.bonds
    705 
    706     @property

MDAnalysis/core/groups.py in getter(self)
    257         """
    258         def getter(self):
--> 259             return attr.__getitem__(self)
    260 
    261         def setter(self, values):

MDAnalysis/core/topologyattrs.py in __getitem__(self, group)
    270         """Accepts an AtomGroup, ResidueGroup or SegmentGroup"""
    271         if isinstance(group, (Atom, AtomGroup)):
--> 272             return self.get_atoms(group)
    273         elif isinstance(group, (Residue, ResidueGroup)):
    274             return self.get_residues(group)

MDAnalysis/core/topologyattrs.py in get_atoms(self, ag)
   1694         except TypeError:
   1695             # maybe we got passed an Atom
-> 1696             unique_bonds = self._bondDict[ag.ix]
   1697         bond_idx, types, guessed, order = np.hsplit(
   1698             np.array(sorted(unique_bonds)), 4)

TypeError: unhashable type: 'numpy.ndarray'

2. Trying a TopologyGroup (doesn't add bonds)

>>> bondgroup = mda.core.topologyobjects.TopologyGroup(np.array(bonds), sol)
>>> bondgroup
<TopologyGroup containing 2000 bonds>
>>> sol.add_TopologyAttr('bonds', bondgroup)
>>> sol.bonds
---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
MDAnalysis/core/universe.py in __getattr__(self, key)
    509         try:
--> 510             segment = self._instant_selectors[key]
    511         except KeyError:

KeyError: 'bonds'

During handling of the above exception, another exception occurred:

AttributeError                            Traceback (most recent call last)
<ipython-input-243-076c3c5d8d4d> in <module>
----> 1 sol.bonds

MDAnalysis/core/universe.py in __getattr__(self, key)
    510             segment = self._instant_selectors[key]
    511         except KeyError:
--> 512             raise AttributeError('No attribute "{}".'.format(key))
    513         else:
    514             warnings.warn("Instant selector Universe.<segid> "

AttributeError: No attribute "bonds".

3. Trying an atom-wise list of lists of tuples of indices (adds empty group)

>>> bonds = []
>>> for o in range(0, n_atoms, 3):
...     bonds.extend([[(o, o+1), (o, o+2)],
...                   [(o, o+1)],
...                   [(o, o+2)]])
...
>>> sol.add_TopologyAttr('bonds', bonds)
>>> sol.bonds
<TopologyGroup containing 0 bonds>

4. Trying a list of tuples of indices (works!)

>>> bonds = []
>>> for o in range(0, n_atoms, 3):
...     bonds.extend([(o, o+1), (o, o+2)])
...
>>> sol.add_TopologyAttr('bonds', bonds)
>>> sol.bonds
<TopologyGroup containing 2000 bonds>

🌟Bonus 5. Trying to raise any kind of error🌟
Turns out add_TopologyAttr will silently accept any iterable, and produce empty TopologyGroups for any iterable of something subscriptable.

>>> sol.add_TopologyAttr('bonds', 'invalid bond')
>>> sol.bonds
<TopologyGroup containing 0 bonds>
>>> sol.add_TopologyAttr('bonds', ['0'] * 10)
>>> sol.bonds
<TopologyGroup containing 0 bonds>

Again, errors only pop up when sol.bonds is called, when I pass an iterable of something not subscriptable.

Comments
Approach 1 and 4 seem very similar to me. Would it be possible to add another except clause that tries to convert it to a tuple here:

mdanalysis/package/MDAnalysis/core/topologyattrs.py

Lines 1690 to 1696 in 9372d09

	def get_atoms(self, ag):
	try:
	unique_bonds = set(itertools.chain(
	*[self._bondDict[a] for a in ag.ix]))
	except TypeError:
	# maybe we got passed an Atom
	unique_bonds = self._bondDict[ag.ix]

Otherwise it would be preferable to raise errors at the add_TopologyAttr stage rather than when the bonds attribute is called. (In the case of Approaches 3 and 5, no errors are ever raised).

Currently version of MDAnalysis

• Which version are you using? (run python -c "import MDAnalysis as mda; print(mda.__version__)") 0.20.1
• Which version of Python (python -V)? 3.7.x
• Which operating system? MacOS Mojave

[lilyminium]

Also: TopologyGroup doesn't accept a list for the bond indices. It must be a numpy array. If it is ever likely that a user might try to instantiate a TopologyGroup, numpy.asarray may be helpful?

>>> bondgroup = mda.core.topologyobjects.TopologyGroup(bonds, sol) # when bonds is a list of lists
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-266-651c19dc9bd9> in <module>
----> 1 bondgroup = mda.core.topologyobjects.TopologyGroup(bonds, sol)

MDAnalysis/core/topologyobjects.py in __init__(self, bondidx, universe, btype, type, guessed, order)
    553 
    554         if nbonds > 0:
--> 555             uniq, uniq_idx = util.unique_rows(bondidx, return_index=True)
    556 
    557             self._bix = uniq

MDAnalysis/lib/util.py in unique_rows(arr, return_index)
   1581     # this can occur when second dimension was created through broadcasting
   1582     # eg: idx = np.array([1, 2])[None, :]
-> 1583     if not arr.flags['OWNDATA']:
   1584         arr = arr.copy()
   1585 

AttributeError: 'list' object has no attribute 'flags'

[orbeckst]

I think we had decided that a bond is internally represented by a tuple that is ordered in increasing atom index order (could have taken a set, I suppose). Using a tuple makes it hashable and is also the fastest container in Python (IIRC).

It seems that the problem is that although topology building is very useful, it's currently operating close to the internals of MDA, with very little shielding. A future project could build a user-friendly API that does the more obvious things by providing add_bond(Atom1, Atom2) and add_bonds(AG1, AG2) (for len(AG1) == len(AG2)).

@richardjgowers @jbarnoud did I correctly summarize? Is it worthwhile making add_TopologyAttr more userfriendly?

[jbarnoud]

I guess we should indeed make add_TopologyAttr fail early, and maybe deal with the tuple conversion. I would be happy with a user-facing API to manipulate bonds, but indeed it would require to change some things at the universe level. From the top of my head, we at least need a way to refresh the cache for the fragments.

[richardjgowers]

Yeah I think they’re strictly tuples for hashing reasons, even n x 2 arrays are a pain (not impossible) to hash for unique rows. But yeah, can likely be made more usable and probably better in most ways.

…

On Sat, Oct 26, 2019 at 09:33, Jonathan Barnoud ***@***.***> wrote: I guess we should indeed make add_TopologyAttr fail early, and maybe deal with the tuple conversion. I would be happy with a user-facing API to manipulate bonds, but indeed it would require to change some things at the universe level. From the top of my head, we at least need a way to refresh the cache for the fragments. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#2373?email_source=notifications&email_token=ACGSGB6ZPXSBLCQN2PAWZ2LQQP6ELA5CNFSM4JE5KIBKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOECKDBHA#issuecomment-546582684>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACGSGB64J4AJ3KV5AWFGU23QQP6ELANCNFSM4JE5KIBA> .

[lilyminium]

I was more hoping that lists of lists could be converted to lists of tuples when bonds are instantiated, as in #2379 . Making a more user-friendly API seems like a bigger project.
del u._cache['fragments'] will refresh the cache for adding an entire new group (u.add_TopologyAttr('bonds', bonds)) as a workaround for now?