Merge branch 'python-improvements' into 'development'

Python improvements

See merge request damask/DAMASK!315
This commit is contained in:
Philip Eisenlohr 2021-01-13 17:25:25 +01:00
commit 8182c9c54e
12 changed files with 290 additions and 110 deletions

4
.gitmodules vendored
View File

@ -1,5 +1,5 @@
[submodule "PRIVATE"] [submodule "PRIVATE"]
path = PRIVATE path = PRIVATE
url = ../PRIVATE.git url = ../PRIVATE.git
branch = master branch = master
shallow = true shallow = true

View File

@ -1,4 +1,4 @@
Copyright 2011-20 Max-Planck-Institut für Eisenforschung GmbH Copyright 2011-21 Max-Planck-Institut für Eisenforschung GmbH
DAMASK is free software: you can redistribute it and/or modify DAMASK is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by it under the terms of the GNU General Public License as published by

@ -1 +1 @@
Subproject commit 7846c71126705cc5d41dd79f2d595f4864434068 Subproject commit 1d0c95c5c1c0e7e6f57bdfc94b695e47a6ad6c60

View File

@ -1,3 +1,4 @@
import copy
from io import StringIO from io import StringIO
import abc import abc
@ -35,6 +36,14 @@ class Config(dict):
output.seek(0) output.seek(0)
return ''.join(output.readlines()) return ''.join(output.readlines())
def __copy__(self):
"""Create deep copy."""
return copy.deepcopy(self)
copy = __copy__
@classmethod @classmethod
def load(cls,fname): def load(cls,fname):
""" """
@ -52,6 +61,7 @@ class Config(dict):
fhandle = fname fhandle = fname
return cls(yaml.safe_load(fhandle)) return cls(yaml.safe_load(fhandle))
def save(self,fname,**kwargs): def save(self,fname,**kwargs):
""" """
Save to yaml file. Save to yaml file.
@ -89,12 +99,37 @@ class Config(dict):
fhandle.write(yaml.dump(self,Dumper=NiceDumper,**kwargs)) fhandle.write(yaml.dump(self,Dumper=NiceDumper,**kwargs))
def add(self,d):
"""
Add dictionary.
d : dict
Dictionary to append.
"""
duplicate = self.copy()
duplicate.update(d)
return duplicate
def delete(self,key):
"""
Delete item.
key : str or scalar
Label of the key to remove.
"""
duplicate = self.copy()
del duplicate[key]
return duplicate
@property @property
@abc.abstractmethod @abc.abstractmethod
def is_complete(self): def is_complete(self):
"""Check for completeness.""" """Check for completeness."""
pass pass
@property @property
@abc.abstractmethod @abc.abstractmethod
def is_valid(self): def is_valid(self):

View File

@ -1,5 +1,3 @@
import copy
import numpy as np import numpy as np
from . import Config from . import Config
@ -13,11 +11,10 @@ class ConfigMaterial(Config):
'homogenization': {}, 'homogenization': {},
'phase': {}} 'phase': {}}
def __init__(self,d={}): def __init__(self,d=_defaults):
"""Initialize object with default dictionary keys.""" """Initialize object with default dictionary keys."""
super().__init__(d) super().__init__(d)
for k,v in self._defaults.items():
if k not in self: self[k] = v
def save(self,fname='material.yaml',**kwargs): def save(self,fname='material.yaml',**kwargs):
""" """
@ -204,7 +201,7 @@ class ConfigMaterial(Config):
Limit renaming to selected constituents. Limit renaming to selected constituents.
""" """
dup = copy.deepcopy(self) dup = self.copy()
for i,m in enumerate(dup['material']): for i,m in enumerate(dup['material']):
if ID and i not in ID: continue if ID and i not in ID: continue
for c in m['constituents']: for c in m['constituents']:
@ -228,7 +225,7 @@ class ConfigMaterial(Config):
Limit renaming to selected homogenization IDs. Limit renaming to selected homogenization IDs.
""" """
dup = copy.deepcopy(self) dup = self.copy()
for i,m in enumerate(dup['material']): for i,m in enumerate(dup['material']):
if ID and i not in ID: continue if ID and i not in ID: continue
try: try:
@ -289,7 +286,7 @@ class ConfigMaterial(Config):
c = [{} for _ in range(length)] if constituents is None else \ c = [{} for _ in range(length)] if constituents is None else \
[{'constituents':u} for u in ConfigMaterial._constituents(**constituents)] [{'constituents':u} for u in ConfigMaterial._constituents(**constituents)]
if len(c) == 1: c = [copy.deepcopy(c[0]) for _ in range(length)] if len(c) == 1: c = [c[0] for _ in range(length)]
if length != 1 and length != len(c): if length != 1 and length != len(c):
raise ValueError('Cannot add entries of different length') raise ValueError('Cannot add entries of different length')
@ -301,7 +298,7 @@ class ConfigMaterial(Config):
else: else:
for i in range(len(c)): for i in range(len(c)):
c[i][k] = v c[i][k] = v
dup = copy.deepcopy(self) dup = self.copy()
dup['material'] = dup['material'] + c if 'material' in dup else c dup['material'] = dup['material'] + c if 'material' in dup else c
return dup return dup

View File

@ -57,13 +57,10 @@ class Grid:
def __copy__(self): def __copy__(self):
"""Copy grid.""" """Create deep copy."""
return copy.deepcopy(self) return copy.deepcopy(self)
copy = __copy__
def copy(self):
"""Copy grid."""
return self.__copy__()
def diff(self,other): def diff(self,other):
@ -766,24 +763,19 @@ class Grid:
if fill is None: fill = np.nanmax(self.material) + 1 if fill is None: fill = np.nanmax(self.material) + 1
dtype = float if np.isnan(fill) or int(fill) != fill or self.material.dtype==np.float else int dtype = float if np.isnan(fill) or int(fill) != fill or self.material.dtype==np.float else int
Eulers = R.as_Euler_angles(degrees=True) material = self.material
material_in = self.material.copy()
# These rotations are always applied in the reference coordinate system, i.e. (z,x,z) not (z,x',z'') # These rotations are always applied in the reference coordinate system, i.e. (z,x,z) not (z,x',z'')
# see https://www.cs.utexas.edu/~theshark/courses/cs354/lectures/cs354-14.pdf # see https://www.cs.utexas.edu/~theshark/courses/cs354/lectures/cs354-14.pdf
for angle,axes in zip(Eulers[::-1], [(0,1),(1,2),(0,1)]): for angle,axes in zip(R.as_Euler_angles(degrees=True)[::-1], [(0,1),(1,2),(0,1)]):
material_out = ndimage.rotate(material_in,angle,axes,order=0, material_temp = ndimage.rotate(material,angle,axes,order=0,prefilter=False,output=dtype,cval=fill)
prefilter=False,output=dtype,cval=fill) # avoid scipy interpolation errors for rotations close to multiples of 90°
if np.prod(material_in.shape) == np.prod(material_out.shape): material = material_temp if np.prod(material_temp.shape) != np.prod(material.shape) else \
# avoid scipy interpolation errors for rotations close to multiples of 90° np.rot90(material,k=np.rint(angle/90.).astype(int),axes=axes)
material_in = np.rot90(material_in,k=np.rint(angle/90.).astype(int),axes=axes)
else:
material_in = material_out
origin = self.origin-(np.asarray(material_in.shape)-self.cells)*.5 * self.size/self.cells origin = self.origin-(np.asarray(material.shape)-self.cells)*.5 * self.size/self.cells
return Grid(material = material_in, return Grid(material = material,
size = self.size/self.cells*np.asarray(material_in.shape), size = self.size/self.cells*np.asarray(material.shape),
origin = origin, origin = origin,
comments = self.comments+[util.execution_stamp('Grid','rotate')], comments = self.comments+[util.execution_stamp('Grid','rotate')],
) )

View File

@ -199,7 +199,7 @@ class Orientation(Rotation):
def __copy__(self,**kwargs): def __copy__(self,**kwargs):
"""Copy.""" """Create deep copy."""
return self.__class__(rotation=kwargs['rotation'] if 'rotation' in kwargs else self.quaternion, return self.__class__(rotation=kwargs['rotation'] if 'rotation' in kwargs else self.quaternion,
lattice =kwargs['lattice'] if 'lattice' in kwargs else self.lattice lattice =kwargs['lattice'] if 'lattice' in kwargs else self.lattice
if self.lattice is not None else self.family, if self.lattice is not None else self.family,
@ -225,30 +225,42 @@ class Orientation(Rotation):
Orientation to check for equality. Orientation to check for equality.
""" """
return super().__eq__(other) \ matching_type = all([hasattr(other,attr) and getattr(self,attr) == getattr(other,attr)
and hasattr(other, 'family') and self.family == other.family \ for attr in ['family','lattice','parameters']])
and hasattr(other, 'lattice') and self.lattice == other.lattice \ return np.logical_and(super().__eq__(other),matching_type)
and hasattr(other, 'parameters') and self.parameters == other.parameters
def __ne__(self,other):
def __matmul__(self,other):
""" """
Rotation of vector, second or fourth order tensor, or rotation object. Not equal to other.
Parameters Parameters
---------- ----------
other : numpy.ndarray, Rotation, or Orientation other : Orientation
Vector, second or fourth order tensor, or rotation object that is rotated. Orientation to check for equality.
"""
return np.logical_not(self==other)
def __mul__(self,other):
"""
Compose this orientation with other.
Parameters
----------
other : Rotation or Orientation
Object for composition.
Returns Returns
------- -------
other_rot : numpy.ndarray or Rotation composition : Orientation
Rotated vector, second or fourth order tensor, or rotation object. Compound rotation self*other, i.e. first other then self rotation.
""" """
return self.copy(rotation=Rotation.__matmul__(self,Rotation(other.quaternion))) \ if isinstance(other,Orientation) or isinstance(other,Rotation):
if isinstance(other,self.__class__) else \ return self.copy(rotation=Rotation.__mul__(self,Rotation(other.quaternion)))
Rotation.__matmul__(self,other) else:
raise TypeError('Use "O@b", i.e. matmul, to apply Orientation "O" to object "b"')
@classmethod @classmethod
@ -429,7 +441,7 @@ class Orientation(Rotation):
raise ValueError('Missing crystal symmetry') raise ValueError('Missing crystal symmetry')
o = self.symmetry_operations.broadcast_to(self.symmetry_operations.shape+self.shape,mode='right') o = self.symmetry_operations.broadcast_to(self.symmetry_operations.shape+self.shape,mode='right')
return self.copy(rotation=o@Rotation(self.quaternion).broadcast_to(o.shape,mode='left')) return self.copy(rotation=o*Rotation(self.quaternion).broadcast_to(o.shape,mode='left'))
@property @property
@ -608,7 +620,7 @@ class Orientation(Rotation):
o,lattice = self.relation_operations(model,return_lattice=True) o,lattice = self.relation_operations(model,return_lattice=True)
target = Orientation(lattice=lattice) target = Orientation(lattice=lattice)
o = o.broadcast_to(o.shape+self.shape,mode='right') o = o.broadcast_to(o.shape+self.shape,mode='right')
return self.copy(rotation=o@Rotation(self.quaternion).broadcast_to(o.shape,mode='left'), return self.copy(rotation=o*Rotation(self.quaternion).broadcast_to(o.shape,mode='left'),
lattice=lattice, lattice=lattice,
b = self.b if target.ratio['b'] is None else self.a*target.ratio['b'], b = self.b if target.ratio['b'] is None else self.a*target.ratio['b'],
c = self.c if target.ratio['c'] is None else self.a*target.ratio['c'], c = self.c if target.ratio['c'] is None else self.a*target.ratio['c'],

View File

@ -35,6 +35,11 @@ class Rotation:
- b = Q @ a - b = Q @ a
- b = np.dot(Q.as_matrix(),a) - b = np.dot(Q.as_matrix(),a)
Compound rotations R1 (first) and R2 (second):
- R = R2 * R1
- R = Rotation.from_matrix(np.dot(R2.as_matrix(),R1.as_matrix())
References References
---------- ----------
D. Rowenhorst et al., Modelling and Simulation in Materials Science and Engineering 23:083501, 2015 D. Rowenhorst et al., Modelling and Simulation in Materials Science and Engineering 23:083501, 2015
@ -65,22 +70,13 @@ class Rotation:
def __repr__(self): def __repr__(self):
"""Represent rotation as unit quaternion, rotation matrix, and Bunge-Euler angles.""" """Represent rotation as unit quaternion(s)."""
if self == Rotation(): return f'Quaternion{" " if self.quaternion.shape == (4,) else "s of shape "+str(self.quaternion.shape)+chr(10)}'\
return 'Rotation()' + str(self.quaternion)
else:
return f'Quaternions {self.shape}:\n'+str(self.quaternion) \
if self.quaternion.shape != (4,) else \
'\n'.join([
'Quaternion: (real={:.3f}, imag=<{:+.3f}, {:+.3f}, {:+.3f}>)'.format(*(self.quaternion)),
'Matrix:\n{}'.format(np.round(self.as_matrix(),8)),
'Bunge Eulers / deg: ({:3.2f}, {:3.2f}, {:3.2f})'.format(*self.as_Euler_angles(degrees=True)),
])
# ToDo: Check difference __copy__ vs __deepcopy__
def __copy__(self,**kwargs): def __copy__(self,**kwargs):
"""Copy.""" """Create deep copy."""
return self.__class__(rotation=kwargs['rotation'] if 'rotation' in kwargs else self.quaternion) return self.__class__(rotation=kwargs['rotation'] if 'rotation' in kwargs else self.quaternion)
copy = __copy__ copy = __copy__
@ -97,6 +93,26 @@ class Rotation:
""" """
Equal to other. Equal to other.
Equality is determined taking limited floating point precision into account.
See numpy.allclose for details.
Parameters
----------
other : Rotation
Rotation to check for equality.
"""
s = self.quaternion
o = other.quaternion
if self.shape == () == other.shape:
return np.allclose(s,o) or (np.isclose(s[0],0.0) and np.allclose(s,-1.0*o))
else:
return np.all(np.isclose(s,o),-1) + np.all(np.isclose(s,-1.0*o),-1) * np.isclose(s[...,0],0.0)
def __ne__(self,other):
"""
Not equal to other.
Equality is determined taking limited floating point precision into Equality is determined taking limited floating point precision into
account. See numpy.allclose for details. account. See numpy.allclose for details.
@ -106,8 +122,7 @@ class Rotation:
Rotation to check for equality. Rotation to check for equality.
""" """
return np.prod(self.shape,dtype=int) == np.prod(other.shape,dtype=int) \ return np.logical_not(self==other)
and np.allclose(self.quaternion,other.quaternion)
@property @property
@ -127,41 +142,46 @@ class Rotation:
return dup return dup
def __pow__(self,pwr): def __pow__(self,exp):
""" """
Raise quaternion to power. Perform the rotation 'exp' times.
Equivalent to performing the rotation 'pwr' times.
Parameters Parameters
---------- ----------
pwr : float exp : float
Power to raise quaternion to. Exponent.
""" """
phi = np.arccos(self.quaternion[...,0:1]) phi = np.arccos(self.quaternion[...,0:1])
p = self.quaternion[...,1:]/np.linalg.norm(self.quaternion[...,1:],axis=-1,keepdims=True) p = self.quaternion[...,1:]/np.linalg.norm(self.quaternion[...,1:],axis=-1,keepdims=True)
return self.copy(rotation=Rotation(np.block([np.cos(pwr*phi),np.sin(pwr*phi)*p]))._standardize()) return self.copy(rotation=Rotation(np.block([np.cos(exp*phi),np.sin(exp*phi)*p]))._standardize())
def __ipow__(self,exp):
def __mul__(self,other):
"""Standard multiplication is not implemented."""
raise NotImplementedError('Use "R@b", i.e. matmul, to apply rotation "R" to object "b"')
def __matmul__(self,other):
""" """
Rotation of vector, second or fourth order tensor, or rotation object. Perform the rotation 'exp' times (in-place).
Parameters Parameters
---------- ----------
other : numpy.ndarray or Rotation exp : float
Vector, second or fourth order tensor, or rotation object that is rotated. Exponent.
"""
return self**exp
def __mul__(self,other):
"""
Compose this rotation with other.
Parameters
----------
other : Rotation of shape(self.shape)
Rotation for composition.
Returns Returns
------- -------
other_rot : numpy.ndarray or Rotation composition : Rotation
Rotated vector, second or fourth order tensor, or rotation object. Compound rotation self*other, i.e. first other then self rotation.
""" """
if isinstance(other,Rotation): if isinstance(other,Rotation):
@ -172,8 +192,71 @@ class Rotation:
q = (q_m*q_o - np.einsum('...i,...i',p_m,p_o).reshape(self.shape+(1,))) q = (q_m*q_o - np.einsum('...i,...i',p_m,p_o).reshape(self.shape+(1,)))
p = q_m*p_o + q_o*p_m + _P * np.cross(p_m,p_o) p = q_m*p_o + q_o*p_m + _P * np.cross(p_m,p_o)
return Rotation(np.block([q,p]))._standardize() return Rotation(np.block([q,p]))._standardize()
else:
raise TypeError('Use "R@b", i.e. matmul, to apply rotation "R" to object "b"')
elif isinstance(other,np.ndarray): def __imul__(self,other):
"""
Compose this rotation with other (in-place).
Parameters
----------
other : Rotation of shape(self.shape)
Rotation for composition.
"""
return self*other
def __truediv__(self,other):
"""
Compose this rotation with inverse of other.
Parameters
----------
other : damask.Rotation of shape (self.shape)
Rotation to inverse composition.
Returns
-------
composition : Rotation
Compound rotation self*(~other), i.e. first inverse of other then self rotation.
"""
if isinstance(other,Rotation):
return self*~other
else:
raise TypeError('Use "R@b", i.e. matmul, to apply rotation "R" to object "b"')
def __itruediv__(self,other):
"""
Compose this rotation with inverse of other (in-place).
Parameters
----------
other : Rotation of shape (self.shape)
Rotation to inverse composition.
"""
return self/other
def __matmul__(self,other):
"""
Rotation of vector, second order tensor, or fourth order tensor.
Parameters
----------
other : numpy.ndarray of shape (...,3), (...,3,3), or (...,3,3,3,3)
Vector or tensor on which to apply the rotation.
Returns
-------
rotated : numpy.ndarray of shape (...,3), (...,3,3), or (...,3,3,3,3)
Rotated vector or tensor, i.e. transformed to frame defined by rotation.
"""
if isinstance(other,np.ndarray):
if self.shape + (3,) == other.shape: if self.shape + (3,) == other.shape:
q_m = self.quaternion[...,0] q_m = self.quaternion[...,0]
p_m = self.quaternion[...,1:] p_m = self.quaternion[...,1:]
@ -193,9 +276,13 @@ class Rotation:
return np.einsum('...im,...jn,...ko,...lp,...mnop',R,R,R,R,other) return np.einsum('...im,...jn,...ko,...lp,...mnop',R,R,R,R,other)
else: else:
raise ValueError('Can only rotate vectors, 2nd order tensors, and 4th order tensors') raise ValueError('Can only rotate vectors, 2nd order tensors, and 4th order tensors')
elif isinstance(other,Rotation):
raise TypeError('Use "R1*R2", i.e. multiplication, to compose rotations "R1" and "R2"')
else: else:
raise TypeError(f'Cannot rotate {type(other)}') raise TypeError(f'Cannot rotate {type(other)}')
apply = __matmul__
def _standardize(self): def _standardize(self):
"""Standardize quaternion (ensure positive real hemisphere).""" """Standardize quaternion (ensure positive real hemisphere)."""
@ -296,7 +383,7 @@ class Rotation:
Rotation to which the misorientation is computed. Rotation to which the misorientation is computed.
""" """
return other@~self return other*~self
################################################################################################ ################################################################################################
@ -819,7 +906,7 @@ class Rotation:
np.sqrt(1-u**2)*np.sin(Theta), np.sqrt(1-u**2)*np.sin(Theta),
u, omega]) u, omega])
return Rotation.from_axis_angle(p) @ center return Rotation.from_axis_angle(p) * center
@staticmethod @staticmethod
@ -870,8 +957,8 @@ class Rotation:
f[::2,:3] *= -1 # flip half the rotation axes to negative sense f[::2,:3] *= -1 # flip half the rotation axes to negative sense
return R_align.broadcast_to(N) \ return R_align.broadcast_to(N) \
@ Rotation.from_axis_angle(p,normalize=True) \ * Rotation.from_axis_angle(p,normalize=True) \
@ Rotation.from_axis_angle(f) * Rotation.from_axis_angle(f)
#################################################################################################### ####################################################################################################
@ -1060,7 +1147,6 @@ class Rotation:
@staticmethod @staticmethod
def _om2ax(om): def _om2ax(om):
"""Rotation matrix to axis angle pair.""" """Rotation matrix to axis angle pair."""
#return Rotation._qu2ax(Rotation._om2qu(om)) # HOTFIX
diag_delta = -_P*np.block([om[...,1,2:3]-om[...,2,1:2], diag_delta = -_P*np.block([om[...,1,2:3]-om[...,2,1:2],
om[...,2,0:1]-om[...,0,2:3], om[...,2,0:1]-om[...,0,2:3],
om[...,0,1:2]-om[...,1,0:1] om[...,0,1:2]-om[...,1,0:1]

View File

@ -42,12 +42,10 @@ class Table:
return len(self.data) return len(self.data)
def __copy__(self): def __copy__(self):
"""Copy Table.""" """Create deep copy."""
return copy.deepcopy(self) return copy.deepcopy(self)
def copy(self): copy = __copy__
"""Copy Table."""
return self.__copy__()
def _label_discrete(self): def _label_discrete(self):

View File

@ -22,6 +22,10 @@ class TestConfig:
with open(tmp_path/'config.yaml') as f: with open(tmp_path/'config.yaml') as f:
assert Config.load(f) == config assert Config.load(f) == config
def test_add_remove(self):
config = Config()
assert config.add({'hello':'world'}).delete('hello') == config
def test_repr(self,tmp_path): def test_repr(self,tmp_path):
config = Config() config = Config()
config['A'] = 1 config['A'] = 1

View File

@ -25,13 +25,16 @@ class TestOrientation:
@pytest.mark.parametrize('shape',[None,5,(4,6)]) @pytest.mark.parametrize('shape',[None,5,(4,6)])
def test_equal(self,lattice,shape): def test_equal(self,lattice,shape):
R = Rotation.from_random(shape) R = Rotation.from_random(shape)
assert Orientation(R,lattice) == Orientation(R,lattice) assert Orientation(R,lattice) == Orientation(R,lattice) if shape is None else \
(Orientation(R,lattice) == Orientation(R,lattice)).all()
@pytest.mark.parametrize('lattice',Orientation.crystal_families) @pytest.mark.parametrize('lattice',Orientation.crystal_families)
@pytest.mark.parametrize('shape',[None,5,(4,6)]) @pytest.mark.parametrize('shape',[None,5,(4,6)])
def test_unequal(self,lattice,shape): def test_unequal(self,lattice,shape):
R = Rotation.from_random(shape) R = Rotation.from_random(shape)
assert not(Orientation(R,lattice) != Orientation(R,lattice)) assert not ( Orientation(R,lattice) != Orientation(R,lattice) if shape is None else \
(Orientation(R,lattice) != Orientation(R,lattice)).any())
@pytest.mark.parametrize('a,b',[ @pytest.mark.parametrize('a,b',[
(dict(rotation=[1,0,0,0]), (dict(rotation=[1,0,0,0]),
@ -403,7 +406,7 @@ class TestOrientation:
def test_relationship_vectorize(self,set_of_quaternions,lattice,model): def test_relationship_vectorize(self,set_of_quaternions,lattice,model):
r = Orientation(rotation=set_of_quaternions[:200].reshape((50,4,4)),lattice=lattice).related(model) r = Orientation(rotation=set_of_quaternions[:200].reshape((50,4,4)),lattice=lattice).related(model)
for i in range(200): for i in range(200):
assert r.reshape((-1,200))[:,i] == Orientation(set_of_quaternions[i],lattice).related(model) assert (r.reshape((-1,200))[:,i] == Orientation(set_of_quaternions[i],lattice).related(model)).all()
@pytest.mark.parametrize('model',['Bain','KS','GT','GT_prime','NW','Pitsch']) @pytest.mark.parametrize('model',['Bain','KS','GT','GT_prime','NW','Pitsch'])
@pytest.mark.parametrize('lattice',['cF','cI']) @pytest.mark.parametrize('lattice',['cF','cI'])

View File

@ -526,7 +526,7 @@ class TestRotation:
o = backward(forward(m)) o = backward(forward(m))
u = np.array([np.pi*2,np.pi,np.pi*2]) u = np.array([np.pi*2,np.pi,np.pi*2])
ok = np.allclose(m,o,atol=atol) ok = np.allclose(m,o,atol=atol)
ok = ok or np.allclose(np.where(np.isclose(m,u),m-u,m),np.where(np.isclose(o,u),o-u,o),atol=atol) ok |= np.allclose(np.where(np.isclose(m,u),m-u,m),np.where(np.isclose(o,u),o-u,o),atol=atol)
if np.isclose(m[1],0.0,atol=atol) or np.isclose(m[1],np.pi,atol=atol): if np.isclose(m[1],0.0,atol=atol) or np.isclose(m[1],np.pi,atol=atol):
sum_phi = np.unwrap([m[0]+m[2],o[0]+o[2]]) sum_phi = np.unwrap([m[0]+m[2],o[0]+o[2]])
ok |= np.isclose(sum_phi[0],sum_phi[1],atol=atol) ok |= np.isclose(sum_phi[0],sum_phi[1],atol=atol)
@ -550,19 +550,22 @@ class TestRotation:
assert ok and np.isclose(np.linalg.norm(o[:3]),1.0) and o[3]<=np.pi+1.e-9, f'{m},{o},{rot.as_quaternion()}' assert ok and np.isclose(np.linalg.norm(o[:3]),1.0) and o[3]<=np.pi+1.e-9, f'{m},{o},{rot.as_quaternion()}'
@pytest.mark.parametrize('forward,backward',[(Rotation._ro2qu,Rotation._qu2ro), @pytest.mark.parametrize('forward,backward',[(Rotation._ro2qu,Rotation._qu2ro),
#(Rotation._ro2om,Rotation._om2ro), (Rotation._ro2om,Rotation._om2ro),
#(Rotation._ro2eu,Rotation._eu2ro), (Rotation._ro2eu,Rotation._eu2ro),
(Rotation._ro2ax,Rotation._ax2ro), (Rotation._ro2ax,Rotation._ax2ro),
(Rotation._ro2ho,Rotation._ho2ro), (Rotation._ro2ho,Rotation._ho2ro),
(Rotation._ro2cu,Rotation._cu2ro)]) (Rotation._ro2cu,Rotation._cu2ro)])
def test_Rodrigues_internal(self,set_of_rotations,forward,backward): def test_Rodrigues_internal(self,set_of_rotations,forward,backward):
"""Ensure invariance of conversion from Rodrigues-Frank vector and back.""" """Ensure invariance of conversion from Rodrigues-Frank vector and back."""
cutoff = np.tan(np.pi*.5*(1.-1e-4)) cutoff = np.tan(np.pi*.5*(1.-1e-5))
for rot in set_of_rotations: for rot in set_of_rotations:
m = rot.as_Rodrigues_vector() m = rot.as_Rodrigues_vector()
o = backward(forward(m)) o = backward(forward(m))
ok = np.allclose(np.clip(m,None,cutoff),np.clip(o,None,cutoff),atol=atol) ok = np.allclose(np.clip(m,None,cutoff),np.clip(o,None,cutoff),atol=atol)
ok = ok or np.isclose(m[3],0.0,atol=atol) ok |= np.isclose(m[3],0.0,atol=atol)
if m[3] > cutoff:
ok |= np.allclose(m[:3],-1*o[:3])
assert ok and np.isclose(np.linalg.norm(o[:3]),1.0), f'{m},{o},{rot.as_quaternion()}' assert ok and np.isclose(np.linalg.norm(o[:3]),1.0), f'{m},{o},{rot.as_quaternion()}'
@pytest.mark.parametrize('forward,backward',[(Rotation._ho2qu,Rotation._qu2ho), @pytest.mark.parametrize('forward,backward',[(Rotation._ho2qu,Rotation._qu2ho),
@ -592,7 +595,7 @@ class TestRotation:
o = backward(forward(m)) o = backward(forward(m))
ok = np.allclose(m,o,atol=atol) ok = np.allclose(m,o,atol=atol)
if np.count_nonzero(np.isclose(np.abs(o),np.pi**(2./3.)*.5)): if np.count_nonzero(np.isclose(np.abs(o),np.pi**(2./3.)*.5)):
ok = ok or np.allclose(m*-1.,o,atol=atol) ok |= np.allclose(m*-1.,o,atol=atol)
assert ok and np.max(np.abs(o)) < np.pi**(2./3.) * 0.5 + 1.e-9, f'{m},{o},{rot.as_quaternion()}' assert ok and np.max(np.abs(o)) < np.pi**(2./3.) * 0.5 + 1.e-9, f'{m},{o},{rot.as_quaternion()}'
@pytest.mark.parametrize('vectorized, single',[(Rotation._qu2om,qu2om), @pytest.mark.parametrize('vectorized, single',[(Rotation._qu2om,qu2om),
@ -719,7 +722,7 @@ class TestRotation:
o = Rotation.from_axis_angle(rot.as_axis_angle()).as_axis_angle() o = Rotation.from_axis_angle(rot.as_axis_angle()).as_axis_angle()
ok = np.allclose(m,o,atol=atol) ok = np.allclose(m,o,atol=atol)
if np.isclose(m[3],np.pi,atol=atol): if np.isclose(m[3],np.pi,atol=atol):
ok = ok or np.allclose(m*np.array([-1.,-1.,-1.,1.]),o,atol=atol) ok |= np.allclose(m*np.array([-1.,-1.,-1.,1.]),o,atol=atol)
assert ok and np.isclose(np.linalg.norm(o[:3]),1.0) \ assert ok and np.isclose(np.linalg.norm(o[:3]),1.0) \
and o[3]<=np.pi+1.e-9, f'{m},{o},{rot.as_quaternion()}' and o[3]<=np.pi+1.e-9, f'{m},{o},{rot.as_quaternion()}'
@ -740,7 +743,7 @@ class TestRotation:
m = rot.as_Rodrigues_vector() m = rot.as_Rodrigues_vector()
o = Rotation.from_homochoric(rot.as_homochoric()*P*-1,P).as_Rodrigues_vector() o = Rotation.from_homochoric(rot.as_homochoric()*P*-1,P).as_Rodrigues_vector()
ok = np.allclose(np.clip(m,None,cutoff),np.clip(o,None,cutoff),atol=atol) ok = np.allclose(np.clip(m,None,cutoff),np.clip(o,None,cutoff),atol=atol)
ok = ok or np.isclose(m[3],0.0,atol=atol) ok |= np.isclose(m[3],0.0,atol=atol)
assert ok and np.isclose(np.linalg.norm(o[:3]),1.0), f'{m},{o},{rot.as_quaternion()}' assert ok and np.isclose(np.linalg.norm(o[:3]),1.0), f'{m},{o},{rot.as_quaternion()}'
@pytest.mark.parametrize('P',[1,-1]) @pytest.mark.parametrize('P',[1,-1])
@ -780,8 +783,22 @@ class TestRotation:
else: else:
assert r.shape == shape assert r.shape == shape
def test_equal(self): @pytest.mark.parametrize('shape',[None,5,(4,6)])
assert Rotation.from_random(rng_seed=1) == Rotation.from_random(rng_seed=1) def test_equal(self,shape):
R = Rotation.from_random(shape,rng_seed=1)
assert R == R if shape is None else (R == R).all()
@pytest.mark.parametrize('shape',[None,5,(4,6)])
def test_unequal(self,shape):
R = Rotation.from_random(shape,rng_seed=1)
assert not (R != R if shape is None else (R != R).any())
def test_equal_ambiguous(self):
qu = np.random.rand(10,4)
qu[:,0] = 0.
qu/=np.linalg.norm(qu,axis=1,keepdims=True)
assert (Rotation(qu) == Rotation(-qu)).all()
def test_inversion(self): def test_inversion(self):
r = Rotation.from_random() r = Rotation.from_random()
@ -798,7 +815,7 @@ class TestRotation:
p = Rotation.from_random(shape=shape) p = Rotation.from_random(shape=shape)
s = r.append(p) s = r.append(p)
print(f'append 2x {shape} --> {s.shape}') print(f'append 2x {shape} --> {s.shape}')
assert s[0,...] == r[0,...] and s[-1,...] == p[-1,...] assert np.logical_and(s[0,...] == r[0,...], s[-1,...] == p[-1,...]).all()
@pytest.mark.parametrize('shape',[None,1,(1,),(4,2),(3,3,2)]) @pytest.mark.parametrize('shape',[None,1,(1,),(4,2),(3,3,2)])
def test_append_list(self,shape): def test_append_list(self,shape):
@ -806,7 +823,7 @@ class TestRotation:
p = Rotation.from_random(shape=shape) p = Rotation.from_random(shape=shape)
s = r.append([r,p]) s = r.append([r,p])
print(f'append 3x {shape} --> {s.shape}') print(f'append 3x {shape} --> {s.shape}')
assert s[0,...] == r[0,...] and s[-1,...] == p[-1,...] assert np.logical_and(s[0,...] == r[0,...], s[-1,...] == p[-1,...]).all()
@pytest.mark.parametrize('quat,standardized',[ @pytest.mark.parametrize('quat,standardized',[
([-1,0,0,0],[1,0,0,0]), ([-1,0,0,0],[1,0,0,0]),
@ -828,7 +845,7 @@ class TestRotation:
@pytest.mark.parametrize('order',['C','F']) @pytest.mark.parametrize('order',['C','F'])
def test_flatten_reshape(self,shape,order): def test_flatten_reshape(self,shape,order):
r = Rotation.from_random(shape=shape) r = Rotation.from_random(shape=shape)
assert r == r.flatten(order).reshape(shape,order) assert (r == r.flatten(order).reshape(shape,order)).all()
@pytest.mark.parametrize('function',[Rotation.from_quaternion, @pytest.mark.parametrize('function',[Rotation.from_quaternion,
Rotation.from_Euler_angles, Rotation.from_Euler_angles,
@ -939,7 +956,7 @@ class TestRotation:
def test_rotate_inverse(self): def test_rotate_inverse(self):
R = Rotation.from_random() R = Rotation.from_random()
assert np.allclose(np.eye(3),(~R@R).as_matrix()) assert np.allclose(np.eye(3),(~R*R).as_matrix())
@pytest.mark.parametrize('data',[np.random.rand(3), @pytest.mark.parametrize('data',[np.random.rand(3),
np.random.rand(3,3), np.random.rand(3,3),
@ -973,6 +990,42 @@ class TestRotation:
R_2 = Rotation.from_Euler_angles([360,0,0],degrees=True) R_2 = Rotation.from_Euler_angles([360,0,0],degrees=True)
assert np.allclose(R_1.misorientation(R_2).as_matrix(),np.eye(3)) assert np.allclose(R_1.misorientation(R_2).as_matrix(),np.eye(3))
def test_composition(self):
a,b = (Rotation.from_random(),Rotation.from_random())
c = a * b
a *= b
assert c == a
def test_composition_invalid(self):
with pytest.raises(TypeError):
Rotation()*np.ones(3)
def test_composition_inverse(self):
a,b = (Rotation.from_random(),Rotation.from_random())
c = a / b
a /= b
assert c == a
def test_composition_inverse_invalid(self):
with pytest.raises(TypeError):
Rotation()/np.ones(3)
def test_power(self):
a = Rotation.from_random()
r = (np.random.rand()-.5)*4
b = a**r
a **= r
assert a == b
def test_invariant(self):
R = Rotation.from_random()
assert R/R == R*R**(-1) == Rotation()
@pytest.mark.parametrize('item',[np.ones(3),np.ones((3,3)), np.ones((3,3,3,3))])
def test_apply(self,item):
r = Rotation.from_random()
assert (r.apply(item) == r@item).all()
@pytest.mark.parametrize('angle',[10,20,30,40,50,60,70,80,90,100,120]) @pytest.mark.parametrize('angle',[10,20,30,40,50,60,70,80,90,100,120])
def test_average(self,angle): def test_average(self,angle):
R = Rotation.from_axis_angle([[0,0,1,10],[0,0,1,angle]],degrees=True) R = Rotation.from_axis_angle([[0,0,1,10],[0,0,1,angle]],degrees=True)