Encore tests failing due to float precision

[richardjgowers]

I'm mostly just going to steal @jbarnoud 's debugging work here, so we can talk about the failures in their own issue...

The value error seems unrelated to the numpy version, indeed. analysis.encore.cutils.PureRMSD expects float64_t coordinates but receives float ones. The first commit to fail the test is e356604 (see travis: https://travis-ci.org/MDAnalysis/mdanalysis/builds/192119548). The pull request (#1132) was passing the tests, but it got merged right after #1145 in 789a96c. The two pull requests are touching parts of ENCORE. My guess is that the MemoryReader is now working with float32 rather than float64 which breaks analysis.encore.cutils.PureRMSD. The simple solution would be to change analysis.encore.cutils.PureRMSD to use float32 instead of float64 but I would like to know what @kain88-de and @wouterboomsma think about it before changing this.

And my 2c, is that regardless of what MemoryReader is storing/supplying, there should probably be a check and/or np.asarray call somewhere to avoid this.

[jbarnoud]

A fancier way to solve the issue would be to have an intermediate layer on top of analysis.encore.cutils.PureRMSD that checks what type the function is called with and then call the float32 or float64 version of analysis.encore.cutils.PureRMSD. I do not know cython, but there may even be a canonical way of solving that issue.

[wouterboomsma]

Yes, I missed this aspect when changing MemoryReader to always use float32. I guess it should be safe to convert the pureRMSD calculation to use float32 as well (at least we can try, and see how many of the encore tests fail). Ultimately, I agree that it would be nice to have support for both float64 and float32, but since float64 does not seem to be generally supported in MDAnalysis currently, it's probably not worth putting in the effort to put in an intermediate layer at this point.

[kain88-de]

I would also like to have it float32. Otherwise the change of dtype will cause a copy of the array in memory and the computation might not succeed because we run out of memory.

[richardjgowers]

http://cython.readthedocs.io/en/latest/src/userguide/fusedtypes.html

@jbarnoud if we want proper float32/float64 support this might be a way to do it

[jbarnoud]

This happens to be th exact page I am reading now ;)

…

On 17-01-17 14:53, Richard Gowers wrote: http://cython.readthedocs.io/en/latest/src/userguide/fusedtypes.html @jbarnoud <https://github.com/jbarnoud> if we want proper float32/float64 support this might be a way to do it — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#1168 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/ABUWujxYdjWnfMoyIQn1CzE1quu8ep3Iks5rTMf0gaJpZM4Lle21>.

[mtiberti]

Thanks for the comments - I agree it would be better to make it float32 as to be consistent with the rest, and then revise it at a later time if we need to introduce support for float64 as well. I'll get on it and let's see how it works.

[jbarnoud]

I investigated the problem a bit. analysis.encore.cutils.PureRMSD is what complains, but it is called by analysis.encore.condistmatrix.set_rmsd_matrix_elements that has two branches.

The first branch is calling PureRMSD without paying attention to the type of the coordinates, this branch is visited when there are not fit coordinates provided, and it is the branch that currently causes the error.

The second branch is visited when there are fir coordinates provided. There the coordinates can be float32 or float64 but they are explicitly casted to float64 before they are fed to PureRMSD.

Having everything homogenized to float32 (with explicit cast to float32 in both branched, and PureRMSD modified) does not generate errors anymore, but one test fails:

======================================================================
FAIL: test_clustering_three_ensembles_two_identical (MDAnalysisTests.analysis.test_encore.TestEncoreClustering)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/home/jon/Envs/oldnp/local/lib/python2.7/site-packages/MDAnalysisTests/analysis/test_encore.py", line 457, in test_clustering_three_ensembles_two_identical
    err_msg="Unexpected result: {0}".format(cluster_collection))
  File "/home/jon/Envs/oldnp/local/lib/python2.7/site-packages/numpy/testing/utils.py", line 381, in assert_equal
    raise AssertionError(msg)
AssertionError: 
Items are not equal:
Unexpected result: 0 (size:2,centroid:2): array([ 2, 43])
1 (size:2,centroid:3): array([ 3, 44])
2 (size:1,centroid:4): array([4])
3 (size:1,centroid:5): array([5])
4 (size:1,centroid:6): array([6])
5 (size:1,centroid:7): array([7])
6 (size:1,centroid:8): array([8])
7 (size:1,centroid:9): array([9])
8 (size:1,centroid:10): array([10])
9 (size:2,centroid:11): array([11, 52])
10 (size:2,centroid:12): array([12, 53])
11 (size:1,centroid:13): array([13])
12 (size:1,centroid:14): array([14])
13 (size:2,centroid:15): array([15, 56])
14 (size:2,centroid:16): array([16, 57])
15 (size:1,centroid:17): array([17])
16 (size:3,centroid:19): array([19, 40, 60])
17 (size:3,centroid:23): array([22, 23, 24])
18 (size:2,centroid:26): array([25, 26])
19 (size:2,centroid:28): array([27, 28])
20 (size:3,centroid:30): array([29, 30, 31])
21 (size:3,centroid:33): array([32, 33, 34])
22 (size:3,centroid:36): array([35, 36, 37])
23 (size:3,centroid:41): array([ 0, 20, 41])
24 (size:3,centroid:42): array([ 1, 21, 42])
25 (size:1,centroid:45): array([45])
26 (size:1,centroid:46): array([46])
27 (size:1,centroid:47): array([47])
28 (size:1,centroid:48): array([48])
29 (size:1,centroid:49): array([49])
30 (size:1,centroid:50): array([50])
31 (size:1,centroid:51): array([51])
32 (size:1,centroid:54): array([54])
33 (size:1,centroid:55): array([55])
34 (size:1,centroid:58): array([58])
35 (size:4,centroid:59): array([18, 38, 39, 59])

 ACTUAL: 36
 DESIRED: 40

Everything homogenized to float64 has every tests passing.

I did not try the fused types, and I won't since @mtiberti stepped in 👍

[mtiberti]

Thanks @jbarnoud - I tried doing the same and got the same behaviour, and it also happens using the floating cython fused type for floating point values. Investigating the issue

[orbeckst]

Small sidenote on the float32/float64 question. When @rbrtdlgd optimized the lib.qcprot RMSD calculation for use with MDAnalysis, he found that switching to 32 bit produced unacceptable errors in the RMSD (PR #930 and https://dx.doi.org/10.6084/m9.figshare.3823293.v1 ).

[jbarnoud]

@orbeckst Isn't the solution to just make sure PureRMSD is always eating float64, then?

[wouterboomsma]

@orbeckst, I don't think we would have the same problem here, since the PureRMSD does not involve a superimposition. I think float32 should suffice - but @mtiberti is looking into it right now.

[richardjgowers]

If this is proving tricky to fix, can we tag the failing tests as known failures so every PR doesn't have failing tests?

[mtiberti]

I had a look yesterday but haven't come up with a solution yet - the differences in RMSD between calculating them with float32 and float64 are small as all of you suggested, but that's apparently enough to change the clustering outcomes in that specific case (all the other tests involving clustering pass). So for the moment I would suggest to change the expected output as the current so that at least other builds pass, and we can have a deeper look in the meanwhile if you think it's worth it.

[kain88-de]

What clustering algorithm do you use? I would be interested in the difference you see if you are able to share the dataset. Also just sharing the code in a notebook would be nice. This is just for my own curiosity for numerical stability.

[jbarnoud]

So for the moment I would suggest to change the expected output as the current so that at least other builds pass, and we can have a deeper look in the meanwhile if you think it's worth it.

Since casting the input for PureRMSD to float64 systematically seems to work, couldn't we do that as a short term solution rather than altering the tests?

[mtiberti]

Sure, that's also a viable option.

@kain88-de I'll get back to you - for now I've been switching manually between float32 and float64. We are using Affinity Propagation.

[kain88-de]

Thanks @miterbi

Also I'm for @jbarnoud solution to do a cast to float64 for now to fix the bug until we understand the numerical implications of using single precision.