Place.rotation_axis_and_angle not robust

[Tristan Croll]

The following bug report has been submitted:
Platform:        Linux-3.10.0-1062.9.1.el7.x86_64-x86_64-with-centos-7.7.1908-Core
ChimeraX Version: 0.92 (2020-01-11)
Description
Place.rotation_axis_and_angle() returns incorrect values when the rotation matrix is symmetric or nearly so. Following the link to http://en.wikipedia.org/wiki/Rotation_matrix from core.geometry.matrix.R_to_axis_angle():

"This does not work if R is symmetric. Above, if R − R T  is zero, then all subsequent steps are invalid. In this case, it is necessary to diagonalize and find the eigenvector corresponding to an eigenvalue of 1. "

Simple case of a flip about x:
{{{
p = Place(axes=[[1,0,0],[0,-1,0],[0,0,-1]])

p.rotation_axis_and_angle()
(array([0., 0., 1.]), 0) # should be [1,0,0],180

# A more "real" case with some numeric noise

p = Place(axes=[[1,0,0],[-1e-16,-1,0],[-1e-16,0,-1]])

p.rotation_axis_and_angle()
(array([ 0.        , -0.70710678,  0.70710678]), 180.0)

}}}

Log:
UCSF ChimeraX version: 0.92 (2020-01-11)  
© 2016-2019 Regents of the University of California. All rights reserved.  
How to cite UCSF ChimeraX  

> open 6ig2 structureFactors true

Summary of feedback from opening 6ig2 fetched from pdb  
---  
warning | WARNING: multiple experimental reflection datasets found:  
F_meas_au, F_meas_sigma_au,  
intensity_meas, intensity_sigma  
Automatically choosing "intensity_meas, intensity_sigma".  
notes | Fetching compressed mmCIF 6ig2 from
http://files.rcsb.org/download/6ig2.cif  
Fetching CCD CTP from http://ligand-expo.rcsb.org/reports/C/CTP/CTP.cif  
Fetching compressed 6ig2 structure factors from
http://files.rcsb.org/download/6ig2-sf.cif  
Resolution: 2.8817619976202815  
Reflection data provided as intensities. Performing French & Wilson scaling to
convert to amplitudes...  
  
6ig2 title:  
Structure of mitochondrial CDP-DAG synthase Tam41 complexed with CTP, δ 74,
F240A [more info...]  
  
Chain information for 6ig2  
---  
Chain | Description  
1.1/A 1.1/B 1.1/C 1.1/D | Phosphatidate cytidylyltransferase, mitochondrial  
  
Non-standard residues in 6ig2 #1.1  
---  
CTP — cytidine-5'-triphosphate  
  
6ig2 mmCIF Assemblies  
---  
1| author_and_software_defined_assembly  
2| author_and_software_defined_assembly  
  
  

> toolshed show Shell

/opt/UCSF/ChimeraX-daily/lib/python3.7/site-
packages/IPython/core/history.py:226: UserWarning: IPython History requires
SQLite, your history will not be saved  
warn("IPython History requires SQLite, your history will not be saved")  




OpenGL version: 3.3.0 NVIDIA 440.33.01
OpenGL renderer: TITAN Xp/PCIe/SSE2
OpenGL vendor: NVIDIA Corporation

[Tristan Croll]

This seems to work robustly - albeit about 4 times slower, but I think 
not so much that it would seriously impact anything (106 vs 25 us on my 
machine).

{{{
def rotation_axis_and_angle_3x3(rot33):
     import numpy
     eigvals, eigvecs = numpy.linalg.eig(rot33)
     unit_eigval = numpy.argwhere(eigvals==1)
     if not len(unit_eigval):
         return (numpy.array([0.,0.,1.]), 0)
     axis = eigvecs[:, unit_eigval[0]].T
     from math import degrees, acos
     angle = degrees(acos((numpy.trace(rot33)-1)/2))
     return axis, angle
}}}


On 2020-01-17 10:45, ChimeraX wrote:

[Tristan Croll]

Ahh... I see this is a more complex task than I thought. The eig() 
function sometimes returns real arrays, sometimes complex... but more 
importantly this approach doesn't determine the sign of the rotation. 
Also turns out that SciPy has it covered:

{{{
def rot33_as_axis_and_angle(rot33):
     from scipy.spatial.transform import Rotation as R
     r = R.from_dcm(rot33)
     from math import degrees
     rvec = r.as_rotvec()
     rnorm = numpy.linalg.norm(rvec)
     angle = degrees(rnorm)
     return (rvec/rnorm, angle)
}}}

... albeit a bit slowly (270 us). Plenty fast enough for my purpose, but 
might be a bit on the slow side if this is used anywhere in big loops? 
For what it's worth, `numpy.allclose(rot33, rot33.T)` costs about 65 us 
- so where it really matters you could use that to decide whether to use 
the fast or the robust method.


On 2020-01-17 11:06, ChimeraX wrote:

[Tom Goddard]

Fixed.

I put in special case code to handle the 180 degree rotation case (R = R_transpose).

[Tom Goddard]

The problem is for any 180 degree rotation (that is the R = Rt condition), and any round-off error leads to a pretty random axis, although the returned angle is 180. The 180 degree rotations are easy to extract the axis very fast from just the diagonal matrix elements

[[2ax2 - 1, 2ax*ay, 2ax*az],
[2*ax*ay, 2ay2-1, 2ay*az],
[2ax*az, 2ay*az, 2az2-1]]

but it could be a bit tricky to make the switch over from the non-180 to the 180 case not glitchy.

[Tristan Croll]

I guess the best approach would be to err on the side of caution, and make the switch when you still have at least a few significant figures to play with. Just messing around, a rotation of 180.0001 degrees about a random axis has differences between equivalent off-axis positions on the order of 1e-6, and the existing code still nails it at that point. Given that's about the limit of float32 precision anyway, why not put the cutoff somewhere around there?

[Tom Goddard]

Double precision (float64) is being used for the 3x4 transformation matrices.  I made the cutoff 1e-14 (rotation matrix elements are less than 1).  I did not study it carefully.  If you run experiments with the current code that suggests a bigger cutoff is better let me know the results.

[Tristan Croll]

OK. A quick run with the following snippet:

def test_rotation_accuracy():
    import numpy
    from chimerax.core.geometry import rotation
    from collections import defaultdict
    axis_errors = defaultdict(lambda: list())
    angle_errors = defaultdict(lambda: list())
    for i in range(20):
        angle = 180-10**-i
        for j in range(50):
            axis = numpy.random.rand(3)-0.5
            axis/=numpy.linalg.norm(axis)
            r = rotation(axis, angle)
            calc_axis, calc_angle = r.rotation_axis_and_angle()
            if sum(calc_axis - axis) > 0.1:
                calc_axis = -calc_axis
                calc_angle = 360-calc_angle
            axis_errors[i].append(numpy.max(numpy.abs(calc_axis-axis)))
            angle_errors[i].append(abs(calc_angle-angle))
    return dict(axis_errors), dict(angle_errors)

... says that the absolute error along any given axis dimension passes 1e-6 at about i=9 (that is, where the angle is 180-(1e-9) degrees), at which point the off-axis elements differ by about 1e-10. So that seems like a reasonable point. The error hits around 1e-3 for i=12, and by i=13 the result is pretty much random.

[Tom Goddard]

Ok I changed the cutoff to 1e-9 after some testing. My previous code for near 180 degrees was simply wrong, had wrong signs on the axis components (always made them positive). And your test was also a bit wrong, the sum(calc_axis - axis) needed to be sum(abs(calc_axis - axis)).

Axis and angle accuracies near 1e-8 aren't to great given float64 is in use, but it will do.

[Tristan Croll]

Still one more bug, I'm afraid:

p = Place()
p.rotation_axis_and_angle()

Out[22]: (array([1., 1., 1.]), 180)



On 2020-01-28 20:36, ChimeraX wrote:

[Tom Goddard]

Yikes.

[Tom Goddard]

Fixed.

The identity rotation is symmetric and was incorrectly being handled by the 180 degree rotation code which was just looking for symmetric matrices.