Implement quaternion arithmetic from numpy-quaternion #281

harripj · 2022-02-14T12:37:48Z

Description of the change

As discussed in #267 and #264, I thought it best to start to implement quaternion arithmetic from numpy-quaternion into orix.quaternion.Quaternion given the performance advantage it has shown, amongst other benefits. The purpose of this PR is to lay some groundwork and discuss before making any concrete changes. To this end I have created a QuaternionNumpy class which inherits from Quaternion to help with testing.

Currently quaternion multiplication (Quaternion and Vector3d), conjugation, and outer products have been overloaded to use numpy-quaternion and give the same results. A possible extension would be to subclass orix.quaternion.Quaternion from numpy-quaternion but this would complicate things further and I don't think this is the right direction to go in. Just by using the numpy-quaternion arithmetic we get a computational speed advantage rather easily. Memory profiling to be done later.

One (potential) downside is that numpy-quaternion will be a further dependency.

@din14970 I achieved the current Quaternion * Vector3d functionality using your logic as written in moble/quaternion#191. Orix currently doesn't perform the outer product when calculating Quaternion * Vector3d so this works well, and quaternion.rotate_vectors is implemented only in Quaternion.outer.

Progress of the PR

[n/a] Docstrings for all functions
Unit tests with pytest for all lines
Clean code style by running black via pre-commit

Minimal example of the bug fix or new feature

from orix.quaternion import Quaternion, Rotation, QuaternionNumpy
import numpy as np
from orix.vector import Vector3d
import time

def test_orix_numpy_quaternion_rotate_vector_shapes():
    q_shape = [(1,), (1,), (4,), (4, 5), (4, 5), (1,), (6,), (4, 5), (4, 2, 5), (4, 2, 5)]
    v_shape = [(1,), (4,), (1,), (1,), (6,), (4, 5), (4, 5), (5,), (2, 5), (3, 6)]

    print('Q shape\t\t V shape\tOrix\t\t Numpy\t\tEqual\t\tOuter V eq\t\tOuter Q eq')

    for qs, vs in zip(q_shape, v_shape):
        q = Rotation.random(qs)
        qn = QuaternionNumpy(q.data)
        assert q.shape == qn.shape

        v = np.random.randn(*vs, 3)
        v3d = Vector3d(v)

        try:
            vo = q * v3d
            passes_orix = True
        except ValueError:
            passes_orix = False

        try:
            vn = qn * v3d
            passes_numpy = True
        except ValueError:
            passes_numpy = False

        assert passes_numpy == passes_orix
        if all([passes_numpy, passes_orix]):
            ac = np.allclose(vo.data, vn.data)
        else:
            ac = 'n/a'
        
        oov = q.outer(v3d)
        onv = qn.outer(v3d)

        q2 = Rotation.random(vs)
        ooq = q.outer(q2)
        onq = qn.outer(q2)
        print(f'{qs}\t\t{vs}\t\t{passes_orix}\t\t{passes_numpy}\t\t{ac}\t\t{np.allclose(oov.data, onv.data)}\t\t\t{np.allclose(ooq.data, onq.data)}')

test_orix_numpy_quaternion_rotate_vector_shapes()

Q shape		 V shape	Orix		 Numpy		Equal		Outer V eq		Outer Q eq
(1,)		(1,)		True		True		True		True			True
(1,)		(4,)		True		True		True		True			True
(4,)		(1,)		True		True		True		True			True
(4, 5)		(1,)		True		True		True		True			True
(4, 5)		(6,)		False		False		n/a		True			True
(1,)		(4, 5)		True		True		True		True			True
(6,)		(4, 5)		False		False		n/a		True			True
(4, 5)		(5,)		True		True		True		True			True
(4, 2, 5)	(2, 5)		True		True		True		True			True
(4, 2, 5)	(3, 6)		False		False		n/a		True			True

Some timings...

Vector operations

qs, vs = (3,), (1,)
q = Quaternion(Rotation.random(qs))
qn = QuaternionNumpy(q)

v = Vector3d(np.random.randn(*vs, 3))
%timeit q * v
%timeit qn * v

71.8 µs ± 5.31 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
37.3 µs ± 1.49 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

qs, vs = (120, 12), (1,)
q = Quaternion(Rotation.random(qs))
qn = QuaternionNumpy(q)

v = Vector3d(np.random.randn(*vs, 3))
%timeit q * v
%timeit qn * v

136 µs ± 8.5 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
66.9 µs ± 3.46 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

qs, vs = (4, 2, 5), (2, 5,)
q = Quaternion(Rotation.random(qs))
qn = QuaternionNumpy(q)

v = Vector3d(np.random.randn(*vs, 3))
%timeit q * v
%timeit qn * v

82.5 µs ± 3.68 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
41.3 µs ± 863 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

Quaternion.outer(Vector3d) seems to be on par for small calculations or quicker for larger ones (using quaternion.rotate_vectors):

qs, vs = (4, 2, 5), (2, 5,)
q = Quaternion(Rotation.random(qs))
qn = QuaternionNumpy(q)

v = Vector3d(np.random.randn(*vs, 3))
%timeit q.outer(v)
%timeit qn.outer(v)

112 µs ± 14.3 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
117 µs ± 4.97 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

qs, vs = (140, 21, 5), (21, 5,)
q = Quaternion(Rotation.random(qs))
qn = QuaternionNumpy(q)

v = Vector3d(np.random.randn(*vs, 3))
%timeit q.outer(v)
%timeit qn.outer(v)

66.1 ms ± 6.97 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
17.9 ms ± 1.78 ms per loop (mean ± std. dev. of 7 runs, 100 loops each)

Quaternion operations

qs, vs = (140, 21, 5), (21, 5,)
q1 = Quaternion(Rotation.random(qs))
q1n = QuaternionNumpy(q1)

q2 = Quaternion(Rotation.random(vs))
q2n = QuaternionNumpy(q2)
%timeit q * q2
%timeit q1n * q2n

574 µs ± 34 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
201 µs ± 11.8 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

qs, vs = (3, 2, 5), (1, 5,)
q1 = Quaternion(Rotation.random(qs))
q1n = QuaternionNumpy(q1)

q2 = Quaternion(Rotation.random(vs))
q2n = QuaternionNumpy(q2)
%timeit q * q2
%timeit q1n * q2n

725 µs ± 31 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
17.5 µs ± 1.15 µs per loop (mean ± std. dev. of 7 runs, 100000 loops each)

qs, vs = (3, 2, 5), (1, 5,)
q1 = Quaternion(Rotation.random(qs))
q1n = QuaternionNumpy(q1)

q2 = Quaternion(Rotation.random(vs))
q2n = QuaternionNumpy(q2)
%timeit q.outer(q2)
%timeit q1n.outer(q2n)

3.19 ms ± 78.5 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
23.1 µs ± 1.7 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

qs, vs = (3, 2, 5), (100, 5,)
q1 = Quaternion(Rotation.random(qs))
q1n = QuaternionNumpy(q1)

q2 = Quaternion(Rotation.random(vs))
q2n = QuaternionNumpy(q2)
%timeit q.outer(q2)
%timeit q1n.outer(q2n)

558 ms ± 18.6 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
169 µs ± 22.6 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

For reviewers

The PR title is short, concise, and will make sense 1 year later.
New functions are imported in corresponding __init__.py.
New features, API changes, and deprecations are mentioned in the
unreleased section in CHANGELOG.rst.

hakonanes · 2022-02-14T13:17:14Z

This PR is great! Smart to write it alongside the current implementation for benchmarking.

I'm already convinced we should replace our own quaternion multiplication and vector rotation with numpy-quaternion. What I'd like to see in this PR is that change, and that we support the missing shapes in the table. This is necessary for pole figure plotting.

One (potential) downside is that numpy-quaternion will be a further dependency.

Their dependencies are SciPy, NumPy and Numba, which we already depend on, so this is isn't a downside in my opinion.

harripj · 2022-02-14T13:58:30Z

and that we support the missing shapes in the table. This is necessary for pole figure plotting.

Do you mean in the testing shape results shown above or something else?

hakonanes · 2022-02-14T14:11:58Z

Do you mean in the testing shape results shown above or something else?

I think the behaviour should be defined. With MTEX

>> qs = {quaternion.rand(4, 5), quaternion.rand(6), quaternion.rand(4, 2, 5)};
>> vs = {vector3d.rand(6), vector3d.rand(4, 5), vector3d.rand(3, 6)};
>> for i=1:length(qs)
..     qs{i} * vs{i}
.. end
ans = vector3d (show methods, plot)
 size: 20 x 6
ans = vector3d (show methods, plot)
 size: 6 x 20
ans = vector3d (show methods, plot)
 size: 40 x 18

So here they flatten all instances. Will that be surprising to a user? Which output shape did the user expect? I think it is OK to do so.

harripj · 2022-02-14T14:48:26Z

The default NumPy behaviour is to flatten both arrays also:

>>> import numpy as np

>>> arr1 = np.random.randn(4, 5)
>>> arr2 = np.random.randn(6, 3)

>>> np.outer(arr1, arr2).shape
(20, 18)

However this is currently not the default behaviour in orix:

>>> q1 = Rotation.random((4, 5))
>>> q2 = Rotation.random((6, 3))

>>> q1.outer(q2).shape
(4, 5, 6, 3)

Personally I like the behaviour in orix, as for each element of q1 we have all the products with q2 (or v2) and I think this makes it easy to interpret. However I appreciate this is also different to the behaviour in NumPy. This can of course be treated with a reshape (assuming correct ordering):

>>> q1 = Rotation.random((4, 5))
>>> q2 = Rotation.random((6, 3))

>>> q12 = q1.outer(q2)

>>> q1n = quaternion.from_float_array(q1.data)
>>> q2n = quaternion.from_float_array(q2.data)

>>> q12n = np.outer(q1n, q2n).reshape(q1.shape + q2.shape)
>>> np.allclose(quaternion.as_float_array(q12n), q12.data)
True

I would favour keeping the behaviour in orix, as it means not changing the current behaviour (and so potentially breaking others' code). If the consensus is that people prefer the NumPy behaviour then we could adopt this and add how to reshape the output from the flattened outer to its docstring to revert back to the current behaviour.

hakonanes · 2022-02-14T16:12:49Z

I would favour keeping the behaviour in orix, as it means not changing the current behaviour (and so potentially breaking others' code).

I see you've given this more thought than I have, as I've only encountered the undesirable ValueError when multiplying 2D quaternion and vector instances with the thought that "this should be allowed". I cannot actually think back to the use case where I encountered this. I guess I could have used outer() in that case... I think you're right in that we should keep current behaviour.

harripj · 2022-02-14T16:54:18Z

I will leave this PR as is for a couple days for discussion if anyone would be interested to pull this PR and check the functionality (Quaternion vs QuaternionNumpy). After that I will tidy it up and finish it off.

harripj · 2022-02-16T17:21:46Z

A note on memory usage between the two implementations, here is the testing script I wrote, run it with python -m memory_profiler script.py:

Code

from memory_profiler import profile
import numpy as np

from orix.quaternion import Quaternion, QuaternionNumpy, Rotation
from orix.vector import Vector3d

S1 = (2200, 1001)
S2 = S1[-1]
S3 = (27,)

# conjugation
@profile
def test_q_conj_orix():
    q = Quaternion(Rotation.random(S1))
    qc = q.conj
    qinv = ~q


@profile
def test_q_conj_numpy():
    q = QuaternionNumpy(Rotation.random(S1))
    qc = q.conj
    qinv = ~q


# straight multiplication
@profile
def test_q_multiplication_orix():
    q1 = Quaternion(Rotation.random(S1))
    q2 = Quaternion(Rotation.random(S2))
    q = q1 * q2
    assert q.shape == q1.shape


@profile
def test_q_multiplication_numpy():
    q1 = QuaternionNumpy(Rotation.random(S1))
    q2 = QuaternionNumpy(Rotation.random(S2))
    q = q1 * q2
    assert q.shape == q1.shape


@profile
def test_v_mutliplication_orix():
    q1 = Quaternion(Rotation.random(S1))
    vdata = np.random.randn(q1.shape[-1], 3)
    v1 = Vector3d(vdata)
    v = q1 * v1
    assert v.shape == q1.shape


@profile
def test_v_multiplication_numpy():
    q1 = QuaternionNumpy(Rotation.random(S1))
    vdata = np.random.randn(q1.shape[-1], 3)
    v1 = Vector3d(vdata)
    v = q1 * v1
    assert v.shape == q1.shape


# outer product
@profile
def test_q_outer_orix():
    q1 = Quaternion(Rotation.random(S1))
    q2 = Quaternion(Rotation.random(S3))
    q = q1.outer(q2)
    assert q.shape == q1.shape + q2.shape


@profile
def test_q_outer_numpy():
    q1 = QuaternionNumpy(Rotation.random(S1))
    q2 = QuaternionNumpy(Rotation.random(S3))
    q = q1.outer(q2)
    assert q.shape == q1.shape + q2.shape


@profile
def test_v_outer_orix():
    q1 = Quaternion(Rotation.random(S1))
    vdata = np.random.randn(*S3, 3)
    v1 = Vector3d(vdata)
    v = q1.outer(v1)
    assert v.shape == q1.shape + v1.shape


@profile
def test_v_outer_numpy():
    q1 = QuaternionNumpy(Rotation.random(S1))
    vdata = np.random.randn(*S3, 3)
    v1 = Vector3d(vdata)
    v = q1.outer(v1)
    assert v.shape == q1.shape + v1.shape


if __name__ == "__main__":
    # test_q_conj_orix()
    # test_q_conj_numpy()
    test_q_multiplication_orix()
    # test_q_multiplication_numpy()
    # test_v_mutliplication_orix()
    # test_v_multiplication_numpy()
    # test_q_outer_orix()
    # test_q_outer_numpy()
    # test_v_outer_orix()
    # test_v_outer_numpy()

NB. run the tests above one by one to avoid intermediate garbage collection noise.

Results

NB. results are in the order above as shown in __main__, ie. same code with orix then numpy-quaternion, orix then numpy-quaternion etc. for direct comparison.

Line #    Mem usage    Increment  Occurrences   Line Contents
=============================================================
    12    124.8 MiB    124.8 MiB           1   @profile
    13                                         def test_q_conj_orix():
    14    479.2 MiB    354.4 MiB           1       q = Quaternion(Rotation.random(S1))
    15    579.5 MiB    100.3 MiB           1       qc = q.conj
    16    714.1 MiB    134.7 MiB           1       qinv = ~q

Line #    Mem usage    Increment  Occurrences   Line Contents
=============================================================
    19    124.6 MiB    124.6 MiB           1   @profile
    20                                         def test_q_conj_numpy():
    21    479.0 MiB    354.4 MiB           1       q = QuaternionNumpy(Rotation.random(S1))
    22    613.2 MiB    134.2 MiB           1       qc = q.conj
    23    697.5 MiB     84.3 MiB           1       qinv = ~q

Line #    Mem usage    Increment  Occurrences   Line Contents
=============================================================
    27    125.3 MiB    125.3 MiB           1   @profile
    28                                         def test_q_multiplication_orix():
    29    479.7 MiB    354.4 MiB           1       q1 = Quaternion(Rotation.random(S1))
    30    479.4 MiB     -0.3 MiB           1       q2 = Quaternion(Rotation.random(S2))
    31    613.6 MiB    134.1 MiB           1       q = q1 * q2
    32    613.6 MiB      0.0 MiB           1       assert q.shape == q1.shape

Line #    Mem usage    Increment  Occurrences   Line Contents
=============================================================
    35    124.7 MiB    124.7 MiB           1   @profile
    36                                         def test_q_multiplication_numpy():
    37    479.1 MiB    354.4 MiB           1       q1 = QuaternionNumpy(Rotation.random(S1))
    38    478.9 MiB     -0.3 MiB           1       q2 = QuaternionNumpy(Rotation.random(S2))
    39    613.1 MiB    134.2 MiB           1       q = q1 * q2
    40    613.1 MiB      0.0 MiB           1       assert q.shape == q1.shape

Line #    Mem usage    Increment  Occurrences   Line Contents
=============================================================
    43    124.7 MiB    124.7 MiB           1   @profile
    44                                         def test_v_mutliplication_orix():
    45    479.1 MiB    354.4 MiB           1       q1 = Quaternion(Rotation.random(S1))
    46    479.1 MiB      0.0 MiB           1       vdata = np.random.randn(q1.shape[-1], 3)
    47    479.1 MiB      0.0 MiB           1       v1 = Vector3d(vdata)
    48    562.7 MiB     83.5 MiB           1       v = q1 * v1
    49    562.7 MiB      0.0 MiB           1       assert v.shape == q1.shape

Line #    Mem usage    Increment  Occurrences   Line Contents
=============================================================
    52    124.7 MiB    124.7 MiB           1   @profile
    53                                         def test_v_multiplication_numpy():
    54    479.2 MiB    354.5 MiB           1       q1 = QuaternionNumpy(Rotation.random(S1))
    55    479.2 MiB      0.0 MiB           1       vdata = np.random.randn(q1.shape[-1], 3)
    56    479.2 MiB      0.0 MiB           1       v1 = Vector3d(vdata)
    57    747.3 MiB    268.1 MiB           1       v = q1 * v1
    58    747.3 MiB      0.0 MiB           1       assert v.shape == q1.shape

Line #    Mem usage    Increment  Occurrences   Line Contents
=============================================================
    62    125.0 MiB    125.0 MiB           1   @profile
    63                                         def test_q_outer_orix():
    64    479.3 MiB    354.3 MiB           1       q1 = Quaternion(Rotation.random(S1))
    65    479.3 MiB      0.0 MiB           1       q2 = Quaternion(Rotation.random(S3))
    66   2294.2 MiB   1814.9 MiB           1       q = q1.outer(q2)
    67   2294.2 MiB      0.0 MiB           1       assert q.shape == q1.shape + q2.shape

Line #    Mem usage    Increment  Occurrences   Line Contents
=============================================================
    70    124.6 MiB    124.6 MiB           1   @profile
    71                                         def test_q_outer_numpy():
    72    478.9 MiB    354.3 MiB           1       q1 = QuaternionNumpy(Rotation.random(S1))
    73    478.9 MiB      0.0 MiB           1       q2 = QuaternionNumpy(Rotation.random(S3))
    74   2360.4 MiB   1881.5 MiB           1       q = q1.outer(q2)
    75   2360.4 MiB      0.0 MiB           1       assert q.shape == q1.shape + q2.shape

Line #    Mem usage    Increment  Occurrences   Line Contents
=============================================================
    78    125.0 MiB    125.0 MiB           1   @profile
    79                                         def test_v_outer_orix():
    80    479.1 MiB    354.1 MiB           1       q1 = Quaternion(Rotation.random(S1))
    81    479.1 MiB      0.0 MiB           1       vdata = np.random.randn(*S3, 3)
    82    479.1 MiB      0.0 MiB           1       v1 = Vector3d(vdata)
    83   1856.9 MiB   1377.8 MiB           1       v = q1.outer(v1)
    84   1856.9 MiB      0.0 MiB           1       assert v.shape == q1.shape + v1.shape

Line #    Mem usage    Increment  Occurrences   Line Contents
=============================================================
    87    124.7 MiB    124.7 MiB           1   @profile
    88                                         def test_v_outer_numpy():
    89    479.2 MiB    354.5 MiB           1       q1 = QuaternionNumpy(Rotation.random(S1))
    90    479.2 MiB      0.0 MiB           1       vdata = np.random.randn(*S3, 3)
    91    479.2 MiB      0.0 MiB           1       v1 = Vector3d(vdata)
    92   1945.6 MiB   1466.4 MiB           1       v = q1.outer(v1)
    93   1945.6 MiB      0.0 MiB           1       assert v.shape == q1.shape + v1.shape

It appears that there is a small memory penalty for using numpy-quaternion which is most apparent when rotating vectors.

Further testing...

@profile
def test_v_multiplication_further():
    q1 = QuaternionNumpy(Rotation.random(S1))
    vdata = np.random.randn(q1.shape[-1], 3)
    v1 = Vector3d(vdata)
    # from source code
    q1 = quaternion.from_float_array(q1.data)
    v_intermediate = quaternion.from_vector_part(v1.data)
    mult1 = q1 * v_intermediate
    q1_inv = ~q1
    mult2 = mult1 * q1_inv
    v = quaternion.as_vector_part(mult2)

Line #    Mem usage    Increment  Occurrences   Line Contents
=============================================================
    97    124.2 MiB    124.2 MiB           1   @profile
    98                                         def test_v_multiplication_further():
    99    479.3 MiB    355.1 MiB           1       q1 = QuaternionNumpy(Rotation.random(S1))
   100    479.4 MiB      0.0 MiB           1       vdata = np.random.randn(q1.shape[-1], 3)
   101    479.4 MiB      0.0 MiB           1       v1 = Vector3d(vdata)
   102                                             # from source code
   103    546.1 MiB     66.7 MiB           1       q1 = quaternion.from_float_array(q1.data)
   104    546.1 MiB      0.0 MiB           1       v_intermediate = quaternion.from_vector_part(v1.data)
   105    546.1 MiB      0.0 MiB           1       mult1 = q1 * v_intermediate
   106    613.3 MiB     67.2 MiB           1       q1_inv = ~q1
   107    680.5 MiB     67.2 MiB           1       mult2 = mult1 * q1_inv
   108    680.5 MiB      0.0 MiB           1       v = quaternion.as_vector_part(mult2)

The q1_inv array is intermediate so I suppose would soon be garbage collected, bringing the overall memory footprint down, but peak memory usage is still increased.

Overall I think memory footprint is outweighed by the performance (CPU) gains and it would be worth including numpy-quaternion in orix.

hakonanes · 2022-02-16T17:34:26Z

Thanks for doing profiling as well! So quaternion/vector multiplication with numpy-quaternion is 2 times as fast, but uses 1.3 times the memory (as seen in profiling of test_v_multiplication_numpy()). I agree that the gain in speed outweighs the increased memory usage, at least for the use cases I can think of.

harripj · 2022-02-17T11:46:26Z

So a note on the failing test:

orix/orix/tests/test_miller.py

Lines 210 to 220 in ff6e541

    
           def test_unique(self): 
        
               # From the "Crystal geometry" notebook 
        
               diamond = Phase(space_group=227) 
        
               m = Miller.from_highest_indices(phase=diamond, uvw=[10, 10, 10]) 
        
               assert m.size == 9260 
        
               m2 = m.unique(use_symmetry=True) 
        
               assert m2.size == 450 
        
               m3, idx = m2.unit.unique(return_index=True) 
        
               assert m3.size == 345 
        
               assert isinstance(m3, Miller) 
        
               assert np.allclose(idx[:10], [74, 448, 434, 409, 374, 330, 278, 447, 412, 432])

orix/orix/tests/test_miller.py

Line 216 in ff6e541

assert m2.size == 450

This line says the test should have 450 unique vectors, however after using the numpy-quaternion implementation, this is not the case and we get 285 vectors.

On further testing this appears to be possibly related to rounding errors (NB I think line 673 should be removed from this function):

orix/orix/vector/miller.py

Lines 645 to 689 in ff6e541

    
               def unique(self, use_symmetry=False, return_index=False): 
        
                   """Unique vectors in `self`. 
        
                   Parameters 
        
                   ---------- 
        
                   use_symmetry : bool, optional 
        
                       Whether to consider equivalent vectors to compute the unique 
        
                       vectors. Default is False. 
        
                   return_index : bool, optional 
        
                       Whether to return the indices of the (flattened) data where 
        
                       the unique entries were found. Default is False. 
        
                   Returns 
        
                   ------- 
        
                   Miller 
        
                       Flattened unique vectors. 
        
                   idx : numpy.ndarray 
        
                       Indices of the unique data in the (flattened) array. 
        
                   """ 
        
                   out = super().unique(return_index=return_index) 
        
                   if return_index: 
        
                       v, idx = out 
        
                   else: 
        
                       v = out 
        
                   if use_symmetry: 
        
                       operations = self.phase.point_group 
        
                       n_v = v.size 
        
                       v2 = operations.outer(v).flatten() 
        
                       v2 = operations.outer(v).flatten().reshape(*(n_v, operations.size)) 
        
                       data = v2.data 
        
                       data_sorted = np.zeros_like(data) 
        
                       for i in range(n_v): 
        
                           a = data[i] 
        
                           order = np.lexsort(a.T)  # Sort by column 1, 2, then 3 
        
                           data_sorted[i] = a[order] 
        
                       _, idx = np.unique(data_sorted, return_index=True, axis=0) 
        
                       v = v[idx[::-1]] 
        
                   m = self.__class__(xyz=v.data, phase=self.phase) 
        
                   m.coordinate_format = self.coordinate_format 
        
                   if return_index: 
        
                       return m, idx 
        
                   else: 
        
                       return m

Using the original implementation without numpy-quaternion if the data is rounded to 14 decimal places in line 675 then we also get only 285 unique vectors. So the question is do we have a ground truth here? @hakonanes git blame informs me you wrote the test, what do you think?

hakonanes · 2022-02-17T12:10:47Z

14 decimals sounds more than enough. I'll have a look the results with rounding if they are as expected.

harripj · 2022-02-17T15:25:21Z

14 decimals sounds more than enough.

Yes, that's what I thought too. So then it follows as to whether we actually expect all 450 unique vectors, or whether some are the same within this margin.

hakonanes · 2022-02-17T15:40:09Z

Tested the return from Miller.unique(use_symmetry=True) with the following example from the user guide:

from orix.crystal_map import Phase
from orix.vector import Miller


ver = "stable"  # 0.8.1
#ver = "paddy"  # harripj:numpy-quaternion

diamond = Phase(space_group=227)

# Figure 1 (both hemispheres in all cases)
h = Miller.from_highest_indices(phase=diamond, uvw=[10, 10, 10])  # stable: 9260; paddy: 9260
fig = h.scatter(hemisphere="both", return_figure=True, figure_kwargs=dict(figsize=(20, 10)))
fig.tight_layout()
fig.savefig(os.path.join(dir_test, f"diamond_pf_{ver}.png"), **savefig_kwargs)

# Figure 2
h2 = h.unique(use_symmetry=True)  # stable: 450; paddy: 285
fig2 = h2.scatter(hemisphere="both", return_figure=True, figure_kwargs=dict(figsize=(20, 10)))
fig2.tight_layout()
fig2.savefig(os.path.join(dir_test, f"diamond_unique_pf_{ver}.png"), **savefig_kwargs)

# Figure 3
h3 = h[h <= diamond.point_group.fundamental_sector]  # stable: 285; paddy: 285
fig3 = h3.scatter(hemisphere="both", return_figure=True, figure_kwargs=dict(figsize=(20, 10)))
fig3.tight_layout()
fig3.savefig(os.path.join(dir_test, f"diamond_sector_pf_{ver}.png"), **savefig_kwargs)

stable	paddy

With this PR, the number of Miller indices obtained by considering only unique indices under symmetry, 285, is equal to the number of indices within the fundamental sector, 285 (comparing row 2 and 3). This is not the case for the current release, which returns 450 indices in the first instance, and 285 in the second. If the fundamental sector of m-3m is defined correctly (which I believe it is), the numbers should be the same. Thus, I conclude that there is an error in the current release which this PR solves.

What do you think, @harripj?

hakonanes · 2022-02-17T15:45:56Z

By changing this line in Miller.unique()

orix/orix/vector/miller.py

Line 675 in ff6e541

data = v2.data

to data = v2.data.round(16), we get the same 450 indices unique under symmetry from our instance of 9260 indices; using 15 we get 449, and using 14 we get 285. Changing the line to data = v2.data.astype(int) (Miller indices aren't restricted to integers, so don't do this!) also returns 285 in this case.

You've fixed a bug here, and I'm very grateful :)

hakonanes · 2022-02-17T15:48:43Z

(NB I think line 673 should be removed from this function):

Yes please!

harripj · 2022-02-17T15:52:24Z

This makes sense to me, in fact it would be great if unique chose the same fundamental sector on the unit sphere as point_group.fundamental_sector.

to data = v2.data.round(16), we get the same 450 indices unique under symmetry from our instance of 9260 indices; using 15 we get 449, and using 14 we get 285. Changing the line to data = v2.data.astype(int) (Miller indices aren't restricted to integers, so don't do this!) also returns 285 in this case.

This test is very informative! So it looks like the correct test value should in fact be 285. On how we got there though, I am not sure that simply implementing numpy-quaternion properly fixes the bug. I will open a new issue to discuss.

hakonanes · 2022-02-17T16:01:00Z

it would be great if unique chose the same fundamental sector on the unit sphere as point_group.fundamental_sector.

We could do that within unique() when use_symmetry=True. Or, we could fix Miller.in_fundamental_sector() so that it works (doesn't at the moment, my bad) and projects the indices, in this case closest to the [100] pole, into the sector.

On how we got there though, I am not sure that simply implementing numpy-quaternion properly fixes the bug. I will open a new issue to discuss.

I'm OK with just changing the test and adding a note to the changelog. orix supports all 32 crystallographic point groups in most cases, only the 11 Laue classes in others. It's on me that the test was wrong in the first place, but there are just so many tests one can write. What's most important to me is that we fix bugs and release patches quickly as we spot them (hint hint #243). This PR changes nothing in the API, I suggest releasing a patch release 0.8.2 after this is merged.

hakonanes · 2022-02-17T16:02:21Z

I will open a new issue to discuss.

I didn't mean to say that you shouldn't open the issue, btw. Sorry if that's what it read like.

harripj · 2022-02-17T16:05:14Z

I didn't mean to say that you shouldn't open the issue, btw. Sorry if that's what it read like.

Don't worry, I didn't take it that way at all, in fact I was writing #285 as the notifications were coming in! Best to redirect this conversation there as these are now two different issues.

hakonanes · 2022-02-18T09:21:56Z

it would be great if unique chose the same fundamental sector on the unit sphere as point_group.fundamental_sector

Actually, I find the assumption that the user wants the vectors within the fundamental sector not strong enough to be used in Miller.unique(). We should instead fix Miller.in_fundamental_sector() (inherited from Vector3d) so that it works, so that the unique ones within the sector can be found by projection followed by rejection.

If you want the direct or reciprocal lattice vectors within the fundamental sector in one call, we should make a convenience method Miller.from_fundamental_sector() calling these different methods. Will that meet your needs @harripj?

orix/quaternion/quaternion.py

harripj · 2022-02-18T10:41:18Z

@hakonanes I have just added the equations for quaternion-quaternion and quaternion-vector multiplication to the Quaternion doctsring. Let me know if this is what you had in mind?

harripj · 2022-02-18T11:13:39Z

Update: this is waiting on #285, after this the tests should pass.

pc494

I'm happy with this, will let @hakonanes merge.

Add QuaternionNumpy class which uses numpy-quaternion

ea00cc6

harripj added discussion dev Package maintenance help-wanted A little help with this would be nice labels Feb 14, 2022

add numpy-quaternion to setup.py

072ee83

hakonanes added this to the v0.9.0 milestone Feb 14, 2022

hakonanes modified the milestones: v0.9.0, v1.0.0 Feb 14, 2022

harripj added 6 commits February 17, 2022 11:15

rm QNumpy class, add np-quat funcs to Quaternion

da700eb

remove QuaternionNumpy import

882cc16

Merge branch 'master' into numpy-quaternion

eadad58

Merge remote-tracking branch 'upstream/master' into numpy-quaternion

67f48fc

update changelog

e332cfc

remove QuaternionNumpy from __init__

855bdeb

harripj mentioned this pull request Feb 17, 2022

Implement rounding consistently to 12 decimal places in orix, eg. in unique() functions #285

Merged

4 tasks

hakonanes mentioned this pull request Feb 17, 2022

0.8.2 patch release #286

Closed

7 tasks

hakonanes reviewed Feb 18, 2022

View reviewed changes

orix/quaternion/quaternion.py Show resolved Hide resolved

hakonanes reviewed Feb 18, 2022

View reviewed changes

orix/quaternion/quaternion.py Show resolved Hide resolved

add explicit q*q and q*v equations to docstring

d8eba91

harripj removed status : wip help-wanted A little help with this would be nice labels Feb 18, 2022

remove QuaternionNumpy class

4be8c88

harripj requested a review from pc494 February 18, 2022 11:52

Merge branch 'master' into numpy-quaternion

8f5bcb7

pc494 approved these changes Feb 18, 2022

View reviewed changes

hakonanes merged commit 25b69ba into pyxem:master Feb 18, 2022

harripj deleted the numpy-quaternion branch February 18, 2022 19:17

harripj mentioned this pull request Apr 21, 2022

Add Dask-based lazy computation of quaternion-vector outer product #321

Merged

6 tasks

harripj mentioned this pull request May 27, 2022

Using scipy's Rotation object #346

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Implement quaternion arithmetic from numpy-quaternion #281

Implement quaternion arithmetic from numpy-quaternion #281

harripj commented Feb 14, 2022 •

edited by hakonanes

Loading

hakonanes commented Feb 14, 2022

harripj commented Feb 14, 2022

hakonanes commented Feb 14, 2022

harripj commented Feb 14, 2022

hakonanes commented Feb 14, 2022

harripj commented Feb 14, 2022

harripj commented Feb 16, 2022

hakonanes commented Feb 16, 2022

harripj commented Feb 17, 2022

hakonanes commented Feb 17, 2022

harripj commented Feb 17, 2022

hakonanes commented Feb 17, 2022

hakonanes commented Feb 17, 2022 •

edited

Loading

hakonanes commented Feb 17, 2022

harripj commented Feb 17, 2022

hakonanes commented Feb 17, 2022

hakonanes commented Feb 17, 2022

harripj commented Feb 17, 2022 •

edited

Loading

hakonanes commented Feb 18, 2022

harripj commented Feb 18, 2022

harripj commented Feb 18, 2022

pc494 left a comment

Implement quaternion arithmetic from numpy-quaternion #281

Implement quaternion arithmetic from numpy-quaternion #281

Conversation

harripj commented Feb 14, 2022 • edited by hakonanes Loading

Description of the change

Progress of the PR

Minimal example of the bug fix or new feature

Vector operations

Quaternion operations

For reviewers

hakonanes commented Feb 14, 2022

harripj commented Feb 14, 2022

hakonanes commented Feb 14, 2022

harripj commented Feb 14, 2022

hakonanes commented Feb 14, 2022

harripj commented Feb 14, 2022

harripj commented Feb 16, 2022

Code

Results

hakonanes commented Feb 16, 2022

harripj commented Feb 17, 2022

hakonanes commented Feb 17, 2022

harripj commented Feb 17, 2022

hakonanes commented Feb 17, 2022

hakonanes commented Feb 17, 2022 • edited Loading

hakonanes commented Feb 17, 2022

harripj commented Feb 17, 2022

hakonanes commented Feb 17, 2022

hakonanes commented Feb 17, 2022

harripj commented Feb 17, 2022 • edited Loading

hakonanes commented Feb 18, 2022

harripj commented Feb 18, 2022

harripj commented Feb 18, 2022

pc494 left a comment

Choose a reason for hiding this comment

harripj commented Feb 14, 2022 •

edited by hakonanes

Loading

hakonanes commented Feb 17, 2022 •

edited

Loading

harripj commented Feb 17, 2022 •

edited

Loading