From 6a8783bab8a7f04ac35cc2dcb4a9130d288778cb Mon Sep 17 00:00:00 2001
From: Martin Diehl <m.diehl@mpie.de>
Date: Fri, 13 Dec 2019 15:46:40 +0100
Subject: [PATCH 01/11] these functions are static

class methods change are allowed to change the attributes of the class.
This is not what we do here
---
 python/damask/orientation.py | 78 ++++++++++++++++--------------------
 1 file changed, 34 insertions(+), 44 deletions(-)

diff --git a/python/damask/orientation.py b/python/damask/orientation.py
index 85a5c1866..6fae5b4c8 100644
--- a/python/damask/orientation.py
+++ b/python/damask/orientation.py
@@ -246,9 +246,8 @@ class Rotation:
     ################################################################################################
     # static constructors. The input data needs to follow the convention, options allow to
     # relax these convections
-    @classmethod
-    def fromQuaternion(cls,
-                       quaternion,
+    @staticmethod
+    def fromQuaternion(quaternion,
                        acceptHomomorph = False,
                        P = -1):
 
@@ -263,11 +262,10 @@ class Rotation:
       if not np.isclose(np.linalg.norm(qu), 1.0):
         raise ValueError('Quaternion is not of unit length.\n{} {} {} {}'.format(*qu))
 
-      return cls(qu)
+      return Rotation(qu)
       
-    @classmethod
-    def fromEulers(cls,
-                   eulers,
+    @staticmethod
+    def fromEulers(eulers,
                    degrees = False):
 
       eu = eulers if isinstance(eulers, np.ndarray) and eulers.dtype == np.dtype(float) \
@@ -276,11 +274,10 @@ class Rotation:
       if np.any(eu < 0.0) or np.any(eu > 2.0*np.pi) or eu[1] > np.pi:
         raise ValueError('Euler angles outside of [0..2π],[0..π],[0..2π].\n{} {} {}.'.format(*eu))
         
-      return cls(eu2qu(eu))
+      return Rotation(eu2qu(eu))
       
-    @classmethod
-    def fromAxisAngle(cls,
-                      angleAxis,
+    @staticmethod
+    def fromAxisAngle(angleAxis,
                       degrees = False,
                       normalise = False,
                       P = -1):
@@ -295,11 +292,10 @@ class Rotation:
       if not np.isclose(np.linalg.norm(ax[0:3]), 1.0):
         raise ValueError('Axis angle rotation axis is not of unit length.\n{} {} {}'.format(*ax[0:3]))
 
-      return cls(ax2qu(ax))
+      return Rotation(ax2qu(ax))
       
-    @classmethod
-    def fromBasis(cls,
-                  basis,
+    @staticmethod
+    def fromBasis(basis,
                   orthonormal = True,
                   reciprocal = False,
                  ):
@@ -318,18 +314,16 @@ class Rotation:
          or not np.isclose(np.dot(om[2],om[0]), 0.0):
         raise ValueError('matrix is not orthogonal.\n{}'.format(om))
 
-      return cls(om2qu(om))
+      return Rotation(om2qu(om))
       
-    @classmethod
-    def fromMatrix(cls,
-                   om,
+    @staticmethod
+    def fromMatrix(om,
                   ):
       
-      return cls.fromBasis(om)
+      return Rotation.fromBasis(om)
       
-    @classmethod
-    def fromRodrigues(cls,
-                      rodrigues,
+    @staticmethod
+    def fromRodrigues(rodrigues,
                       normalise = False,
                       P = -1):
         
@@ -342,35 +336,32 @@ class Rotation:
       if ro[3] < 0.0:
         raise ValueError('Rodriques rotation angle not positive.\n'.format(ro[3]))
         
-      return cls(ro2qu(ro))
+      return Rotation(ro2qu(ro))
       
-    @classmethod
-    def fromHomochoric(cls,
-                      homochoric,
-                      P = -1):
+    @staticmethod
+    def fromHomochoric(homochoric,
+                       P = -1):
         
       ho = homochoric if isinstance(homochoric, np.ndarray) and homochoric.dtype == np.dtype(float) \
                       else np.array(homochoric,dtype=float)
       if P > 0: ho *= -1                                                                            # convert from P=1 to P=-1
 
-      return cls(ho2qu(ho))
+      return Rotation(ho2qu(ho))
       
-    @classmethod
-    def fromCubochoric(cls,
-                      cubochoric,
-                      P = -1):
+    @staticmethod
+    def fromCubochoric(cubochoric,
+                       P = -1):
       
       cu = cubochoric if isinstance(cubochoric, np.ndarray) and cubochoric.dtype == np.dtype(float) \
                       else np.array(cubochoric,dtype=float)
       ho = cu2ho(cu)
       if P > 0: ho *= -1                                                                            # convert from P=1 to P=-1
         
-      return cls(ho2qu(ho))
+      return Rotation(ho2qu(ho))
 
 
-    @classmethod
-    def fromAverage(cls,
-                    rotations,
+    @staticmethod
+    def fromAverage(rotations,
                     weights = []):
       """
       Average rotation.
@@ -400,11 +391,11 @@ class Rotation:
             else M + r.asM() * n                                                                    # noqa add (multiples) of this rotation to average noqa
       eig, vec = np.linalg.eig(M/N)
 
-      return cls.fromQuaternion(np.real(vec.T[eig.argmax()]),acceptHomomorph = True)
+      return Rotation.fromQuaternion(np.real(vec.T[eig.argmax()]),acceptHomomorph = True)
 
 
-    @classmethod
-    def fromRandom(cls):
+    @staticmethod
+    def fromRandom():
       r = np.random.random(3) 
       A = np.sqrt(r[2]) 
       B = np.sqrt(1.0-r[2]) 
@@ -412,7 +403,7 @@ class Rotation:
       x = np.sin(2.0*np.pi*r[1])*B 
       y = np.cos(2.0*np.pi*r[1])*B 
       z = np.sin(2.0*np.pi*r[0])*A 
-      return cls.fromQuaternion([w,x,y,z],acceptHomomorph=True) 
+      return Rotation.fromQuaternion([w,x,y,z],acceptHomomorph=True) 
 
 
 
@@ -1180,9 +1171,8 @@ class Orientation:
     return color                                                                        
 
 
-  @classmethod
-  def fromAverage(cls,
-                  orientations,
+  @staticmethod
+  def fromAverage(orientations,
                   weights = []):
     """Create orientation from average of list of orientations."""
     if not all(isinstance(item, Orientation) for item in orientations):

From ecb8510217dbb3401f62988304575d9d63d6e4b9 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Tue, 14 Jan 2020 05:32:42 +0100
Subject: [PATCH 02/11] leaner (and hopefully faster) code

---
 python/damask/orientation.py | 54 +++++++++++++++++++-----------------
 1 file changed, 28 insertions(+), 26 deletions(-)

diff --git a/python/damask/orientation.py b/python/damask/orientation.py
index a86ba9331..b158495c8 100644
--- a/python/damask/orientation.py
+++ b/python/damask/orientation.py
@@ -3,7 +3,6 @@ import math
 import numpy as np
 
 from . import Lambert
-from .quaternion import Quaternion
 from .quaternion import P
 
 
@@ -49,10 +48,7 @@ class Rotation:
           Unit quaternion that follows the conventions. Use .fromQuaternion to perform a sanity check.
 
       """
-      if isinstance(quaternion,Quaternion):
-        self.quaternion = quaternion.copy()
-      else:
-        self.quaternion = Quaternion(q=quaternion[0],p=quaternion[1:4])
+      self.quaternion = quaternion.copy()
         
     def __copy__(self):
       """Copy."""
@@ -64,7 +60,7 @@ class Rotation:
     def __repr__(self):
       """Orientation displayed as unit quaternion, rotation matrix, and Bunge-Euler angles."""
       return '\n'.join([
-             '{}'.format(self.quaternion),
+             'Quaternion: (real={q:+.6f}, imag=<{p[0]:+.6f}, {p[1]:+.6f}, {p[2]:+.6f}>)'.format(*(self.quaternion)),
              'Matrix:\n{}'.format( '\n'.join(['\t'.join(list(map(str,self.asMatrix()[i,:]))) for i in range(3)]) ),
              'Bunge Eulers / deg: {}'.format('\t'.join(list(map(str,self.asEulers(degrees=True)))) ),
               ])
@@ -85,14 +81,20 @@ class Rotation:
 
       """
       if isinstance(other, Rotation):                                                               # rotate a rotation
-        return self.__class__(self.quaternion * other.quaternion).standardize()
+        self_q  = self.quaternion[0]
+        self_p  = self.quaternion[1:]
+        other_q = other.quaternion[0]
+        other_p = other.quaternion[1:]
+        r = np.append(self_q*other_q - np.dot(self_p,other_p),
+                      self_q*other_p + other_q*self_p + P * np.cross(self_p,other_p))
+        return __class__.fromQuaternion(r,acceptHomomorph=True)
       elif isinstance(other, np.ndarray):
         if other.shape == (3,):                                                                     # rotate a single (3)-vector
-          ( x, y, z) = self.quaternion.p
+          ( x, y, z) = self.quaternion[1:]
           (Vx,Vy,Vz) = other[0:3]
-          A = self.quaternion.q*self.quaternion.q - np.dot(self.quaternion.p,self.quaternion.p)
+          A = self.quaternion[0]**2.0 - np.dot(self.quaternion[1:],self.quaternion[1:])
           B = 2.0 * (x*Vx + y*Vy + z*Vz)
-          C = 2.0 * P*self.quaternion.q
+          C = 2.0 * P*self.quaternion[0]
 
           return np.array([
             A*Vx + B*x + C*(y*Vz - z*Vy),
@@ -106,11 +108,12 @@ class Rotation:
         else:
           return NotImplemented
       elif isinstance(other, tuple):                                                                # used to rotate a meshgrid-tuple
-          ( x, y, z) = self.quaternion.p
+          ( x, y, z) = self.quaternion[1:]
           (Vx,Vy,Vz) = other[0:3]
-          A = self.quaternion.q*self.quaternion.q - np.dot(self.quaternion.p,self.quaternion.p)
+          A = self.quaternion[0]**2.0 - np.dot(self.quaternion[1:],self.quaternion[1:])
           B = 2.0 * (x*Vx + y*Vy + z*Vz)
-          C = 2.0 * P*self.quaternion.q
+          C = 2.0 * P*self.quaternion[0]
+
 
           return np.array([
             A*Vx + B*x + C*(y*Vz - z*Vy),
@@ -123,7 +126,7 @@ class Rotation:
     
     def inverse(self):
       """In-place inverse rotation/backward rotation."""
-      self.quaternion.conjugate()
+      self.quaternion[1:] *= -1
       return self
       
     def inversed(self):
@@ -133,7 +136,7 @@ class Rotation:
 
     def standardize(self):
       """In-place quaternion representation with positive q."""
-      if self.quaternion.q < 0.0: self.quaternion.homomorph()
+      if self.quaternion[0] < 0.0: self.quaternion*=-1
       return self
 
     def standardized(self):
@@ -170,8 +173,7 @@ class Rotation:
     ################################################################################################
     # convert to different orientation representations (numpy arrays)
 
-    def asQuaternion(self,
-                     quaternion = False):
+    def asQuaternion(self):
       """
       Unit quaternion [q, p_1, p_2, p_3] unless quaternion == True: damask.quaternion object.
 
@@ -181,7 +183,7 @@ class Rotation:
           return quaternion as DAMASK object.
       
       """
-      return self.quaternion if quaternion else self.quaternion.asArray()
+      return self.quaternion
       
     def asEulers(self,
                  degrees = False):
@@ -194,7 +196,7 @@ class Rotation:
           return angles in degrees.
 
       """
-      eu = qu2eu(self.quaternion.asArray())
+      eu = qu2eu(self.quaternion)
       if degrees: eu = np.degrees(eu)
       return eu
     
@@ -212,13 +214,13 @@ class Rotation:
           return tuple of axis and angle.
       
       """
-      ax = qu2ax(self.quaternion.asArray())
+      ax = qu2ax(self.quaternion)
       if degrees: ax[3] = np.degrees(ax[3])
       return (ax[:3],np.degrees(ax[3])) if pair else ax
       
     def asMatrix(self):
       """Rotation matrix."""
-      return qu2om(self.quaternion.asArray())
+      return qu2om(self.quaternion)
 
     def asRodrigues(self,
                     vector = False):
@@ -231,16 +233,16 @@ class Rotation:
           return as actual Rodrigues--Frank vector, i.e. rotation axis scaled by tan(ω/2).
 
       """
-      ro = qu2ro(self.quaternion.asArray())
+      ro = qu2ro(self.quaternion)
       return ro[:3]*ro[3] if vector else ro
   
     def asHomochoric(self):
       """Homochoric vector: (h_1, h_2, h_3)."""
-      return qu2ho(self.quaternion.asArray())
+      return qu2ho(self.quaternion)
          
     def asCubochoric(self):
       """Cubochoric vector: (c_1, c_2, c_3)."""
-      return qu2cu(self.quaternion.asArray())
+      return qu2cu(self.quaternion)
       
     def asM(self):
       """
@@ -252,7 +254,7 @@ class Rotation:
       https://doi.org/10.2514/1.28949
 
       """
-      return self.quaternion.asM()
+      return np.outer(self.quaternion,self.quaternion)
       
 
     ################################################################################################
@@ -264,7 +266,7 @@ class Rotation:
                        acceptHomomorph = False,
                        P = -1):
 
-      qu =   quaternion if isinstance(quaternion, np.ndarray) and quaternion.dtype == np.dtype(float) \
+      qu =   quaternion if isinstance(quaternion,np.ndarray) and quaternion.dtype == np.dtype(float) \
                         else np.array(quaternion,dtype=float)
       if P > 0: qu[1:4] *= -1                                                                       # convert from P=1 to P=-1
       if qu[0] < 0.0:

From c7e92c672dd3f4c7b2fe82139e573745f8c2609a Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Tue, 14 Jan 2020 07:30:17 +0100
Subject: [PATCH 03/11] avoid code duplication

---
 python/damask/orientation.py | 34 ++++++++++------------------------
 1 file changed, 10 insertions(+), 24 deletions(-)

diff --git a/python/damask/orientation.py b/python/damask/orientation.py
index b158495c8..29cded8c3 100644
--- a/python/damask/orientation.py
+++ b/python/damask/orientation.py
@@ -85,21 +85,20 @@ class Rotation:
         self_p  = self.quaternion[1:]
         other_q = other.quaternion[0]
         other_p = other.quaternion[1:]
-        r = np.append(self_q*other_q - np.dot(self_p,other_p),
-                      self_q*other_p + other_q*self_p + P * np.cross(self_p,other_p))
-        return __class__.fromQuaternion(r,acceptHomomorph=True)
-      elif isinstance(other, np.ndarray):
-        if other.shape == (3,):                                                                     # rotate a single (3)-vector
-          ( x, y, z) = self.quaternion[1:]
-          (Vx,Vy,Vz) = other[0:3]
+        return __class__(np.append(self_q*other_q - np.dot(self_p,other_p),
+                                   self_q*other_p + other_q*self_p + P * np.cross(self_p,other_p))).standardize()
+      elif isinstance(other, (tuple,np.ndarray)):
+        if isinstance(other,tuple) or other.shape == (3,):                                          # rotate a single (3)-vector or meshgrid
           A = self.quaternion[0]**2.0 - np.dot(self.quaternion[1:],self.quaternion[1:])
-          B = 2.0 * (x*Vx + y*Vy + z*Vz)
+          B = 2.0 * (  self.quaternion[1]*other[0]
+                     + self.quaternion[2]*other[1]
+                     + self.quaternion[3]*other[2])
           C = 2.0 * P*self.quaternion[0]
 
           return np.array([
-            A*Vx + B*x + C*(y*Vz - z*Vy),
-            A*Vy + B*y + C*(z*Vx - x*Vz),
-            A*Vz + B*z + C*(x*Vy - y*Vx),
+            A*Vx + B*x + C*(self.quaternion[2]*other[2] - self.quaternion[3]*other[1]),
+            A*Vy + B*y + C*(self.quaternion[3]*other[0] - self.quaternion[1]*other[2]),
+            A*Vz + B*z + C*(self.quaternion[1]*other[1] - self.quaternion[2]*other[0]),
             ])
         elif other.shape == (3,3,):                                                                 # rotate a single (3x3)-matrix
            return np.dot(self.asMatrix(),np.dot(other,self.asMatrix().T))
@@ -107,19 +106,6 @@ class Rotation:
           raise NotImplementedError
         else:
           return NotImplemented
-      elif isinstance(other, tuple):                                                                # used to rotate a meshgrid-tuple
-          ( x, y, z) = self.quaternion[1:]
-          (Vx,Vy,Vz) = other[0:3]
-          A = self.quaternion[0]**2.0 - np.dot(self.quaternion[1:],self.quaternion[1:])
-          B = 2.0 * (x*Vx + y*Vy + z*Vz)
-          C = 2.0 * P*self.quaternion[0]
-
-
-          return np.array([
-            A*Vx + B*x + C*(y*Vz - z*Vy),
-            A*Vy + B*y + C*(z*Vx - x*Vz),
-            A*Vz + B*z + C*(x*Vy - y*Vx),
-            ])
       else:
         return NotImplemented
         

From 26c242b997438b7fb7804129db47fda30c21fa9a Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Tue, 14 Jan 2020 09:03:50 +0100
Subject: [PATCH 04/11] not needed anymore

conceptually, having the quaternion class is nice but performance-wise
it is a nightmare and prevents vectorization
---
 python/damask/orientation.py |   2 +-
 python/damask/quaternion.py  | 212 -----------------------------------
 2 files changed, 1 insertion(+), 213 deletions(-)
 delete mode 100644 python/damask/quaternion.py

diff --git a/python/damask/orientation.py b/python/damask/orientation.py
index d116cf044..12173cbfc 100644
--- a/python/damask/orientation.py
+++ b/python/damask/orientation.py
@@ -3,8 +3,8 @@ import math
 import numpy as np
 
 from . import Lambert
-from .quaternion import P
 
+P = -1
 
 ####################################################################################################
 class Rotation:
diff --git a/python/damask/quaternion.py b/python/damask/quaternion.py
deleted file mode 100644
index 203a398ef..000000000
--- a/python/damask/quaternion.py
+++ /dev/null
@@ -1,212 +0,0 @@
-import numpy as np
-
-P = -1                                                                                              # convention (see 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.
-    Definition 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 = float(q)
-      self.p = np.array(p,dtype=float)
-
-
-    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 negative operator"""  
-      return self * -1.0
-      
-      
-    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.q *= other
-        self.p *= 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)):
-        return self.__class__(q=self.q / other,
-                              p=self.p / other)
-      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
-        self.p /= 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)
-        return self
-      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):
-      """As list"""
-      return [self.q]+list(self.p)
-
-
-    def normalize(self):
-      """Normalizes in-place"""
-      d = self.magnitude()
-      if d > 0.0: self /= d
-      return self
-      
-    def normalized(self):
-      """Returns normalized copy"""
-      return self.copy().normalize()
-
-
-    def conjugate(self):
-      """Conjugates in-place"""
-      self.p = -self.p
-      return self
-      
-    def conjugated(self):
-      """Returns conjugated copy"""
-      return self.copy().conjugate()
-
-
-    def homomorph(self):
-      """Homomorphs in-place"""
-      self *= -1.0
-      return self
-      
-    def homomorphed(self):
-      """Returns homomorphed copy"""
-      return self.copy().homomorph()

From 8e5aa60f34f88f109046aae7900ab86d3ee41f8c Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Tue, 14 Jan 2020 09:37:59 +0100
Subject: [PATCH 05/11] fixed (and polished) printing

---
 python/damask/orientation.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/python/damask/orientation.py b/python/damask/orientation.py
index 12173cbfc..71c8efef9 100644
--- a/python/damask/orientation.py
+++ b/python/damask/orientation.py
@@ -60,9 +60,9 @@ class Rotation:
     def __repr__(self):
       """Orientation displayed as unit quaternion, rotation matrix, and Bunge-Euler angles."""
       return '\n'.join([
-             'Quaternion: (real={q:+.6f}, imag=<{p[0]:+.6f}, {p[1]:+.6f}, {p[2]:+.6f}>)'.format(*(self.quaternion)),
-             'Matrix:\n{}'.format( '\n'.join(['\t'.join(list(map(str,self.asMatrix()[i,:]))) for i in range(3)]) ),
-             'Bunge Eulers / deg: {}'.format('\t'.join(list(map(str,self.asEulers(degrees=True)))) ),
+             'Quaternion: (real={:.3f}, imag=<{:+.3f}, {:+.3f}, {:+.3f}>)'.format(*(self.quaternion)),
+             'Matrix:\n{}'.format(self.asMatrix()),
+             'Bunge Eulers / deg: ({:3.2f}, {:3.2f}, {:3.2f})'.format(*self.asEulers(degrees=True)),
               ])
               
     def __mul__(self, other):

From 06fc374359e5d1dfe73782d8350283013c52f7a8 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Tue, 14 Jan 2020 10:22:51 +0100
Subject: [PATCH 06/11] avoid overhead

length of random quaternion is one, no need to go via fromQuaternion but
direct initialization and standardization suffices
---
 python/damask/orientation.py | 9 ++++-----
 1 file changed, 4 insertions(+), 5 deletions(-)

diff --git a/python/damask/orientation.py b/python/damask/orientation.py
index 71c8efef9..16f5d66b1 100644
--- a/python/damask/orientation.py
+++ b/python/damask/orientation.py
@@ -399,11 +399,10 @@ class Rotation:
       r = np.random.random(3) 
       A = np.sqrt(r[2]) 
       B = np.sqrt(1.0-r[2]) 
-      w = np.cos(2.0*np.pi*r[0])*A 
-      x = np.sin(2.0*np.pi*r[1])*B 
-      y = np.cos(2.0*np.pi*r[1])*B 
-      z = np.sin(2.0*np.pi*r[0])*A 
-      return Rotation.fromQuaternion([w,x,y,z],acceptHomomorph=True) 
+      return Rotation(np.array([np.cos(2.0*np.pi*r[0])*A,
+                                np.sin(2.0*np.pi*r[1])*B,
+                                np.cos(2.0*np.pi*r[1])*B,
+                                np.sin(2.0*np.pi*r[0])*A])).standardize()
 
 
 

From 3afc011a58215003edcb93b58e8d7f491d8c9687 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Tue, 14 Jan 2020 10:34:13 +0100
Subject: [PATCH 07/11] fix for untested functionality

thanks to prospector
---
 python/damask/orientation.py | 11 ++++++-----
 1 file changed, 6 insertions(+), 5 deletions(-)

diff --git a/python/damask/orientation.py b/python/damask/orientation.py
index 16f5d66b1..8ea6499bf 100644
--- a/python/damask/orientation.py
+++ b/python/damask/orientation.py
@@ -85,8 +85,9 @@ class Rotation:
         self_p  = self.quaternion[1:]
         other_q = other.quaternion[0]
         other_p = other.quaternion[1:]
-        return __class__(np.append(self_q*other_q - np.dot(self_p,other_p),
-                                   self_q*other_p + other_q*self_p + P * np.cross(self_p,other_p))).standardize()
+        R = self.__class__(np.append(self_q*other_q - np.dot(self_p,other_p),
+                                     self_q*other_p + other_q*self_p + P * np.cross(self_p,other_p)))
+        return R.standardize()
       elif isinstance(other, (tuple,np.ndarray)):
         if isinstance(other,tuple) or other.shape == (3,):                                          # rotate a single (3)-vector or meshgrid
           A = self.quaternion[0]**2.0 - np.dot(self.quaternion[1:],self.quaternion[1:])
@@ -96,9 +97,9 @@ class Rotation:
           C = 2.0 * P*self.quaternion[0]
 
           return np.array([
-            A*Vx + B*x + C*(self.quaternion[2]*other[2] - self.quaternion[3]*other[1]),
-            A*Vy + B*y + C*(self.quaternion[3]*other[0] - self.quaternion[1]*other[2]),
-            A*Vz + B*z + C*(self.quaternion[1]*other[1] - self.quaternion[2]*other[0]),
+            A*self.quaternion[1] + B*other[0] + C*(self.quaternion[2]*other[2] - self.quaternion[3]*other[1]),
+            A*self.quaternion[2] + B*other[1] + C*(self.quaternion[3]*other[0] - self.quaternion[1]*other[2]),
+            A*self.quaternion[3] + B*other[2] + C*(self.quaternion[1]*other[1] - self.quaternion[2]*other[0]),
             ])
         elif other.shape == (3,3,):                                                                 # rotate a single (3x3)-matrix
            return np.dot(self.asMatrix(),np.dot(other,self.asMatrix().T))

From 18d6d364ec23fef0f6dd2edad0c1daca603aa6ae Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Tue, 14 Jan 2020 10:57:50 +0100
Subject: [PATCH 08/11] adjustment to tolerances

Rodrigues test is 'prone' to fail due to diverging tangens(pi/2).
Happens only with a very low chance
---
 python/tests/test_Rotation.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py
index e352e1c26..bb45c2906 100644
--- a/python/tests/test_Rotation.py
+++ b/python/tests/test_Rotation.py
@@ -7,7 +7,7 @@ import damask
 from damask import Rotation
 from damask import Orientation
    
-n = 1000
+n = 100000
 
 @pytest.fixture
 def default():
@@ -45,17 +45,17 @@ class TestRotation:
     def test_Homochoric(self,default):
         for rot in default:
             assert np.allclose(rot.asRodrigues(),
-                               Rotation.fromHomochoric(rot.asHomochoric()).asRodrigues(),rtol=5.e-5)
+                               Rotation.fromHomochoric(rot.asHomochoric()).asRodrigues(),rtol=1.e-4)
 
     def test_Cubochoric(self,default):
         for rot in default:
             assert np.allclose(rot.asHomochoric(),
-                               Rotation.fromCubochoric(rot.asCubochoric()).asHomochoric(),rtol=5.e-5)
+                               Rotation.fromCubochoric(rot.asCubochoric()).asHomochoric())
 
     def test_Quaternion(self,default):
         for rot in default:
             assert np.allclose(rot.asCubochoric(),
-                               Rotation.fromQuaternion(rot.asQuaternion()).asCubochoric(),rtol=5.e-5)
+                               Rotation.fromQuaternion(rot.asQuaternion()).asCubochoric())
 
 
     @pytest.mark.parametrize('model',['Bain','KS','GT','GT_prime','NW','Pitsch'])

From 8a674d02792dfc666cf6ff12b54f4acf2053f586 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Tue, 14 Jan 2020 22:30:08 +0100
Subject: [PATCH 09/11] 2 bugfixes:

- mixed up other and me when removing quaternion class
- 1000 random orientations are enough for testing
---
 python/damask/orientation.py  | 6 +++---
 python/tests/test_Rotation.py | 2 +-
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/python/damask/orientation.py b/python/damask/orientation.py
index 8ea6499bf..f830024f8 100644
--- a/python/damask/orientation.py
+++ b/python/damask/orientation.py
@@ -97,9 +97,9 @@ class Rotation:
           C = 2.0 * P*self.quaternion[0]
 
           return np.array([
-            A*self.quaternion[1] + B*other[0] + C*(self.quaternion[2]*other[2] - self.quaternion[3]*other[1]),
-            A*self.quaternion[2] + B*other[1] + C*(self.quaternion[3]*other[0] - self.quaternion[1]*other[2]),
-            A*self.quaternion[3] + B*other[2] + C*(self.quaternion[1]*other[1] - self.quaternion[2]*other[0]),
+            A*other[0] + B*self.quaternion[1] + C*(self.quaternion[2]*other[2] - self.quaternion[3]*other[1]),
+            A*other[1] + B*self.quaternion[2] + C*(self.quaternion[3]*other[0] - self.quaternion[1]*other[2]),
+            A*other[2] + B*self.quaternion[3] + C*(self.quaternion[1]*other[1] - self.quaternion[2]*other[0]),
             ])
         elif other.shape == (3,3,):                                                                 # rotate a single (3x3)-matrix
            return np.dot(self.asMatrix(),np.dot(other,self.asMatrix().T))
diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py
index bb45c2906..6c9cdb8a5 100644
--- a/python/tests/test_Rotation.py
+++ b/python/tests/test_Rotation.py
@@ -7,7 +7,7 @@ import damask
 from damask import Rotation
 from damask import Orientation
    
-n = 100000
+n = 1000
 
 @pytest.fixture
 def default():

From 2bc93b643aa6f96d83a64562fb95f62729dfc043 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Tue, 14 Jan 2020 23:17:47 +0100
Subject: [PATCH 10/11] testing IPF color values for cube

---
 python/tests/test_Rotation.py | 11 +++++++++++
 1 file changed, 11 insertions(+)

diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py
index 6c9cdb8a5..e0a454590 100644
--- a/python/tests/test_Rotation.py
+++ b/python/tests/test_Rotation.py
@@ -1,4 +1,5 @@
 import os
+from itertools import permutations
 
 import pytest
 import numpy as np
@@ -58,6 +59,16 @@ class TestRotation:
                                Rotation.fromQuaternion(rot.asQuaternion()).asCubochoric())
 
 
+    @pytest.mark.parametrize('color',[{'label':'red',  'RGB':[1,0,0],'direction':[0,0,1]},
+                                      {'label':'green','RGB':[0,1,0],'direction':[0,1,1]},
+                                      {'label':'blue', 'RGB':[0,0,1],'direction':[1,1,1]}])
+    @pytest.mark.parametrize('lattice',['fcc','bcc'])
+    def test_IPF_cubic(self,default,color,lattice):
+        cube = damask.Orientation(damask.Rotation(),lattice)
+        for direction in set(permutations(np.array(color['direction']))):
+            assert np.allclose(cube.IPFcolor(direction),np.array(color['RGB']))
+
+
     @pytest.mark.parametrize('model',['Bain','KS','GT','GT_prime','NW','Pitsch'])
     @pytest.mark.parametrize('lattice',['fcc','bcc'])
     def test_relationship_forward_backward(self,model,lattice):

From 2ba0752b71dd8fc69e03b3ba96a51e0ab932f730 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Tue, 14 Jan 2020 23:56:52 +0100
Subject: [PATCH 11/11] ensure that symmetrically equivalent rotations have the
 same IPF color

---
 python/tests/test_Rotation.py | 9 +++++++++
 1 file changed, 9 insertions(+)

diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py
index e0a454590..012c95469 100644
--- a/python/tests/test_Rotation.py
+++ b/python/tests/test_Rotation.py
@@ -7,6 +7,7 @@ import numpy as np
 import damask
 from damask import Rotation
 from damask import Orientation
+from damask import Lattice
    
 n = 1000
 
@@ -68,6 +69,14 @@ class TestRotation:
         for direction in set(permutations(np.array(color['direction']))):
             assert np.allclose(cube.IPFcolor(direction),np.array(color['RGB']))
 
+    @pytest.mark.parametrize('lattice',Lattice.lattices)
+    def test_IPF(self,lattice):
+        direction = np.random.random(3)*2.0-1
+        for rot in [Rotation.fromRandom() for r in range(n//100)]:
+            R = damask.Orientation(rot,lattice)
+            color = R.IPFcolor(direction)
+            for equivalent in R.equivalentOrientations():
+                assert np.allclose(color,R.IPFcolor(direction))
 
     @pytest.mark.parametrize('model',['Bain','KS','GT','GT_prime','NW','Pitsch'])
     @pytest.mark.parametrize('lattice',['fcc','bcc'])