Division by 0 in alinging two lines with rotation_matrix(a, b)

rzu512 commented 6 years ago

Expected behaviour

no error

Actual behaviour

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "<python3.6 dir>/site-packages/MDAnalysis/analysis/align.py", line 276, in rotation_matrix
    rmsd = qcp.CalcRMSDRotationalMatrix(a, b, N, rot, weights)
  File "MDAnalysis/lib/qcprot.pyx", line 286, in MDAnalysis.lib.qcprot.CalcRMSDRotationalMatrix (MDAnalysis/lib/qcprot.c:2764)
  File "MDAnalysis/lib/qcprot.pyx", line 387, in MDAnalysis.lib.qcprot.FastCalcRMSDAndRotation (MDAnalysis/lib/qcprot.c:3625)
ZeroDivisionError: float division

Code to reproduce the behaviour

import numpy as np
from MDAnalysis.analysis.align import rotation_matrix
dtype = np.float32  # np.float64 gives same error
a = np.array([[0.13662023, -0.39526145, 0.75835566],
              [-0.13662023, 0.39526145, -0.75835566]],
             dtype=dtype)
b = np.array([[-0.30452963, 0.78899221, -0.18642158],
              [0.30452963, -0.78899221, 0.18642158]],
             dtype=dtype)
R, rmsd = rotation_matrix(a, b)  # division by zero

Current version of MDAnalysis:

$ python3 -c "import MDAnalysis as mda; print(mda.__version__)"
Warning: #####
MDAnalysis on python 3 is highly experimental!
It is mostly non functional and dramatically untested.
Use at your own risks!!!

  ''')
0.16.2

I tried python 2, and the test code gives the same error.

kain88-de commented 6 years ago

This is still present in develop.

tylerjereddy commented 6 years ago

Observations:

the Newton-Raphson method code block where this happens appears to be identical in the C & Cython versions of the code
changing just one of the x coordinates in just one of the arrays by 0.0000001 alleviates the division by 0 Error

Here's some debug values from the Cython Newton's Method loop using the above input:

('x2:', 2.250000172635994)
('mxEigenV:', 1.5000000575453303)
('b:', -3.375000388430994)
('a:', -3.375000388430994)
('first term:', 6.750000776861988)
('remaining terms:', -6.750000776861988)
('numerator:', -8.881784197001252e-16)

Some more observations based on the above debug output:

that's a lot of floating point values being printed !!
it is clear how the zero in the denominator arises based on the first and remaining terms in the denom
the numerator is incredibly close to zero -- placing an if statement that checks for numerator below a certain threshold value (i.e., 1e-8) and continues the loop if achieved rather than allowing 0 / denominator to happen resolves the Error in this particular case, but not sure if that would hold for all cases or be appropriate (break is probably more sensible if that does end up working)

Something about the dot or cross product of the two intersecting lines or the fact that they are both EXACTLY at the origin to a very high level of precision may be causing some kind of issue. Here's the plot of the lines: abplot

orbeckst commented 6 years ago

Is this a duplicate of the mysterious issue #945 ?

rzu512 commented 6 years ago

Does that mean the division by 0 is caused by having a 0 derivative in the Newton-Raphson method? https://www.quora.com/When-does-Newton-Raphson-fail

MDAnalysis / mdanalysis