From ea00cc6ebd498685447a0122f12da571802fddcb Mon Sep 17 00:00:00 2001 From: Paddy Harrison Date: Mon, 14 Feb 2022 12:53:22 +0100 Subject: [PATCH 1/8] Add QuaternionNumpy class which uses numpy-quaternion --- orix/quaternion/__init__.py | 3 +- orix/quaternion/quaternion.py | 81 ++++++++++++++++++++++++ orix/tests/quaternion/test_quaternion.py | 2 +- 3 files changed, 84 insertions(+), 2 deletions(-) diff --git a/orix/quaternion/__init__.py b/orix/quaternion/__init__.py index e5d0e9c7..d856ef69 100644 --- a/orix/quaternion/__init__.py +++ b/orix/quaternion/__init__.py @@ -28,7 +28,7 @@ orientations. """ -from orix.quaternion.quaternion import check_quaternion, Quaternion +from orix.quaternion.quaternion import check_quaternion, Quaternion, QuaternionNumpy from orix.quaternion.rotation import Rotation, von_mises from orix.quaternion.orientation import Misorientation, Orientation from orix.quaternion.orientation_region import get_proper_groups, OrientationRegion @@ -38,6 +38,7 @@ __all__ = [ "check_quaternion", "Quaternion", + "QuaternionNumpy", "Rotation", "von_mises", "Misorientation", diff --git a/orix/quaternion/quaternion.py b/orix/quaternion/quaternion.py index aab06c5b..2103fff7 100644 --- a/orix/quaternion/quaternion.py +++ b/orix/quaternion/quaternion.py @@ -20,6 +20,7 @@ import dask.array as da import numpy as np +import quaternion from orix.base import check, Object3d from orix.scalar import Scalar @@ -339,3 +340,83 @@ def _outer_dask(self, other, chunk_size=20): new_chunks = tuple(chunks1[:-1]) + tuple(chunks2[:-1]) + (-1,) return da.stack((a, b, c, d), axis=-1).rechunk(new_chunks) + + +class QuaternionNumpy(Quaternion): + r"""Basic quaternion object. + + Quaternions support the following mathematical operations: + + - Unary negation. + - Inversion. + - Multiplication with other quaternions and vectors. + + Attributes + ---------- + data : numpy.ndarray + The numpy array containing the quaternion data. + a, b, c, d : Scalar + The individual elements of each vector. + conj : Quaternion + The conjugate of this quaternion :math:`q^* = a - bi - cj - dk`. + """ + + @property + def conj(self): + q = quaternion.from_float_array(self.data).conj() + return Quaternion(quaternion.as_float_array(q)) + + def __mul__(self, other): + if isinstance(other, Quaternion): + q1 = quaternion.from_float_array(self.data) + q2 = quaternion.from_float_array(other.data) + return other.__class__(quaternion.as_float_array(q1 * q2)) + elif isinstance(other, Vector3d): + # check broadcast shape is correct before calculation, as + # quaternion.rotat_vectors will perform outer product + # this keeps current __mul__ broadcast behaviour + q1 = quaternion.from_float_array(self.data) + v = quaternion.as_vector_part( + (q1 * quaternion.from_vector_part(other.data)) * ~q1 + ) + if isinstance(other, Miller): + m = other.__class__(xyz=v, phase=other.phase) + m.coordinate_format = other.coordinate_format + return m + else: + return other.__class__(v) + return NotImplemented + + def outer(self, other): + """Compute the outer product of this quaternion and the other + quaternion or vector. + + Parameters + ---------- + other : orix.quaternion.Quaternion or orix.vector.Vector3d + + Returns + ------- + orix.quaternion.Quaternion or orix.vector.Vector3d + """ + + if isinstance(other, Quaternion): + q1 = quaternion.from_float_array(self.data) + q2 = quaternion.from_float_array(other.data) + # np.outer works with flattened array + q = np.outer(q1, q2).reshape(q1.shape + q2.shape) + return other.__class__(quaternion.as_float_array(q)) + elif isinstance(other, Vector3d): + q = quaternion.from_float_array(self.data) + v = quaternion.rotate_vectors(q, other.data) + if isinstance(other, Miller): + m = other.__class__(xyz=v, phase=other.phase) + m.coordinate_format = other.coordinate_format + return m + else: + return other.__class__(v) + else: + raise NotImplementedError( + "This operation is currently not avaliable in orix, please use outer " + "with `other` of type `Quaternion` or `Vector3d`" + ) diff --git a/orix/tests/quaternion/test_quaternion.py b/orix/tests/quaternion/test_quaternion.py index 516a0c85..74ff2dd1 100644 --- a/orix/tests/quaternion/test_quaternion.py +++ b/orix/tests/quaternion/test_quaternion.py @@ -21,7 +21,7 @@ import pytest from orix.base import DimensionError -from orix.quaternion import Quaternion, check_quaternion +from orix.quaternion import QuaternionNumpy as Quaternion, check_quaternion from orix.vector import Vector3d # fmt: off From 072ee836b2cdbf8123a26326fe710f56f4da7afd Mon Sep 17 00:00:00 2001 From: Paddy Harrison Date: Mon, 14 Feb 2022 13:42:35 +0100 Subject: [PATCH 2/8] add numpy-quaternion to setup.py --- setup.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/setup.py b/setup.py index 5631402a..d57d7037 100644 --- a/setup.py +++ b/setup.py @@ -56,8 +56,9 @@ "matplotlib-scalebar", "numba", "numpy", + "numpy-quaternion", "scipy", - "tqdm", + "tqdm" ], # fmt: on package_data={"": ["LICENSE", "README.rst", "readthedocs.yml"], "orix": ["*.py"]}, From da700eb9f9e5226d7cb258ccb4987ad201bfd34c Mon Sep 17 00:00:00 2001 From: Paddy Harrison Date: Thu, 17 Feb 2022 11:15:35 +0100 Subject: [PATCH 3/8] rm QNumpy class, add np-quat funcs to Quaternion --- orix/quaternion/quaternion.py | 148 +++++----------------------------- 1 file changed, 18 insertions(+), 130 deletions(-) diff --git a/orix/quaternion/quaternion.py b/orix/quaternion/quaternion.py index 2103fff7..a25921dc 100644 --- a/orix/quaternion/quaternion.py +++ b/orix/quaternion/quaternion.py @@ -90,39 +90,25 @@ def antipodal(self): @property def conj(self): - a = self.a.data - b, c, d = -self.b.data, -self.c.data, -self.d.data - q = np.stack((a, b, c, d), axis=-1) - return Quaternion(q) + q = quaternion.from_float_array(self.data).conj() + return Quaternion(quaternion.as_float_array(q)) def __invert__(self): return self.__class__(self.conj.data / (self.norm.data**2)[..., np.newaxis]) def __mul__(self, other): if isinstance(other, Quaternion): - sa, oa = self.a.data, other.a.data - sb, ob = self.b.data, other.b.data - sc, oc = self.c.data, other.c.data - sd, od = self.d.data, other.d.data - a = sa * oa - sb * ob - sc * oc - sd * od - b = sb * oa + sa * ob - sd * oc + sc * od - c = sc * oa + sd * ob + sa * oc - sb * od - d = sd * oa - sc * ob + sb * oc + sa * od - q = np.stack((a, b, c, d), axis=-1) - return other.__class__(q) + q1 = quaternion.from_float_array(self.data) + q2 = quaternion.from_float_array(other.data) + return other.__class__(quaternion.as_float_array(q1 * q2)) elif isinstance(other, Vector3d): - a, b, c, d = self.a.data, self.b.data, self.c.data, self.d.data - x, y, z = other.x.data, other.y.data, other.z.data - x_new = (a**2 + b**2 - c**2 - d**2) * x + 2 * ( - (a * c + b * d) * z + (b * c - a * d) * y - ) - y_new = (a**2 - b**2 + c**2 - d**2) * y + 2 * ( - (a * d + b * c) * x + (c * d - a * b) * z - ) - z_new = (a**2 - b**2 - c**2 + d**2) * z + 2 * ( - (a * b + c * d) * y + (b * d - a * c) * x + # check broadcast shape is correct before calculation, as + # quaternion.rotat_vectors will perform outer product + # this keeps current __mul__ broadcast behaviour + q1 = quaternion.from_float_array(self.data) + v = quaternion.as_vector_part( + (q1 * quaternion.from_vector_part(other.data)) * ~q1 ) - v = np.stack((x_new, y_new, z_new), axis=-1) if isinstance(other, Miller): m = other.__class__(xyz=v, phase=other.phase) m.coordinate_format = other.coordinate_format @@ -225,33 +211,15 @@ def outer(self, other): orix.quaternion.Quaternion or orix.vector.Vector3d """ - def e(x, y): - return np.multiply.outer(x, y) - if isinstance(other, Quaternion): - q = np.zeros(self.shape + other.shape + (4,), dtype=float) - sa, oa = self.data[..., 0], other.data[..., 0] - sb, ob = self.data[..., 1], other.data[..., 1] - sc, oc = self.data[..., 2], other.data[..., 2] - sd, od = self.data[..., 3], other.data[..., 3] - q[..., 0] = e(sa, oa) - e(sb, ob) - e(sc, oc) - e(sd, od) - q[..., 1] = e(sb, oa) + e(sa, ob) - e(sd, oc) + e(sc, od) - q[..., 2] = e(sc, oa) + e(sd, ob) + e(sa, oc) - e(sb, od) - q[..., 3] = e(sd, oa) - e(sc, ob) + e(sb, oc) + e(sa, od) - return other.__class__(q) + q1 = quaternion.from_float_array(self.data) + q2 = quaternion.from_float_array(other.data) + # np.outer works with flattened array + q = np.outer(q1, q2).reshape(q1.shape + q2.shape) + return other.__class__(quaternion.as_float_array(q)) elif isinstance(other, Vector3d): - a, b, c, d = self.a.data, self.b.data, self.c.data, self.d.data - x, y, z = other.x.data, other.y.data, other.z.data - x_new = e(a**2 + b**2 - c**2 - d**2, x) + 2 * ( - e(a * c + b * d, z) + e(b * c - a * d, y) - ) - y_new = e(a**2 - b**2 + c**2 - d**2, y) + 2 * ( - e(a * d + b * c, x) + e(c * d - a * b, z) - ) - z_new = e(a**2 - b**2 - c**2 + d**2, z) + 2 * ( - e(a * b + c * d, y) + e(b * d - a * c, x) - ) - v = np.stack((x_new, y_new, z_new), axis=-1) + q = quaternion.from_float_array(self.data) + v = quaternion.rotate_vectors(q, other.data) if isinstance(other, Miller): m = other.__class__(xyz=v, phase=other.phase) m.coordinate_format = other.coordinate_format @@ -340,83 +308,3 @@ def _outer_dask(self, other, chunk_size=20): new_chunks = tuple(chunks1[:-1]) + tuple(chunks2[:-1]) + (-1,) return da.stack((a, b, c, d), axis=-1).rechunk(new_chunks) - - -class QuaternionNumpy(Quaternion): - r"""Basic quaternion object. - - Quaternions support the following mathematical operations: - - - Unary negation. - - Inversion. - - Multiplication with other quaternions and vectors. - - Attributes - ---------- - data : numpy.ndarray - The numpy array containing the quaternion data. - a, b, c, d : Scalar - The individual elements of each vector. - conj : Quaternion - The conjugate of this quaternion :math:`q^* = a - bi - cj - dk`. - """ - - @property - def conj(self): - q = quaternion.from_float_array(self.data).conj() - return Quaternion(quaternion.as_float_array(q)) - - def __mul__(self, other): - if isinstance(other, Quaternion): - q1 = quaternion.from_float_array(self.data) - q2 = quaternion.from_float_array(other.data) - return other.__class__(quaternion.as_float_array(q1 * q2)) - elif isinstance(other, Vector3d): - # check broadcast shape is correct before calculation, as - # quaternion.rotat_vectors will perform outer product - # this keeps current __mul__ broadcast behaviour - q1 = quaternion.from_float_array(self.data) - v = quaternion.as_vector_part( - (q1 * quaternion.from_vector_part(other.data)) * ~q1 - ) - if isinstance(other, Miller): - m = other.__class__(xyz=v, phase=other.phase) - m.coordinate_format = other.coordinate_format - return m - else: - return other.__class__(v) - return NotImplemented - - def outer(self, other): - """Compute the outer product of this quaternion and the other - quaternion or vector. - - Parameters - ---------- - other : orix.quaternion.Quaternion or orix.vector.Vector3d - - Returns - ------- - orix.quaternion.Quaternion or orix.vector.Vector3d - """ - - if isinstance(other, Quaternion): - q1 = quaternion.from_float_array(self.data) - q2 = quaternion.from_float_array(other.data) - # np.outer works with flattened array - q = np.outer(q1, q2).reshape(q1.shape + q2.shape) - return other.__class__(quaternion.as_float_array(q)) - elif isinstance(other, Vector3d): - q = quaternion.from_float_array(self.data) - v = quaternion.rotate_vectors(q, other.data) - if isinstance(other, Miller): - m = other.__class__(xyz=v, phase=other.phase) - m.coordinate_format = other.coordinate_format - return m - else: - return other.__class__(v) - else: - raise NotImplementedError( - "This operation is currently not avaliable in orix, please use outer " - "with `other` of type `Quaternion` or `Vector3d`" - ) From 882cc164802a594b8309b2ef707db2586f588f83 Mon Sep 17 00:00:00 2001 From: Paddy Harrison Date: Thu, 17 Feb 2022 11:15:49 +0100 Subject: [PATCH 4/8] remove QuaternionNumpy import --- orix/tests/quaternion/test_quaternion.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/orix/tests/quaternion/test_quaternion.py b/orix/tests/quaternion/test_quaternion.py index 74ff2dd1..516a0c85 100644 --- a/orix/tests/quaternion/test_quaternion.py +++ b/orix/tests/quaternion/test_quaternion.py @@ -21,7 +21,7 @@ import pytest from orix.base import DimensionError -from orix.quaternion import QuaternionNumpy as Quaternion, check_quaternion +from orix.quaternion import Quaternion, check_quaternion from orix.vector import Vector3d # fmt: off From e332cfcf1a40cb7980d65df68b2a714241dd6e91 Mon Sep 17 00:00:00 2001 From: Paddy Harrison Date: Thu, 17 Feb 2022 11:34:20 +0100 Subject: [PATCH 5/8] update changelog --- CHANGELOG.rst | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 440df4ea..2542cad6 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -10,6 +10,12 @@ this project adheres to `Semantic Versioning `_ + for quaternion conjugation, quaternion-quaternion and quaternion-vector multiplication, + and quaternion-quaternion and quaternion-vector outer products. + 2022-02-14 - version 0.8.1 ========================== From 855bdebc94109fe568860c1a7ce9671275f460b2 Mon Sep 17 00:00:00 2001 From: Paddy Harrison Date: Thu, 17 Feb 2022 11:36:14 +0100 Subject: [PATCH 6/8] remove QuaternionNumpy from __init__ --- orix/quaternion/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/orix/quaternion/__init__.py b/orix/quaternion/__init__.py index d856ef69..3481472e 100644 --- a/orix/quaternion/__init__.py +++ b/orix/quaternion/__init__.py @@ -28,7 +28,7 @@ orientations. """ -from orix.quaternion.quaternion import check_quaternion, Quaternion, QuaternionNumpy +from orix.quaternion.quaternion import check_quaternion, Quaternion from orix.quaternion.rotation import Rotation, von_mises from orix.quaternion.orientation import Misorientation, Orientation from orix.quaternion.orientation_region import get_proper_groups, OrientationRegion From d8eba9116b1ee1a316b255d0fbd5066198d9f377 Mon Sep 17 00:00:00 2001 From: Paddy Harrison Date: Fri, 18 Feb 2022 11:39:44 +0100 Subject: [PATCH 7/8] add explicit q*q and q*v equations to docstring --- orix/quaternion/quaternion.py | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) diff --git a/orix/quaternion/quaternion.py b/orix/quaternion/quaternion.py index a25921dc..57520a12 100644 --- a/orix/quaternion/quaternion.py +++ b/orix/quaternion/quaternion.py @@ -40,6 +40,31 @@ class Quaternion(Object3d): - Inversion. - Multiplication with other quaternions and vectors. + Quaternion-quaternion multiplication for two quaternions + :math:`q_1 = (a_1, b_1, c_1, d_1)` + and :math:`q_2 = (a_2, b_2, c_2, d_2)` + with :math:`q_3 = (a_3, b_3, c_3, d_3) = q_1 * q_2` follows as: + + .. math:: + a_3 = (a_1 * a_2 - b_1 * b_2 - c_1 * c_2 - d_1 * d_2) + + b_3 = (a_1 * b_2 + b_1 * a_2 + c_1 * d_2 - d_1 * c_2) + + c_3 = (a_1 * c_2 - b_1 * d_2 + c_1 * a_2 + d_1 * b_2) + + d_3 = (a_1 * d_2 + b_1 * c_2 - c_1 * b_2 + d_1 * a_2) + + Quaternion-vector multiplication with a three-dimensional vector + :math:`v = (x, y, z)` calculates a rotated vector + :math:`v' = q * v * q^{-1}` and follows as: + + .. math:: + v'_x = x(a^2 + b^2 - c^2 - d^2) + 2z(a * c + b * d) + y(b * c - a * d) + + v'_y = y(a^2 - b^2 + c^2 - d^2) + 2x(a * d + b * c) + z(c * d - a * b) + + v'_z = z(a^2 - b^2 - c^2 + d^2) + 2y(a * b + c * d) + x(b * d - a * c) + Attributes ---------- data : numpy.ndarray From 4be8c88d9b7afecc15c3d30912e31d4a91263199 Mon Sep 17 00:00:00 2001 From: Paddy Harrison Date: Fri, 18 Feb 2022 12:18:02 +0100 Subject: [PATCH 8/8] remove QuaternionNumpy class --- orix/quaternion/__init__.py | 1 - setup.py | 2 +- 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/orix/quaternion/__init__.py b/orix/quaternion/__init__.py index 3481472e..e5d0e9c7 100644 --- a/orix/quaternion/__init__.py +++ b/orix/quaternion/__init__.py @@ -38,7 +38,6 @@ __all__ = [ "check_quaternion", "Quaternion", - "QuaternionNumpy", "Rotation", "von_mises", "Misorientation", diff --git a/setup.py b/setup.py index d57d7037..25cbd152 100644 --- a/setup.py +++ b/setup.py @@ -58,7 +58,7 @@ "numpy", "numpy-quaternion", "scipy", - "tqdm" + "tqdm", ], # fmt: on package_data={"": ["LICENSE", "README.rst", "readthedocs.yml"], "orix": ["*.py"]},