Issue 2046 lib distances rework continued by zemanj · Pull Request #2114 · MDAnalysis/mdanalysis

zemanj · 2018-10-19T21:14:24Z

Fixes (partially) #2046. This is the fifth of a series of related PRs following PRs #2048, #2062, #2070, and #2083.

Changes made in this Pull Request:

Moved lib.distances._check_box() to lib.util.check_box()
Ensured that all methods used in capped_distance() have the same cut-off criterion (distances <= max_cutoff). According to the scipy.spatial.cKDTree.query_ball_tree docs and according to some tests I ran, this should also be the case for the pKDTree method. In practice, this is rather pathological (comparing floating point values for equality, you know...) but should still be consistent.
Ensured that none of the functions in lib.distances and lib.nsgrid make arbitrary assumptions about when a box is orthogonal. Since the number 90 has an exact float representation, we should refrain from defining arbitrary margins around that value. This is definitely a case where comparing floating point values for equality makes perfect sense (as it has always been done by lib.mdamath.triclinic_vectors and the box checks in lib.distances (now in lib.utils)).
Just like all other distance-related functions, lib.nsgridnow uses double precision for distance computation.
Prefixed all non-user-level classes in lib.nsgrid with an underscore.
Adjusted all affected tests accordingly.
Fixed up half-broken docs of lib.nsgrid.
Minor doc fix-ups in lib.utils I stumbled accross.

PR Checklist

Tests?
Docs?
CHANGELOG updated?
Issue raised/referenced?

codecov · 2018-10-19T23:24:56Z

Codecov Report

Merging #2114 into develop will increase coverage by <.01%.
The diff coverage is 100%.

@@             Coverage Diff             @@
##           develop    #2114      +/-   ##
===========================================
+ Coverage    89.41%   89.41%   +<.01%     
===========================================
  Files          157      157              
  Lines        18765    18766       +1     
  Branches      2711     2711              
===========================================
+ Hits         16778    16779       +1     
  Misses        1382     1382              
  Partials       605      605

Impacted Files	Coverage Δ
package/MDAnalysis/lib/util.py	`88.31% <100%> (+0.13%)`	⬆️
package/MDAnalysis/lib/distances.py	`97.4% <100%> (-0.06%)`	⬇️

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update e08e817...42b856a. Read the comment docs.

zemanj · 2018-10-22T09:00:32Z

testsuite/MDAnalysisTests/lib/test_util.py

        assert line == line.lstrip()
+
+
+class TestCheckBox(object):


I just copied this over from lib/test_distances.py. The code itself remained unchanged.

zemanj · 2018-10-22T09:02:43Z

testsuite/MDAnalysisTests/lib/test_util.py

        assert line == line.lstrip()
+
+
+class TestCheckBox(object):


I just copied this over from lib/test_distances.py. The code itself remained unchanged.

zemanj · 2018-10-22T09:02:58Z

package/MDAnalysis/lib/util.py

    return lines[0].lstrip() + "\n" + textwrap.dedent("\n".join(lines[1:]))
+
+
+def check_box(box):


I just cut & pasted this from lib/distances.py. The code itself is unchanged.

zemanj

Some explanatory comments

zemanj · 2018-10-22T09:08:57Z

package/MDAnalysis/lib/distances.py

 from .c_distances_openmp import OPENMP_ENABLED as USED_OPENMP


-def _check_box(box):


This function is not specific to distance computation and might be useful elsewhere --> moved to lib.util.check_box()

zemanj · 2018-10-22T09:14:44Z

package/MDAnalysis/lib/distances.py

    if N > 1:
        distvec = self_distance_array(reference, box=box)
-        dist = np.full((N, N), max_cutoff, dtype=np.float64)
+        dist = np.full((N, N), np.finfo(np.float64).max, dtype=np.float64)


Old fill value max_cutoff doesn't work with cutoff criterion distance <= max_cutoff.
np.finfo(np.float64).max is the largest number representable by np.float64, so this fill value is guaranteed to work for all possible values of max_cutoff

zemanj · 2018-10-22T09:16:02Z

package/MDAnalysis/lib/nsgrid.pyx

 cdef real rvec_norm2(const rvec a) nogil:
    return a[XX]*a[XX] + a[YY]*a[YY] + a[ZZ]*a[ZZ]

-cdef void rvec_clear(rvec a) nogil:


zemanj · 2018-10-22T09:18:22Z

package/MDAnalysis/lib/nsgrid.pyx

+                    # mdamath.triclinic_vectors explicitly sets the off-diagonal
+                    # elements to zero if the box is orthogonal, so we can
+                    # safely check floating point values for equality here
+                    if box[i, j] != 0.0:


This is to correspond to lib.util.check_box() and therefore all functions in lib.distances.

zemanj · 2018-10-22T09:23:36Z

package/MDAnalysis/lib/util.py

    return lines[0].lstrip() + "\n" + textwrap.dedent("\n".join(lines[1:]))
+
+
+def check_box(box):


This is just cut & pasted from lib/distances.py, no changes in the code.

zemanj · 2018-10-22T09:24:00Z

testsuite/MDAnalysisTests/lib/test_distances.py

 from MDAnalysis.lib import mdamath
 from MDAnalysis.tests.datafiles import PSF, DCD, TRIC

-class TestCheckBox(object):


Moved to lib/test_util.py

zemanj · 2018-10-22T09:24:20Z

testsuite/MDAnalysisTests/lib/test_util.py

        assert line == line.lstrip()
+
+
+class TestCheckBox(object):


This is just cut & pasted from lib/test_distances.py, no changes in the code.

zemanj

Some explanatory comments

zemanj · 2018-10-22T09:33:56Z

Sorry for the comment mess. GitHub is extremely unresponsive for me today. It took half an hour until I could finally see that all my seemingly unsuccessful attempts to add comments did in fact work. 😠
Hope you can see the comments, they only appear for me on every fifth page reload or so.

richardjgowers · 2018-10-22T15:13:14Z

package/MDAnalysis/lib/nsgrid.pyx


 ctypedef np.int_t ns_int
 ctypedef np.float32_t real
+ctypedef np.float64_t dreal


Is there a performance different for moving to double precision? (probably yes?) If so, maybe we should open an issue about precision in our distance calculations. Especially as for this Cython it's easy to used a fused type (aka template)

Good call, I'll run some benchmarks.

richardjgowers · 2018-10-22T15:14:01Z

@zemanj I'll find some time to go over this in detail, but it looks good at a glance

zemanj · 2018-10-23T12:49:48Z

@richardjgowers I tested the performance impact of switching from single to double precision distance computation in lib.nsgrid with the following code:

>>> import MDAnalysis as mda
>>> from MDAnalysisTests.datafiles import TRR, TPR
>>> from MDAnalysis.lib.distances import capped_distance, self_capped_distance
>>> u = mda.Universe(TPR, TRR)
>>> water_oxygen = u.select_atoms('resname SOL and name OW')
>>> print(water_oxygen.n_atoms)
11084
>>> print(len(u.trajectory))
10

>>> %%timeit
>>> for ts in u.trajectory:
>>>     pairs, dists = capped_distance(water_oxygen.positions,
...                                   water_oxygen.positions,
...                                   15, box=ts.dimensions,
...                                   method='nsgrid')
>>> %%timeit
>>> for ts in u.trajectory:
>>>     pairs, dists = self_capped_distance(water_oxygen.positions,
...                                         15, box=ts.dimensions,
...                                         method='nsgrid')

Output for capped_distance:
1 loop, best of 3: 24.3 s per loop (single precision)
1 loop, best of 3: 25.6 s per loop (double precision)
➡️ Performance loss: ~5%

Output for self_capped_distance:
1 loop, best of 3: 21.1 s per loop (single precision)
1 loop, best of 3: 21.9 s per loop (double precision)
➡️ Performance loss: ~4%

So I'd say that the performance impact is rather negligible.

richardjgowers

Cool thanks @zemanj

Interesting that the performance loss is so small, but I guess a good thing for now

richardjgowers · 2018-10-24T21:59:13Z

Looks like it needs a rebase onto develop

…d_distances() methods to <= max_cutoff

…cutoff; added underscore prefix to non-user classes

zemanj · 2018-10-25T06:44:09Z

@richardjgowers I think the small performance impact indicates one or both of the following:

the algorithm is strongly memory-bound
the implementation might have some inefficiencies

I guess it's a bit of both since I already identified some rather crude thingies in nsgrid.pyx (such as the way neighboring cells are determined). However, I'll need to do some profiling before I can tell more.

zemanj · 2018-10-25T11:13:25Z

I recently noticed an increasing amount of CI hiccups where jobs fail due to errors unrelated to the actual PR. Usually, those are I/O-related errors. I started collecting these cases and will open an issue once I gathered enough information.

zemanj · 2018-10-25T11:14:38Z

@richardjgowers Thank you for reviewing!

richardjgowers · 2018-10-25T13:25:59Z

@zemanj yeah probably because it's memory bound.

WRT the I/O hiccups, it might be possible to fix some of these using longer scoped fixtures.... There's currently the encore HES thing and badzipfile I can think of :)

richardjgowers self-assigned this Oct 19, 2018

zemanj commented Oct 22, 2018

View reviewed changes

richardjgowers reviewed Oct 22, 2018

View reviewed changes

richardjgowers approved these changes Oct 24, 2018

View reviewed changes

richardjgowers added this to the 0.19.x milestone Oct 24, 2018

zemanj added 5 commits October 25, 2018 07:49

moved check_box() to lib.util; changed cutoff criterion of all *cappe…

5518d0b

…d_distances() methods to <= max_cutoff

double precision distances in nsgrid; changed cutoff criterion to <= …

225dfc0

…cutoff; added underscore prefix to non-user classes

changed tests accordingly

767c5b2

nsgrid docs fixup

145490d

minor doc fix-ups for lib.util

42b856a

zemanj force-pushed the issue-2046-lib-distances-rework branch from 81e3f5a to 42b856a Compare October 25, 2018 05:52

zemanj merged commit c9a37e3 into MDAnalysis:develop Oct 25, 2018

zemanj deleted the issue-2046-lib-distances-rework branch October 25, 2018 11:15

zemanj mentioned this pull request Oct 25, 2018

MDAnalysis.lib.distances needs rework #2046

Open

51 tasks

		return lines[0].lstrip() + "\n" + textwrap.dedent("\n".join(lines[1:]))


		def check_box(box):

		from .c_distances_openmp import OPENMP_ENABLED as USED_OPENMP


		def _check_box(box):

Conversation

zemanj commented Oct 19, 2018 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

PR Checklist

Uh oh!

codecov bot commented Oct 19, 2018 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

zemanj left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

zemanj left a comment

Choose a reason for hiding this comment

Uh oh!

zemanj commented Oct 22, 2018

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

richardjgowers commented Oct 22, 2018

Uh oh!

zemanj commented Oct 23, 2018

Uh oh!

richardjgowers left a comment

Choose a reason for hiding this comment

Uh oh!

richardjgowers commented Oct 24, 2018

Uh oh!

zemanj commented Oct 25, 2018

Uh oh!

zemanj commented Oct 25, 2018

Uh oh!

zemanj commented Oct 25, 2018

Uh oh!

richardjgowers commented Oct 25, 2018

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

zemanj commented Oct 19, 2018 •

edited

Loading

codecov bot commented Oct 19, 2018 •

edited

Loading