rotation_matrix issues

[MLopez-Ibanez]

🐞 Problem

Perhaps I'm missing some detail but the implementation of rotation_matrix looks strange and it doesn't match the one in astropy (but the one of astropy also looks strange).

According to Wikipedia: https://en.wikipedia.org/wiki/Rotation_matrix#Basic_rotations
the rotation for axis=1 should be

[[c, 0.0, s], [0.0, 1.0, 0.0], [-s, 0.0, c]]

instead of

poliastro/src/poliastro/core/util.py

Line 33 in 18a390d

return np.array([[c, 0.0, s], [0.0, 1.0, 0.0], [s, 0.0, c]])

Moreover, the version in astropy:

https://github.com/astropy/astropy/blob/39cb148a4bd69368f6f30b9abfe32112a42f58e6/astropy/coordinates/matrix_utilities.py#L92-L99

        a1 = (i + 1) % 3
        a2 = (i + 2) % 3
        R = np.zeros(getattr(angle, 'shape', ()) + (3, 3))
        R[..., i, i] = 1.
        R[..., a1, a1] = c
        R[..., a1, a2] = s
        R[..., a2, a1] = -s
        R[..., a2, a2] = c

Actually creates:

which probably is due to the fact that "the rotation sense is counterclockwise looking down the + axis (e.g. positive rotations obey left-hand-rule)." But I don't understand the matrix generated by poliastro.

[astrojuanlu]

The Astropy result could an alias vs alibi kind of confusion https://en.m.wikipedia.org/wiki/Rotation_matrix#Ambiguities

For our code, this test is not very convincing:

poliastro/tests/tests_core/test_core_util.py

Line 18 in 18a390d

def test_rotation_matrix_y():

Maybe this is a genuine bug? That's embarrassing...

[MLopez-Ibanez]

The Astropy result could an alias vs alibi kind of confusion https://en.m.wikipedia.org/wiki/Rotation_matrix#Ambiguities

Surprisingly, you do use the astropy function in quite a number of places in poliastro, hopefully, aware of the differences...

By the way, after a few experiments, I do believe that the astropy function can be both vectorized and jitted with a few modifications and the latest numba version.

Maybe this is a genuine bug? That's embarrassing...

I don't see axis=1 ever used in the code, so it wouldn't be surprising that nobody noticed until now. Is there some testsuite that could be used to validate the functions?

[MLopez-Ibanez]

@astrojuanlu

I'm happy to send a pull request but I would like to make sure what the correct behavior should be.

[MLopez-Ibanez]

Just for the record in case I forget about this. I created a vectorized by angle and JITable version of rotation_matrix, however, it does not match the code currently implemented in poliastro nor astropy. Whatever is the correct code, it should be easy to adjust for it:

@jit
def rotation_matrix_vec(angle, axis):
    assert axis in (0,1,2)
    angle = np.asarray(angle)
    c = cos(angle)
    s = sin(angle)
    i = axis
    a1 = (i + 1) % 3
    a2 = (i + 2) % 3
    R = np.zeros(angle.shape + (3, 3))
    R[..., i, i] = 1.
    R[..., a1, a1] = c
    R[..., a1, a2] = -s
    R[..., a2, a1] = s
    R[..., a2, a2] = c
    return R

Moreover, I believe the above function is nice to have as a utility, however, given that there are only 3 possible values of axis, it would be better to create 3 specializations that are simpler and faster:

import numpy as np
from numba import njit as jit
from numpy import cos, sin

# Original from poliastro
@jit
def rotation_matrix(angle, axis):
    c = cos(angle)
    s = sin(angle)
    if axis == 0:
        return np.array([[1.0, 0.0, 0.0], [0.0, c, -s], [0.0, s, c]])
    elif axis == 1:
        # Original: np.array([[c, 0.0, s], [0.0, 1.0, 0.0], [s, 0.0, c]])
        # Fixed ???
        return np.array([[c, 0.0, s], [0.0, 1.0, 0.0], [-s, 0.0, c]])
    elif axis == 2:
        return np.array([[c, -s, 0.0], [s, c, 0.0], [0.0, 0.0, 1.0]])
    else:
        raise ValueError("Invalid axis: must be one of 'x', 'y' or 'z'")

@jit
def rotation_matrix_vec(angle, axis):
    assert axis in (0,1,2)
    angle = np.asarray(angle)
    c = cos(angle)
    s = sin(angle)

    i = axis
    a1 = (i + 1) % 3
    a2 = (i + 2) % 3
    R = np.zeros(angle.shape + (3, 3))
    R[..., i, i] = 1.
    R[..., a1, a1] = c
    R[..., a1, a2] = -s
    R[..., a2, a1] = s
    R[..., a2, a2] = c
    return R

@jit
def rotation_matrix_x(angle):
    angle = np.asarray(angle)
    c = cos(angle)
    s = sin(angle)
    R = np.zeros(angle.shape + (3, 3))
    R[..., 0, 0] = 1.
    R[..., 1, 1] = c
    R[..., 1, 2] = -s
    R[..., 2, 1] = s
    R[..., 2, 2] = c
    return R

@jit
def rotation_matrix_y(angle):
    angle = np.asarray(angle)
    c = cos(angle)
    s = sin(angle)
    a1 = 2
    a2 = 0
    R = np.zeros(angle.shape + (3, 3))
    R[..., 1, 1] = 1.
    R[..., 2, 2] = c
    R[..., 2, 0] = -s
    R[..., 0, 2] = s
    R[..., 0, 0] = c
    return R

@jit
def rotation_matrix_z(angle):
    angle = np.asarray(angle)
    c = cos(angle)
    s = sin(angle)

    a1 = 0
    a2 = 1
    R = np.zeros(angle.shape + (3, 3))
    R[..., 2, 2] = 1.
    R[..., 0, 0] = c
    R[..., 0, 1] = -s
    R[..., 1, 0] = s
    R[..., 1, 1] = c
    return R


def _test_rotation_matrix(angle, axis):
    angle = np.asarray(angle)
    expected = rotation_matrix(angle, axis)
    result = rotation_matrix_vec(angle, axis)
    assert np.allclose(expected, result),f'exp={expected}\nres={result}'
    if axis==0:
        expected = rotation_matrix_x(angle)
    elif axis==1:
        expected = rotation_matrix_y(angle)
    elif axis==2:
        expected = rotation_matrix_z(angle)
    assert np.allclose(expected, result), f'exp={expected}\nres={result}'
    
    
def test_rotation_matrix_x():
    _test_rotation_matrix(0.218,0)

def test_rotation_matrix_y():
    _test_rotation_matrix(0.218, 1)

def test_rotation_matrix_z():
    _test_rotation_matrix(0.218, 2)

test_rotation_matrix_x()
test_rotation_matrix_y()
test_rotation_matrix_z()

[astrojuanlu]

Progress:

• As I suspected, Astropy uses the alias or passive transformation, whereas poliastro uses the alibi or active transformation:

In [33]: u = np.array([1, 0, 0])

In [34]: rotation_matrix_astropy(45, "z") @ u  # Alias/passive, coordinate system is rotated
Out[34]: array([ 0.70710678, -0.70710678,  0.        ])

In [35]: rotation_matrix_poliastro(np.radians(45), 2) @ u  # Alibi/active, point is rotated
Out[35]: array([0.70710678, 0.70710678, 0.        ])

This is more clearly explained in the docstring of the function https://github.com/astropy/astropy/blob/0eb1bf2ab367c3203ba7d98e385da906880b8d3d/astropy/coordinates/matrix_utilities.py#L51-L53

the rotation sense is counterclockwise looking down the + axis (e.g. positive rotations obey left-hand-rule).

I'm used to the right-hand rule and pre-multiplication of column vectors. We should clarify this in the docs.

• Rotations around X and Z axis give the same results:

In [41]: v = np.random.randn(3); v
Out[41]: array([-0.30387748, -1.4202498 ,  0.24305655])

In [42]: rotation_matrix_astropy(np.degrees(-0.5), "x") @ v
Out[42]: array([-0.30387748, -1.36291397, -0.46760183])

In [43]: rotation_matrix_poliastro(0.5, 0) @ v
Out[43]: array([-0.30387748, -1.36291397, -0.46760183])

In [46]: rotation_matrix_astropy(np.degrees(-0.5), "z") @ v
Out[46]: array([ 0.41422644, -1.39207308,  0.24305655])

In [47]: rotation_matrix_poliastro(0.5, 2) @ v
Out[47]: array([ 0.41422644, -1.39207308,  0.24305655])

But around Y, they give different results:

In [44]: rotation_matrix_astropy(np.degrees(-0.5), "y") @ v
Out[44]: array([-0.15015006, -1.4202498 ,  0.35898882])

In [45]: rotation_matrix_poliastro(0.5, 1) @ v
Out[45]: array([-0.15015006, -1.4202498 ,  0.06761556])

and if I change the sign as you suggested initially, everything works:

In [5]: rotation_matrix_astropy(np.degrees(-0.5), "y") @ v
Out[5]: array([-0.15015006, -1.4202498 ,  0.35898881])

In [6]: rotation_matrix_poliastro(0.5, 1) @ v
Out[6]: array([-0.15015006, -1.4202498 ,  0.35898881])

And I absolutely trust Astropy to be correct on this one, plus it's very clear that we have a typo here and that our tests are dumb.

So:

• Happy to have PRs to correct the typo and the tests
• Bonus points if we clarify that our rotation_matrix assumes right-hand rule (positive rotations are counterclockwise), pre-multiplication, and alibi/active rotation (point is rotated while coordinate system is left fixed).
• I will look for uses of the Astropy rotation_matrix to double-check that we're using it correctly
• About your proposal @MLopez-Ibanez , is it really much faster to have all the specialized functions rather than to have only one?

[astrojuanlu]

We are using the Astropy version in 2 cases:

poliastro/src/poliastro/frames/ecliptic.py

Lines 93 to 107 in 0e14e52

	def _make_rotation_matrix_from_reprs(start_representation, end_representation):
	"""
	Return the matrix for the direct rotation from one representation to a second representation.
	The representations need not be normalized first.
	"""
	A = start_representation.to_cartesian()
	B = end_representation.to_cartesian()
	rotation_axis = A.cross(B)
	rotation_angle = -np.arccos(
	A.dot(B) / (A.norm() * B.norm())
	) # Negation is required
	# This line works around some input/output quirks of Astropy's rotation_matrix()
	matrix = np.array(rotation_matrix(rotation_angle, rotation_axis.xyz.value.tolist()))
	return matrix

This one is okay, because the rotation axis is not (x, y, z), but a custom one. We don't support that (yet?).

poliastro/src/poliastro/frames/fixed.py

Line 92 in 95cf3ff

r = r.transform(rotation_matrix(-W, "z"))

(both to_equatorial and from_equatorial)

This is probably because we copy-pasted part of this code from Astropy itself. The signs look correct.

[MLopez-Ibanez]

* About your proposal @MLopez-Ibanez , is it really much faster to have all the specialized functions rather than to have only one?

It has to be faster since I far as I can tell, numba cannot do specialization by itself (that would require interprocedural optimization). But we have to define "really much faster". I will try to test it with my code for GTOC11, but that code is too massive to become a test (or even to share it).