223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

Merge branch 'consistent-orientation-from' into 'development' 

consistent "shape" keyword in from_X

Closes #165

See merge request damask/DAMASK!546
This commit is contained in:
Philip Eisenlohr 2022-03-20 00:00:25 +00:00
parent b5d6f2dc83 52d595ff62
commit 5b87fafcae
4 changed files with 56 additions and 60 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
PRIVATE

@ -1 +1 @@
Subproject commit 459326e9840c843ade72b04cf28e50889c9779f1
Subproject commit 00b3eb79ee6f8df2ca50276d2111008e5f79b3e1

75
python/damask/_rotation.py
View File

@ -418,15 +418,15 @@ class Rotation:
def reshape(self: MyType,
shape: Union[int, Tuple[int, ...]],
shape: Union[int, IntSequence],
order: Literal['C','F','A'] = 'C') -> MyType:
"""
Reshape array.
Parameters
----------
shape : int or tuple of ints
The new shape should be compatible with the original shape.
shape : int or sequence of ints
New shape, number of elements needs to match the original shape.
If an integer is supplied, then the result will be a 1-D array of that length.
order : {'C', 'F', 'A'}, optional
'C' flattens in row-major (C-style) order.
@ -446,15 +446,15 @@ class Rotation:
def broadcast_to(self: MyType,
shape: Union[int, Tuple[int, ...]],
shape: Union[int, IntSequence],
mode: Literal['left', 'right'] = 'right') -> MyType:
"""
Broadcast array.
Parameters
----------
shape : int or tuple of ints
Shape of broadcasted array.
shape : int or sequence of ints
Shape of broadcasted array, needs to be compatible with the original shape.
mode : str, optional
Where to preferentially locate missing dimensions.
Either 'left' or 'right' (default).
@ -465,9 +465,9 @@ class Rotation:
Rotation broadcasted to given shape.
"""
if isinstance(shape,(int,np.integer)): shape = (shape,)
return self.copy(np.broadcast_to(self.quaternion.reshape(util.shapeshifter(self.shape,shape,mode)+(4,)),
shape+(4,)))
shape_ = (shape,) if isinstance(shape,(int,np.integer)) else tuple(shape)
return self.copy(np.broadcast_to(self.quaternion.reshape(util.shapeshifter(self.shape,shape_,mode)+(4,)),
shape_+(4,)))
def average(self: MyType,
@ -979,17 +979,15 @@ class Rotation:
@staticmethod
def from_random(shape: Tuple[int, ...] = None,
def from_random(shape: Union[int, IntSequence] = None,
rng_seed: NumpyRngSeed = None) -> 'Rotation':
"""
Initialize with random rotation.
Rotations are uniformly distributed.
Initialize with samples from a uniform distribution.
Parameters
----------
shape : tuple of ints, optional
Shape of the sample. Defaults to None, which gives a single rotation.
shape : int or sequence of ints, optional
Shape of the returned array. Defaults to None, which gives a scalar.
rng_seed : {None, int, array_like[ints], SeedSequence, BitGenerator, Generator}, optional
A seed to initialize the BitGenerator.
Defaults to None, i.e. unpredictable entropy will be pulled from the OS.
@ -1011,12 +1009,12 @@ class Rotation:
@staticmethod
def from_ODF(weights: np.ndarray,
phi: np.ndarray,
N: int = 500,
shape: Union[int, IntSequence] = None,
degrees: bool = True,
fractions: bool = True,
rng_seed: NumpyRngSeed = None) -> 'Rotation':
"""
Sample discrete values from a binned orientation distribution function (ODF).
Initialize with samples from a binned orientation distribution function (ODF).
Parameters
----------
@ -1024,9 +1022,8 @@ class Rotation:
Texture intensity values (probability density or volume fraction) at Euler space grid points.
phi : numpy.ndarray, shape (n,3)
Grid coordinates in Euler space at which weights are defined.
N : integer, optional
Number of discrete orientations to be sampled from the given ODF.
Defaults to 500.
shape : int or sequence of ints, optional
Shape of the returned array. Defaults to None, which gives a scalar.
degrees : bool, optional
Euler space grid coordinates are in degrees. Defaults to True.
fractions : bool, optional
@ -1036,11 +1033,6 @@ class Rotation:
A seed to initialize the BitGenerator.
Defaults to None, i.e. unpredictable entropy will be pulled from the OS.
Returns
-------
samples : damask.Rotation, shape (N)
Array of sampled rotations that approximate the input ODF.
Notes
-----
Due to the distortion of Euler space in the vicinity of ϕ = 0, probability densities, p, defined on
@ -1063,26 +1055,27 @@ class Rotation:
dg = 1.0 if fractions else _dg(phi,degrees)
dV_V = dg * np.maximum(0.0,weights.squeeze())
return Rotation.from_Euler_angles(phi[util.hybrid_IA(dV_V,N,rng_seed)],degrees)
N = 1 if shape is None else np.prod(shape)
return Rotation.from_Euler_angles(phi[util.hybrid_IA(dV_V,N,rng_seed)],degrees).reshape(() if shape is None else shape)
@staticmethod
def from_spherical_component(center: 'Rotation',
sigma: float,
N: int = 500,
shape: Union[int, IntSequence] = None,
degrees: bool = True,
rng_seed: NumpyRngSeed = None) -> 'Rotation':
"""
Calculate set of rotations with Gaussian distribution around center.
Initialize with samples from a Gaussian distribution around a given center.
Parameters
----------
center : Rotation
Central Rotation.
center : Rotation or Orientation
Central rotation.
sigma : float
Standard deviation of (Gaussian) misorientation distribution.
N : int, optional
Number of samples. Defaults to 500.
shape : int or sequence of ints, optional
Shape of the returned array. Defaults to None, which gives a scalar.
degrees : bool, optional
sigma is given in degrees. Defaults to True.
rng_seed : {None, int, array_like[ints], SeedSequence, BitGenerator, Generator}, optional
@ -1092,24 +1085,25 @@ class Rotation:
"""
rng = np.random.default_rng(rng_seed)
sigma = np.radians(sigma) if degrees else sigma
N = 1 if shape is None else np.prod(shape)
u,Theta = (rng.random((N,2)) * 2.0 * np.array([1,np.pi]) - np.array([1.0, 0])).T
omega = abs(rng.normal(scale=sigma,size=N))
p = np.column_stack([np.sqrt(1-u**2)*np.cos(Theta),
np.sqrt(1-u**2)*np.sin(Theta),
u, omega])
return Rotation.from_axis_angle(p) * center
return Rotation.from_axis_angle(p).reshape(() if shape is None else shape) * center
@staticmethod
def from_fiber_component(alpha: IntSequence,
beta: IntSequence,
sigma: float = 0.0,
N: int = 500,
shape: Union[int, IntSequence] = None,
degrees: bool = True,
rng_seed: NumpyRngSeed = None):
"""
Calculate set of rotations with Gaussian distribution around direction.
Initialize with samples from a Gaussian distribution around a given direction.
Parameters
----------
@ -1120,8 +1114,8 @@ class Rotation:
sigma : float, optional
Standard deviation of (Gaussian) misorientation distribution.
Defaults to 0.
N : int, optional
Number of samples. Defaults to 500.
shape : int or sequence of ints, optional
Shape of the returned array. Defaults to None, which gives a scalar.
degrees : bool, optional
sigma, alpha, and beta are given in degrees.
rng_seed : {None, int, array_like[ints], SeedSequence, BitGenerator, Generator}, optional
@ -1139,6 +1133,7 @@ class Rotation:
if np.isclose(ax_align[3],0.0): ax_align[:3] = np.array([1,0,0])
R_align = Rotation.from_axis_angle(ax_align if ax_align[3] > 0.0 else -ax_align,normalize=True) # rotate fiber axis from sample to crystal frame
N = 1 if shape is None else np.prod(shape)
u,Theta = (rng.random((N,2)) * 2.0 * np.array([1,np.pi]) - np.array([1.0, 0])).T
omega = abs(rng.normal(scale=sigma_,size=N))
p = np.column_stack([np.sqrt(1-u**2)*np.cos(Theta),
@ -1148,9 +1143,9 @@ class Rotation:
f = np.column_stack((np.broadcast_to(d_lab,(N,3)),rng.random(N)*np.pi))
f[::2,:3] *= -1 # flip half the rotation axes to negative sense
return R_align.broadcast_to(N) \
* Rotation.from_axis_angle(p,normalize=True) \
* Rotation.from_axis_angle(f)
return (R_align.broadcast_to(N)
* Rotation.from_axis_angle(p,normalize=True)
* Rotation.from_axis_angle(f)).reshape(() if shape is None else shape)
####################################################################################################

6
python/tests/test_Orientation.py
View File

@ -146,14 +146,14 @@ class TestOrientation:
def test_from_spherical_component(self):
assert np.all(Orientation.from_spherical_component(center=Rotation(),
sigma=0.0,N=1,family='triclinic').as_matrix()
sigma=0.0,shape=1,family='triclinic').as_matrix()
== np.eye(3))
def test_from_fiber_component(self):
r = Rotation.from_fiber_component(alpha=np.zeros(2),beta=np.zeros(2),
sigma=0.0,N=1,rng_seed=0)
sigma=0.0,shape=1,rng_seed=0)
assert np.all(Orientation.from_fiber_component(alpha=np.zeros(2),beta=np.zeros(2),
sigma=0.0,N=1,rng_seed=0,family='triclinic').quaternion
sigma=0.0,shape=None,rng_seed=0,family='triclinic').quaternion
== r.quaternion)
@pytest.mark.parametrize('kwargs',[

33
python/tests/test_Rotation.py
View File

@ -6,6 +6,7 @@ from damask import Rotation
from damask import Table
from damask import _rotation
from damask import grid_filters
from damask import util
n = 1000
atol=1.e-4
@ -1055,12 +1056,12 @@ class TestRotation:
@pytest.mark.parametrize('sigma',[5,10,15,20])
@pytest.mark.parametrize('N',[1000,10000,100000])
def test_spherical_component(self,N,sigma):
@pytest.mark.parametrize('shape',[1000,10000,100000,(10,100)])
def test_spherical_component(self,sigma,shape):
p = []
for run in range(5):
c = Rotation.from_random()
o = Rotation.from_spherical_component(c,sigma,N)
o = Rotation.from_spherical_component(c,sigma,shape)
_, angles = c.misorientation(o).as_axis_angle(pair=True,degrees=True)
angles[::2] *= -1 # flip angle for every second to symmetrize distribution
@ -1072,8 +1073,8 @@ class TestRotation:
@pytest.mark.parametrize('sigma',[5,10,15,20])
@pytest.mark.parametrize('N',[1000,10000,100000])
def test_from_fiber_component(self,N,sigma):
@pytest.mark.parametrize('shape',[1000,10000,100000,(10,100)])
def test_from_fiber_component(self,sigma,shape):
p = []
for run in range(5):
alpha = np.random.random()*2*np.pi,np.arccos(np.random.random())
@ -1084,9 +1085,9 @@ class TestRotation:
ax = np.append(np.cross(f_in_C,f_in_S), - np.arccos(np.dot(f_in_C,f_in_S)))
n = Rotation.from_axis_angle(ax if ax[3] > 0.0 else ax*-1.0 ,normalize=True) # rotation to align fiber axis in crystal and sample system
o = Rotation.from_fiber_component(alpha,beta,np.radians(sigma),N,False)
angles = np.arccos(np.clip(np.dot(o@np.broadcast_to(f_in_S,(N,3)),n@f_in_S),-1,1))
dist = np.array(angles) * (np.random.randint(0,2,N)*2-1)
o = Rotation.from_fiber_component(alpha,beta,np.radians(sigma),shape,False)
angles = np.arccos(np.clip(np.dot(o@np.broadcast_to(f_in_S,tuple(util.aslist(shape))+(3,)),n@f_in_S),-1,1))
dist = np.array(angles) * (np.random.randint(0,2,util.aslist(shape))*2-1)
p.append(stats.normaltest(dist)[1])
@ -1097,8 +1098,8 @@ class TestRotation:
@pytest.mark.parametrize('fractions',[True,False])
@pytest.mark.parametrize('degrees',[True,False])
@pytest.mark.parametrize('N',[2**13,2**14,2**15])
def test_ODF_cell(self,ref_path,fractions,degrees,N):
@pytest.mark.parametrize('shape',[2**13,2**14,2**15,(2**8,2**6)])
def test_ODF_cell(self,ref_path,fractions,degrees,shape):
steps = np.array([144,36,36])
limits = np.array([360.,90.,90.])
rng = tuple(zip(np.zeros(3),limits))
@ -1107,14 +1108,14 @@ class TestRotation:
Eulers = grid_filters.coordinates0_point(steps,limits)
Eulers = np.radians(Eulers) if not degrees else Eulers
Eulers_r = Rotation.from_ODF(weights,Eulers.reshape(-1,3,order='F'),N,degrees,fractions).as_Euler_angles(True)
weights_r = np.histogramdd(Eulers_r,steps,rng)[0].flatten(order='F')/N * np.sum(weights)
Eulers_r = Rotation.from_ODF(weights,Eulers.reshape(-1,3,order='F'),shape,degrees,fractions).as_Euler_angles(True)
weights_r = np.histogramdd(Eulers_r.reshape(-1,3),steps,rng)[0].flatten(order='F')/np.prod(shape) * np.sum(weights)
if fractions: assert np.sqrt(((weights_r - weights) ** 2).mean()) < 4
@pytest.mark.parametrize('degrees',[True,False])
@pytest.mark.parametrize('N',[2**13,2**14,2**15])
def test_ODF_node(self,ref_path,degrees,N):
@pytest.mark.parametrize('shape',[2**13,2**14,2**15,(2**8,2**6)])
def test_ODF_node(self,ref_path,degrees,shape):
steps = np.array([144,36,36])
limits = np.array([360.,90.,90.])
rng = tuple(zip(-limits/steps*.5,limits-limits/steps*.5))
@ -1125,8 +1126,8 @@ class TestRotation:
Eulers = grid_filters.coordinates0_node(steps,limits)[:-1,:-1,:-1]
Eulers = np.radians(Eulers) if not degrees else Eulers
Eulers_r = Rotation.from_ODF(weights,Eulers.reshape(-1,3,order='F'),N,degrees).as_Euler_angles(True)
weights_r = np.histogramdd(Eulers_r,steps,rng)[0].flatten(order='F')/N * np.sum(weights)
Eulers_r = Rotation.from_ODF(weights,Eulers.reshape(-1,3,order='F'),shape,degrees).as_Euler_angles(True)
weights_r = np.histogramdd(Eulers_r.reshape(-1,3),steps,rng)[0].flatten(order='F')/np.prod(shape) * np.sum(weights)
assert np.sqrt(((weights_r - weights) ** 2).mean()) < 5