Use of numpy quaternion

[din14970]

I recently had a very good experience using the numpy-quaternion package for dealing with rotations/orientations. It may be an idea to use this package as a dependency in orix instead of maintaining a custom quaternion class; perhaps even a thin wrapper could be used.

[pc494]

I agree in principle, but given our current shortage of developers this might have to go in the "would be nice to have" column.

[hakonanes]

The package looks nice and easy to work with, and you bringing it up here confirms that it is good to list it among the projects users of orix might find useful in our documentation! Can you give some examples where using numpy-quaternion was better than using orix? That might motivate someone to work on this.

Here are my thoughts regarding wrapping numpy-quaternion:

Pros:

• We get access to a lot of numpy.ndarray class methods for free.
• Quaternion multiplication is about 50x faster than in orix!

Cons:

• Not really a "con", but I would say maintaining our Quaternion class is easy, since it hasn't changed much since I started contributing in 2020. It should be noted here that some of the functionality in the numpy-quaternion class we have in our Rotation class.
• Vector multiplication is about 4.3x slower than in orix.

---

orix uses slightly more memory in quaternion and vector multiplication. I tested timings with

import numpy as np
import quaternion
from orix.quaternion import Rotation
from orix.vector import Vector3d

v_size = 10
q_size = 1000
v = Vector3d(np.random.random(v_size * 3).reshape((v_size, 3)))
q_orix = Rotation.random((q_size,))  # Not available in Quaternion
q_numpy = quaternion.as_quat_array(q_orix.data)

%timeit q_orix.outer(v)
# 321 µs ± 7.3 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

%timeit quaternion.rotate_vectors(q_numpy, v.data)
# 1.39 ms ± 8.83 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

np.allclose(
    q_orix.outer(v).data,
    quaternion.rotate_vectors(q_numpy, v.data)
)
# True

%timeit q_orix * q_orix
# 149 µs ± 1.21 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

%timeit q_numpy * q_numpy
# 3.01 µs ± 35.2 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

np.allclose(
    (q_orix * q_orix).data,
    quaternion.as_float_array(q_numpy * q_numpy)
)
# True

[harripj]

It would be interesting to rewrite the Quaternion class as a wrapper around this package to see the lay of the land!

[hakonanes]

If a wrapper is attempted, I would like to see #198 be adressed at the same time :)

[din14970]

among the projects users of orix might find useful in our documentation!

I actually became aware of this project because I saw you starred in on the github feed. So thanks for that :D.

some examples where using numpy-quaternion was better

In the case of my project it was mainly because I wanted minimal and non GPL licensed dependencies, as I don't trust the company I made it for to respect open source licenses. Purely functionality wise, for most of the basic things I would say it's pretty similar, but as you say is more tightly integrated with numpy and also offers some support for calculus related operations.

It should be noted here that some of the functionality in the numpy-quaternion class we have in our Rotation class.

This I noticed and honestly now looking at it find it more confusing than helpful to have this split. In the end a quaternion is pretty much equal to a rotation.

Vector multiplication is about 4.3x slower than in orix.

This I found very curious but I did see that the documentation had a big asterisk around this function, see the bottom of: https://quaternion.readthedocs.io/en/latest/Package%20API%3A/quaternion/#rotate_vectors. Actually found a better (though ugly looking) way to do the exact same thing and get a slightly better performance than orix, I'm curious why they do it in the slow way. I should make a small PR. On my system:

%timeit quaternion.rotate_vectors(q_numpy, v.data)
# 730 µs ± 2.62 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

%timeit quaternion.as_vector_part((np.outer(q_numpy, quaternion.from_vector_part(v.data)).T*q_numpy.conj()).T)
# 173 µs ± 941 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

%timeit q_orix.outer(v)
# 237 µs ± 1.45 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

orix = q_orix.outer(v).data
qnumpy1 = quaternion.rotate_vectors(q_numpy, v.data)
qnumpy2 = quaternion.as_vector_part((np.outer(q_numpy, quaternion.from_vector_part(v.data)).T*q_numpy.conj()).T)

print(np.allclose(orix, qnumpy1))  #True
print(np.allclose(orix, qnumpy2))  #True

[din14970]

PR submitted with minor modifications to the code above moble/quaternion#191

[hakonanes]

This I noticed and honestly now looking at it find it more confusing than helpful to have this split. In the end a quaternion is pretty much equal to a rotation.

I see from the docstrings (Rotation, Quaternion) and the code that rotations can be improper, while there is no notion of proper or improper quaternions. Rotation is the only class inheriting from Quaternion, so I don't see any technical reason for the split.

PR submitted with minor modifications to the code above moble/quaternion#191

Thank you for making the PR, a 1.5x faster vector multiplication compared to orix is much better than 4.3x slower one. If a release is made with your change, I see no drawbacks to wrapping numpy-quaternion in our Quaternion class. We do lots of quaternion-vector multiplications (not so many quaternion-quaternion multiplications, actually) in kikuchipy, for example in orientation/projection center refinements, so kikuchipy will surely benefit from these speed-ups.

[hakonanes]

Might quaternionic be a better option? Developed by the same people behind numpy-quaternion, but in pure Python, no C. Would be good to do the same timings as above.

[din14970]

This guy really wants to be known as the quaternions guy it seems - he also maintains a package for quaternions in Julia.

[din14970]

Update on this: apparently the maintainer was not very happy with my PR but it did prompt him to look into speeding up the function himself. So he now made a new release where rotate_vectors is now substantially faster, taking only 90 microseconds for the example above. Win either way. So I would say the numpy-quaternion package is something to consider. Of course quaternionic may also be fine, just haven't looked at any benchmarks there.

[hakonanes]

apparently the maintainer was not very happy with my PR but it did prompt him to look into speeding up the function himself. So he now made a new release where rotate_vectors is now substantially faster, taking only 90 microseconds for the example above. Win either way.

I see that, but as you say, it might be that he hadn't improved rotate_vectors() if not for your PR. Win either way!

The speed-up makes the use of numpy-quaternion in orix even more benefitial, hopefully someone has time to look into this...

[hakonanes]

Quaternion arithmetic with numpy-quaternion is now in master, and soon released in 0.8.2.

I suggest we leave this open, because there might be other places in the code where we can integrate more tightly with numpy-quaternion.