ENH: Implemented Rotation.from_davenport and Rotation.as_davenport

[evbernardes]

Reference issue

This is a follow-up to #17392.

What does this implement/fix?

Davenport angles are a generalization of Euler angles.
Implemented the function Rotation.as_davenport, and for convenience, an inverse function called Rotation.from_davenport that makes internal calls to Rotation.from_rotvec but using the same type of parameters as Rotation.as_davenport.

Additional information

Not any generalization is possible, basically we must have dot(n1, n2) = dot(n2, n3) = 0, like in Euler angles.
Moreover, for Euler angles: dot(n1, n3) = 0, 1 or -1. For Davenport angles, this last can be any value between -1 and 1 instead (Source).

This is a follow-up to my implementation of a new algorithm for Rotation.as_euler() based on this paper (open access). I now implemented a slightly modified version of the same algorithm to compute the Davenport angles for a set of input axes in the function.

In this implementation, there may be some set differences between the angles computed for the same axes by Rotation.as_davenport and Rotation.as_euler for some asymmetrical cases, but these different angles represent the same final rotation.

[evbernardes]

@nmayorov As we discussed last time, I made this PR for the generalization of Euler angles.
Also posted on the mail list, to see if there is any interest.

[evbernardes]

Alright, finally managed to fix all of the Sphinx problems.

If there is any interest for this algorithm, I think it is ready.

[scottshambaugh]

I don't know enough to evaluate the correctness of the algorithm, but I think it's certainly very interesting and I can see uses for it! My main suggestion was going to be to have tests that show equivalence between the euler angle functions and the davenport functions, but you've already got it in there! The only other thing I would say to add is a simple (non-Euler) test case where the target rotation matrix is hard-coded and something that could be verified by hand, rather than relying on programmatic test cases for everything.

[evbernardes]

Hey there @scottshambaugh, thanks for your input! I think I might do it indeed when I have some time!

I ended up not implementing any hard-coded case because I tried to keep it not too different from the tests for as_euler, and the only hard-coded tests there are from the degenerate cases. Maybe I'll try to think of how to integrate one there…

[evbernardes]

@nmayorov in case you have some time to review it, I implemented some time ago in this PR the generalization of my conversion algorithm to Davenport angles, in case you think this could be useful in Scipy.

[nmayorov]

@evbernardes Yes, sorry I somehow forgot about this activity, I'll check it out.

[nmayorov]

@evbernardes My thoughts on quick scan:

1. I think it's fine to include it into scipy, albeit borderline. Can view scipy as a universal reference source for attitude algorithms.
2. But it's beneficial to reduce the amount of code if possible. Here I see that the algorithms are almost identical to as_euler, from_euler. Definitely we should try to factor out the common logic (like as_euler will call _compute_davenport_from_quat)
3. It's better not to introduce minor unrelated changes, they have created a bit too much noise here.

So I'd like to see a better generalization with Euler angle algorithms to reduce conceptual and technical duplication.

[evbernardes]

Hello @nmayorov, thanks for your review!

Using only the Davenport implementation is indeed a possibility I thought about. The only problem with that is that _compute_euler_from_quat runs considerably faster than _compute_davenport_from_quat, since it only does a permutation instead of multiple dot products for each rotation.

Do you think it is worth it anyways?

PS.: For the minor unrelated change, I think this is a fault of my automatic linter, maybe a separate PR should be done in order to remove these spaces altogether...

[nmayorov]

Using only the Davenport implementation is indeed a possibility I thought about. The only problem with that is that _compute_euler_from_quat runs considerably faster than _compute_davenport_from_quat, since it only does a permutation instead of multiple dot products for each rotation.

I see, I don't know how to judge it. Can you run %timeit in jupyter to compare as_euler vs as_davenport as it now?

[nmayorov]

My test has shown, that as_davenport as it currently implemented is 10 times slower.

As there were efforts to speed up everything by using Cython (which reduced readability quite a bit) we probably should not pessimise as_euler because it is supposedly used often.

Is there a way to reduce as_davenport problem to as_euler problem doing some preliminary rotation transform?

[nmayorov]

I think that will be all from me.

Fix refguide checker complaints:

=======================
scipy.spatial.transform
=======================

scipy.spatial.transform.Rotation
--------------------------------

File "build-install/lib/python3.9/site-packages/scipy/spatial/transform/_rotation.cpython-39-x86_64-linux-gnu.so", line ?, in Rotation.as_euler
Failed example:
    r = R.from_rotvec([def as_davenport
    [0, 0, np.pi/2],
    [0, -np.pi/3, 0],
    [np.pi/4, 0, 0]])
Exception raised:
    Traceback (most recent call last):
      File "/home/circleci/.pyenv/versions/3.9.18/lib/python3.9/doctest.py", line 1334, in __run
        exec(compile(example.source, filename, "single",
      File "<doctest Rotation.as_euler[8]>", line 1
        r = R.from_rotvec([def as_davenport
                           ^
    SyntaxError: invalid syntax


File "build-install/lib/python3.9/site-packages/scipy/spatial/transform/_rotation.cpython-39-x86_64-linux-gnu.so", line ?, in Rotation.as_euler
Failed example:
    r.as_euler('zxy', degrees=True)
Expected:
    array([[ 90.,   0.,   0.],
           [  0.,   0., -60.],
           [  0.,  45.,   0.]])
Got:
    array([[90.,  0.,  0.]])


File "build-install/lib/python3.9/site-packages/scipy/spatial/transform/_rotation.cpython-39-x86_64-linux-gnu.so", line ?, in Rotation.as_euler
Failed example:
    r.as_euler('zxy', degrees=True).shape
Expected:
    (3, 3)
Got:
    (1, 3)


File "build-install/lib/python3.9/site-packages/scipy/spatial/transform/_rotation.cpython-39-x86_64-linux-gnu.so", line ?, in Rotation.from_davenport
Failed example:
    r.as_quat().shape
Expected:
    [ 0.701057,  0.430459, -0.092296,  0.560986]
Got:
    (4,)


ERROR: refguide or doctests have errors

Exited with code exit status 1

[evbernardes]

Fix refguide checker complaints:

Oops! I think it's all good now.

[nmayorov]

Some pending approve required.

[evbernardes]

Some pending approve required.

Is there something I should do, or just wait until someone approves it?

[nmayorov]

Some pending approve required.

Is there something I should do, or just wait until someone approves it?

I wish I knew 😄 I see that I can merge it now. I've send one approving review.

I've left the last comment and it in the excellent state in my opinion. So I plan to merge it today or tomorrow.

[evbernardes]

Some pending approve required.

Is there something I should do, or just wait until someone approves it?

I wish I knew 😄 I see that I can merge it now. I've send one approving review.

I've left the last comment and it in the excellent state in my opinion. So I plan to merge it today or tomorrow.

Alright, awesome! Once again, thank you for your time :)

[nmayorov]

@evbernardes Great working with you. Come back with new ideas!

[Carreau]

Would you mind looking at #19512, I'm doing to modification to where/when the Gimbal Lock warning is emitted and as this touches code that was added by this recently thought that might be of interest.