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

Open
opened this Issue Dec 19, 2017 · 4 comments

Projects
None yet
4 participants

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.

''')
0.16.2```

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

Member

### kain88-de commented Dec 19, 2017

 This is still present in develop.
Member

### tylerjereddy commented Dec 21, 2017

 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:
Member

### orbeckst commented Dec 22, 2017

 Is this a duplicate of the mysterious issue #945 ?

### rzu512 commented Dec 29, 2017

 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