why is (cached) fragments slow?

[orbeckst]

Is your feature request related to a problem? Please describe.
The "make whole" (unwrap) on the fly transformation PR #2038 needs the fragments at every time step. Apparently, the fragments are recalculated every time (EDIT: They are not – see discussion below)

>>> u = mda.Universe(data.TPR, data.XTC)
>>> %timeit u.atoms.fragments
556 ms ± 71.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Recalculating fragments takes 0.5 seconds per step.

Describe the solution you'd like
There is no need to recalculate fragments as they are based on bonds. Fragments can be created when the topology is created.

Or only recalculate if the bonds change, or just tell users that if they change bonds then they also need to recalculate fragments.

EDIT: Recalculating fragments is not the problem but building the actual data structure appears to be slow.

Describe alternatives you've considered
Recalculating fragments every time the managed attribute is accessed guards against stale data, which is good. However, user expectation is that attributes behave as instantaneous value stores (and the fact that we use managed attributes should not concern users). Furthermore, I contend that the case of accessing fragments for different time steps is more common than the case of changing the topology.

Additional context

Doing the unwrap is slow no matter what, ca 4s/iteration

In [81]: u = mda.Universe(TPR, XTC)
In [82]: u.trajectory.add_transformations(MDAnalysis.transformations.wrap.unwrap(u.atoms))
In [83]: with mda.Writer("unwrap.xtc", n_atoms=u.atoms.n_atoms) as W:
    ...:     for ts in tqdm.tqdm(u.trajectory):
    ...:         W.write(u.atoms)
    ...:
100%|██████████████████████████████████████████████████████████| 10/10 [00:39<00:00,  3.90s/it]

but at least shaving off 0.5 s is a start.

[richardjgowers]

https://github.com/MDAnalysis/mdanalysis/blob/develop/package/MDAnalysis/core/universe.py#L998 I think they are cached at a universe level... so maybe something else weird is happening here. Maybe it’s the construction of the fragments for this AG from the master list?

…

On Fri, Oct 25, 2019 at 22:21, Oliver Beckstein ***@***.***> wrote: *Is your feature request related to a problem? Please describe.* The "make whole" (unwrap) on the fly transformation PR #2038 <#2038> needs the fragments at every time step. Apparently, the fragments are recalculated every time >>> u = mda.Universe(data.TPR, data.XTC) >>> %timeit u.atoms.fragments 556 ms ± 71.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each) Recalculating fragments takes 0.5 seconds per step. *Describe the solution you'd like* There is no need to recalculate fragments as they are based on bonds. Fragments can be created when the topology is created. Or only recalculate if the bonds change, or just tell users that if they change bonds then they also need to recalculate fragments. *Describe alternatives you've considered* Recalculating fragments every time the managed attribute is accessed guards against stale data, which is good. However, user expectation is that attributes behave as instantaneous value stores (and the fact that we use managed attributes should not concern users). Furthermore, I contend that the case of accessing fragments for different time steps is more common than the case of changing the topology. *Additional context* Doing the unwrap is slow no matter what, ca 4s/iteration In [83]: with mda.Writer("unwrap.xtc", n_atoms=u.atoms.n_atoms) as W: ...: for ts in tqdm.tqdm(u.trajectory): ...: W.write(u.atoms) ...: 100%|██████████████████████████████████████████████████████████| 10/10 [00:39<00:00, 3.90s/it] but at least shaving off 0.5 s is a start. — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#2376?email_source=notifications&email_token=ACGSGB65L6BPM7BKBCWNO43QQNPORA5CNFSM4JFJIAGKYY3PNVWWK3TUL52HS4DFUVEXG43VMWVGG33NNVSW45C7NFSM4HUPSTAQ>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACGSGB4XLLVXD5XAINFUZI3QQNPORANCNFSM4JFJIAGA> .

[orbeckst]

As they are cached then something else appears to be slow. Needs some line profiling, I think.

[orbeckst]

Can we cache

mdanalysis/package/MDAnalysis/core/topologyattrs.py

Line 1830 in e110108

def fragments(self):

?

ADD: The comprehension

mdanalysis/package/MDAnalysis/core/topologyattrs.py

Lines 1857 to 1858 in e110108

	return tuple(sorted(set(fragdict[aix].fragment for aix in self.ix),
	key=lambda x: x[0].ix))

looks slow for large systems.

[zemanj]

Caching AtomGroup.fragments is not straightforward since it requires an AtomGroup to know if the underlying Universe._fragdict changed.

I worked on a solution (without caching) to the performance problem a while ago based on 0.19.3, I just didn't have time to write tests and create a PR for it, plus some of the stuff I added probably needs to get thrown out again. It currently lives in my faster_unwrapping branch.

Regarding the performance of AtomGroup.fragments, I get the following:

import MDAnalysis as mda
from MDAnalysisTests.datafiles import TPR, XTC
%timeit u.atoms.fragments

Result for current dev branch (0.20.2):

328 ms ± 6.05 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Results for faster_unwrapping branch (0.19.3):

537 µs ± 5.53 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

That's 610 times faster and should make the problem go away.

Anyway, it won't help if it keeps growing mold in my MDA fork...

[orbeckst]

Maybe it would motivate you to see a massive drop in the benchmarks? (I realized that we were not timing fragments because the naming of the benchmark was incorrect – thus PR #2240 )

Making fragments faster would be an important first step to making make_whole() faster. 4s per frame is just way too slow to be really useful.

[zemanj]

Gamification FTW! 😉
Unless I missed something in the AtomGroup methods benchmarks, we still need to add AtomGroup.unwrap() there for all supported compounds.

[orbeckst]

Hmph – the benchmark https://www.mdanalysis.org/benchmarks/#ag_methods.FragmentFinding.time_find_fragments is failing. Need to look into it.