From fa182004470eb67f0e7f8a4a253b4b43787b6dce Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 17 Apr 2019 08:46:03 +0200 Subject: [PATCH] quaternion is now in separate module avoid long modules with multiple, only loosely related classes --- python/damask/orientation.py | 212 +---------------------------------- python/damask/quaternion.py | 210 ++++++++++++++++++++++++++++++++++ 2 files changed, 213 insertions(+), 209 deletions(-) create mode 100644 python/damask/quaternion.py diff --git a/python/damask/orientation.py b/python/damask/orientation.py index 007aecfd9..349c368fa 100644 --- a/python/damask/orientation.py +++ b/python/damask/orientation.py @@ -3,214 +3,8 @@ import math import numpy as np from . import Lambert - -P = -1 - -#################################################################################################### -class Quaternion: - u""" - Quaternion with basic operations - - q is the real part, p = (x, y, z) are the imaginary parts. - Defintion of multiplication depends on variable P, P ∉ {-1,1}. - """ - - def __init__(self, - q = 0.0, - p = np.zeros(3,dtype=float)): - """Initializes to identity unless specified""" - self.q = q - self.p = np.array(p) - - - def __copy__(self): - """Copy""" - return self.__class__(q=self.q, - p=self.p.copy()) - - copy = __copy__ - - - def __iter__(self): - """Components""" - return iter(self.asList()) - - def __repr__(self): - """Readable string""" - return 'Quaternion: (real={q:+.6f}, imag=<{p[0]:+.6f}, {p[1]:+.6f}, {p[2]:+.6f}>)'.format(q=self.q,p=self.p) - - - def __add__(self, other): - """Addition""" - if isinstance(other, Quaternion): - return self.__class__(q=self.q + other.q, - p=self.p + other.p) - else: - return NotImplemented - - def __iadd__(self, other): - """In-place addition""" - if isinstance(other, Quaternion): - self.q += other.q - self.p += other.p - return self - else: - return NotImplemented - - def __pos__(self): - """Unary positive operator""" - return self - - - def __sub__(self, other): - """Subtraction""" - if isinstance(other, Quaternion): - return self.__class__(q=self.q - other.q, - p=self.p - other.p) - else: - return NotImplemented - - def __isub__(self, other): - """In-place subtraction""" - if isinstance(other, Quaternion): - self.q -= other.q - self.p -= other.p - return self - else: - return NotImplemented - - def __neg__(self): - """Unary positive operator""" - self.q *= -1.0 - self.p *= -1.0 - return self - - - def __mul__(self, other): - """Multiplication with quaternion or scalar""" - if isinstance(other, Quaternion): - return self.__class__(q=self.q*other.q - np.dot(self.p,other.p), - p=self.q*other.p + other.q*self.p + P * np.cross(self.p,other.p)) - elif isinstance(other, (int, float)): - return self.__class__(q=self.q*other, - p=self.p*other) - else: - return NotImplemented - - def __imul__(self, other): - """In-place multiplication with quaternion or scalar""" - if isinstance(other, Quaternion): - self.q = self.q*other.q - np.dot(self.p,other.p) - self.p = self.q*other.p + other.q*self.p + P * np.cross(self.p,other.p) - return self - elif isinstance(other, (int, float)): - self *= other - return self - else: - return NotImplemented - - - def __truediv__(self, other): - """Divsion with quaternion or scalar""" - if isinstance(other, Quaternion): - s = other.conjugate()/abs(other)**2. - return self.__class__(q=self.q * s, - p=self.p * s) - elif isinstance(other, (int, float)): - self.q /= other - self.p /= other - return self - else: - return NotImplemented - - def __itruediv__(self, other): - """In-place divsion with quaternion or scalar""" - if isinstance(other, Quaternion): - s = other.conjugate()/abs(other)**2. - self *= s - return self - elif isinstance(other, (int, float)): - self.q /= other - return self - else: - return NotImplemented - - - def __pow__(self, exponent): - """Power""" - if isinstance(exponent, (int, float)): - omega = np.acos(self.q) - return self.__class__(q= np.cos(exponent*omega), - p=self.p * np.sin(exponent*omega)/np.sin(omega)) - else: - return NotImplemented - - def __ipow__(self, exponent): - """In-place power""" - if isinstance(exponent, (int, float)): - omega = np.acos(self.q) - self.q = np.cos(exponent*omega) - self.p *= np.sin(exponent*omega)/np.sin(omega) - else: - return NotImplemented - - - def __abs__(self): - """Norm""" - return math.sqrt(self.q ** 2 + np.dot(self.p,self.p)) - - magnitude = __abs__ - - - def __eq__(self,other): - """Equal (sufficiently close) to each other""" - return np.isclose(( self-other).magnitude(),0.0) \ - or np.isclose((-self-other).magnitude(),0.0) - - def __ne__(self,other): - """Not equal (sufficiently close) to each other""" - return not self.__eq__(other) - - - def asM(self): - """Intermediate representation useful for quaternion averaging (see F. Landis Markley et al.)""" - return np.outer(self.asArray(),self.asArray()) - - def asArray(self): - """As numpy array""" - return np.array((self.q,self.p[0],self.p[1],self.p[2])) - - def asList(self): - return [self.q]+list(self.p) - - def normalize(self): - d = self.magnitude() - if d > 0.0: - self.q /= d - self.p /= d - return self - - def normalized(self): - return self.copy().normalize() - - - def conjugate(self): - self.p = -self.p - return self - - def conjugated(self): - return self.copy().conjugate() - - - def homomorph(self): - if self.q < 0.0: - self.q = -self.q - self.p = -self.p - return self - - def homomorphed(self): - return self.copy().homomorph() - +from quaternion import Quaternion +from quaternion import P as P #################################################################################################### @@ -488,7 +282,7 @@ class Rotation: for i,(r,n) in enumerate(zip(rotations,weights)): M = r.asM() * n if i == 0 \ - else M + r.asM() * n # noqa add (multiples) of this rotation to average noqa + else M + r.asM() * n # noqa add (multiples) of this rotation to average noqa eig, vec = np.linalg.eig(M/N) return Rotation.fromQuaternion(np.real(vec.T[eig.argmax()]),acceptHomomorph = True) diff --git a/python/damask/quaternion.py b/python/damask/quaternion.py new file mode 100644 index 000000000..7c98a6d1c --- /dev/null +++ b/python/damask/quaternion.py @@ -0,0 +1,210 @@ +# -*- coding: UTF-8 no BOM -*- + +import numpy as np + +P = -1 # convention (sed DOI:10.1088/0965-0393/23/8/083501) + +#################################################################################################### +class Quaternion: + u""" + Quaternion with basic operations + + q is the real part, p = (x, y, z) are the imaginary parts. + Defintion of multiplication depends on variable P, P ∈ {-1,1}. + """ + + def __init__(self, + q = 0.0, + p = np.zeros(3,dtype=float)): + """Initializes to identity unless specified""" + self.q = q + self.p = np.array(p) + + + def __copy__(self): + """Copy""" + return self.__class__(q=self.q, + p=self.p.copy()) + + copy = __copy__ + + + def __iter__(self): + """Components""" + return iter(self.asList()) + + def __repr__(self): + """Readable string""" + return 'Quaternion: (real={q:+.6f}, imag=<{p[0]:+.6f}, {p[1]:+.6f}, {p[2]:+.6f}>)'.format(q=self.q,p=self.p) + + + def __add__(self, other): + """Addition""" + if isinstance(other, Quaternion): + return self.__class__(q=self.q + other.q, + p=self.p + other.p) + else: + return NotImplemented + + def __iadd__(self, other): + """In-place addition""" + if isinstance(other, Quaternion): + self.q += other.q + self.p += other.p + return self + else: + return NotImplemented + + def __pos__(self): + """Unary positive operator""" + return self + + + def __sub__(self, other): + """Subtraction""" + if isinstance(other, Quaternion): + return self.__class__(q=self.q - other.q, + p=self.p - other.p) + else: + return NotImplemented + + def __isub__(self, other): + """In-place subtraction""" + if isinstance(other, Quaternion): + self.q -= other.q + self.p -= other.p + return self + else: + return NotImplemented + + def __neg__(self): + """Unary positive operator""" + self.q *= -1.0 + self.p *= -1.0 + return self + + + def __mul__(self, other): + """Multiplication with quaternion or scalar""" + if isinstance(other, Quaternion): + return self.__class__(q=self.q*other.q - np.dot(self.p,other.p), + p=self.q*other.p + other.q*self.p + P * np.cross(self.p,other.p)) + elif isinstance(other, (int, float)): + return self.__class__(q=self.q*other, + p=self.p*other) + else: + return NotImplemented + + def __imul__(self, other): + """In-place multiplication with quaternion or scalar""" + if isinstance(other, Quaternion): + self.q = self.q*other.q - np.dot(self.p,other.p) + self.p = self.q*other.p + other.q*self.p + P * np.cross(self.p,other.p) + return self + elif isinstance(other, (int, float)): + self *= other + return self + else: + return NotImplemented + + + def __truediv__(self, other): + """Divsion with quaternion or scalar""" + if isinstance(other, Quaternion): + s = other.conjugate()/abs(other)**2. + return self.__class__(q=self.q * s, + p=self.p * s) + elif isinstance(other, (int, float)): + self.q /= other + self.p /= other + return self + else: + return NotImplemented + + def __itruediv__(self, other): + """In-place divsion with quaternion or scalar""" + if isinstance(other, Quaternion): + s = other.conjugate()/abs(other)**2. + self *= s + return self + elif isinstance(other, (int, float)): + self.q /= other + return self + else: + return NotImplemented + + + def __pow__(self, exponent): + """Power""" + if isinstance(exponent, (int, float)): + omega = np.acos(self.q) + return self.__class__(q= np.cos(exponent*omega), + p=self.p * np.sin(exponent*omega)/np.sin(omega)) + else: + return NotImplemented + + def __ipow__(self, exponent): + """In-place power""" + if isinstance(exponent, (int, float)): + omega = np.acos(self.q) + self.q = np.cos(exponent*omega) + self.p *= np.sin(exponent*omega)/np.sin(omega) + else: + return NotImplemented + + + def __abs__(self): + """Norm""" + return np.sqrt(self.q ** 2 + np.dot(self.p,self.p)) + + magnitude = __abs__ + + + def __eq__(self,other): + """Equal (sufficiently close) to each other""" + return np.isclose(( self-other).magnitude(),0.0) \ + or np.isclose((-self-other).magnitude(),0.0) + + def __ne__(self,other): + """Not equal (sufficiently close) to each other""" + return not self.__eq__(other) + + + def asM(self): + """Intermediate representation useful for quaternion averaging (see F. Landis Markley et al.)""" + return np.outer(self.asArray(),self.asArray()) + + def asArray(self): + """As numpy array""" + return np.array((self.q,self.p[0],self.p[1],self.p[2])) + + def asList(self): + return [self.q]+list(self.p) + + def normalize(self): + d = self.magnitude() + if d > 0.0: + self.q /= d + self.p /= d + return self + + def normalized(self): + return self.copy().normalize() + + + def conjugate(self): + self.p = -self.p + return self + + def conjugated(self): + return self.copy().conjugate() + + + def homomorph(self): + if self.q < 0.0: + self.q = -self.q + self.p = -self.p + return self + + def homomorphed(self): + return self.copy().homomorph()