use consistent "shape" keyword when shape not clear from input data

This commit is contained in:
Philip Eisenlohr 2022-03-14 16:24:05 -04:00
parent a925a99f43
commit 378b8b2396
3 changed files with 41 additions and 35 deletions

View File

@ -1011,7 +1011,7 @@ class Rotation:
@staticmethod @staticmethod
def from_ODF(weights: np.ndarray, def from_ODF(weights: np.ndarray,
phi: np.ndarray, phi: np.ndarray,
N: int = 500, shape: Tuple[int, ...] = None,
degrees: bool = True, degrees: bool = True,
fractions: bool = True, fractions: bool = True,
rng_seed: NumpyRngSeed = None) -> 'Rotation': rng_seed: NumpyRngSeed = None) -> 'Rotation':
@ -1024,9 +1024,9 @@ class Rotation:
Texture intensity values (probability density or volume fraction) at Euler space grid points. Texture intensity values (probability density or volume fraction) at Euler space grid points.
phi : numpy.ndarray, shape (n,3) phi : numpy.ndarray, shape (n,3)
Grid coordinates in Euler space at which weights are defined. Grid coordinates in Euler space at which weights are defined.
N : integer, optional shape : tuple of ints, optional
Number of discrete orientations to be sampled from the given ODF. Shape of the array of discrete rotations sampled from the given ODF.
Defaults to 500. Defaults to None, which gives a single rotation.
degrees : bool, optional degrees : bool, optional
Euler space grid coordinates are in degrees. Defaults to True. Euler space grid coordinates are in degrees. Defaults to True.
fractions : bool, optional fractions : bool, optional
@ -1038,7 +1038,7 @@ class Rotation:
Returns Returns
------- -------
samples : damask.Rotation, shape (N) samples : damask.Rotation
Array of sampled rotations that approximate the input ODF. Array of sampled rotations that approximate the input ODF.
Notes Notes
@ -1063,13 +1063,14 @@ class Rotation:
dg = 1.0 if fractions else _dg(phi,degrees) dg = 1.0 if fractions else _dg(phi,degrees)
dV_V = dg * np.maximum(0.0,weights.squeeze()) 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((-1,) if shape is None else shape)
@staticmethod @staticmethod
def from_spherical_component(center: 'Rotation', def from_spherical_component(center: 'Rotation',
sigma: float, sigma: float,
N: int = 500, shape: Tuple[int, ...] = None,
degrees: bool = True, degrees: bool = True,
rng_seed: NumpyRngSeed = None) -> 'Rotation': rng_seed: NumpyRngSeed = None) -> 'Rotation':
""" """
@ -1081,8 +1082,9 @@ class Rotation:
Central Rotation. Central Rotation.
sigma : float sigma : float
Standard deviation of (Gaussian) misorientation distribution. Standard deviation of (Gaussian) misorientation distribution.
N : int, optional shape : tuple of ints, optional
Number of samples. Defaults to 500. Shape of the array of sampled rotations.
Defaults to None, which gives a single rotation.
degrees : bool, optional degrees : bool, optional
sigma is given in degrees. Defaults to True. sigma is given in degrees. Defaults to True.
rng_seed : {None, int, array_like[ints], SeedSequence, BitGenerator, Generator}, optional rng_seed : {None, int, array_like[ints], SeedSequence, BitGenerator, Generator}, optional
@ -1092,20 +1094,21 @@ class Rotation:
""" """
rng = np.random.default_rng(rng_seed) rng = np.random.default_rng(rng_seed)
sigma = np.radians(sigma) if degrees else sigma 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 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)) omega = abs(rng.normal(scale=sigma,size=N))
p = np.column_stack([np.sqrt(1-u**2)*np.cos(Theta), p = np.column_stack([np.sqrt(1-u**2)*np.cos(Theta),
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).reshape((-1,) if shape is None else shape) * center
@staticmethod @staticmethod
def from_fiber_component(alpha: IntSequence, def from_fiber_component(alpha: IntSequence,
beta: IntSequence, beta: IntSequence,
sigma: float = 0.0, sigma: float = 0.0,
N: int = 500, shape: Tuple[int, ...] = None,
degrees: bool = True, degrees: bool = True,
rng_seed: NumpyRngSeed = None): rng_seed: NumpyRngSeed = None):
""" """
@ -1120,8 +1123,9 @@ class Rotation:
sigma : float, optional sigma : float, optional
Standard deviation of (Gaussian) misorientation distribution. Standard deviation of (Gaussian) misorientation distribution.
Defaults to 0. Defaults to 0.
N : int, optional shape : tuple of ints, optional
Number of samples. Defaults to 500. Shape of the array of sampled rotations.
Defaults to None, which gives a single rotation.
degrees : bool, optional degrees : bool, optional
sigma, alpha, and beta are given in degrees. sigma, alpha, and beta are given in degrees.
rng_seed : {None, int, array_like[ints], SeedSequence, BitGenerator, Generator}, optional rng_seed : {None, int, array_like[ints], SeedSequence, BitGenerator, Generator}, optional
@ -1139,6 +1143,7 @@ class Rotation:
if np.isclose(ax_align[3],0.0): ax_align[:3] = np.array([1,0,0]) 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 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 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)) omega = abs(rng.normal(scale=sigma_,size=N))
p = np.column_stack([np.sqrt(1-u**2)*np.cos(Theta), p = np.column_stack([np.sqrt(1-u**2)*np.cos(Theta),
@ -1148,9 +1153,9 @@ class Rotation:
f = np.column_stack((np.broadcast_to(d_lab,(N,3)),rng.random(N)*np.pi)) 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 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)).reshape((-1,) if shape is None else shape)
#################################################################################################### ####################################################################################################

View File

@ -146,14 +146,14 @@ class TestOrientation:
def test_from_spherical_component(self): def test_from_spherical_component(self):
assert np.all(Orientation.from_spherical_component(center=Rotation(), 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)) == np.eye(3))
def test_from_fiber_component(self): def test_from_fiber_component(self):
r = Rotation.from_fiber_component(alpha=np.zeros(2),beta=np.zeros(2), 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), 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) == r.quaternion)
@pytest.mark.parametrize('kwargs',[ @pytest.mark.parametrize('kwargs',[

View File

@ -6,6 +6,7 @@ from damask import Rotation
from damask import Table from damask import Table
from damask import _rotation from damask import _rotation
from damask import grid_filters from damask import grid_filters
from damask import util
n = 1000 n = 1000
atol=1.e-4 atol=1.e-4
@ -1055,12 +1056,12 @@ class TestRotation:
@pytest.mark.parametrize('sigma',[5,10,15,20]) @pytest.mark.parametrize('sigma',[5,10,15,20])
@pytest.mark.parametrize('N',[1000,10000,100000]) @pytest.mark.parametrize('shape',[1000,10000,100000,(10,100)])
def test_spherical_component(self,N,sigma): def test_spherical_component(self,sigma,shape):
p = [] p = []
for run in range(5): for run in range(5):
c = Rotation.from_random() 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 = c.misorientation(o).as_axis_angle(pair=True,degrees=True)
angles[::2] *= -1 # flip angle for every second to symmetrize distribution 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('sigma',[5,10,15,20])
@pytest.mark.parametrize('N',[1000,10000,100000]) @pytest.mark.parametrize('shape',[1000,10000,100000,(10,100)])
def test_from_fiber_component(self,N,sigma): def test_from_fiber_component(self,sigma,shape):
p = [] p = []
for run in range(5): for run in range(5):
alpha = np.random.random()*2*np.pi,np.arccos(np.random.random()) 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))) 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 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) 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,(N,3)),n@f_in_S),-1,1)) 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,N)*2-1) dist = np.array(angles) * (np.random.randint(0,2,util.aslist(shape))*2-1)
p.append(stats.normaltest(dist)[1]) p.append(stats.normaltest(dist)[1])
@ -1097,8 +1098,8 @@ class TestRotation:
@pytest.mark.parametrize('fractions',[True,False]) @pytest.mark.parametrize('fractions',[True,False])
@pytest.mark.parametrize('degrees',[True,False]) @pytest.mark.parametrize('degrees',[True,False])
@pytest.mark.parametrize('N',[2**13,2**14,2**15]) @pytest.mark.parametrize('shape',[2**13,2**14,2**15,(2**8,2**6)])
def test_ODF_cell(self,ref_path,fractions,degrees,N): def test_ODF_cell(self,ref_path,fractions,degrees,shape):
steps = np.array([144,36,36]) steps = np.array([144,36,36])
limits = np.array([360.,90.,90.]) limits = np.array([360.,90.,90.])
rng = tuple(zip(np.zeros(3),limits)) rng = tuple(zip(np.zeros(3),limits))
@ -1107,14 +1108,14 @@ class TestRotation:
Eulers = grid_filters.coordinates0_point(steps,limits) Eulers = grid_filters.coordinates0_point(steps,limits)
Eulers = np.radians(Eulers) if not degrees else Eulers 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) 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,steps,rng)[0].flatten(order='F')/N * np.sum(weights) 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 if fractions: assert np.sqrt(((weights_r - weights) ** 2).mean()) < 4
@pytest.mark.parametrize('degrees',[True,False]) @pytest.mark.parametrize('degrees',[True,False])
@pytest.mark.parametrize('N',[2**13,2**14,2**15]) @pytest.mark.parametrize('shape',[2**13,2**14,2**15,(2**8,2**6)])
def test_ODF_node(self,ref_path,degrees,N): def test_ODF_node(self,ref_path,degrees,shape):
steps = np.array([144,36,36]) steps = np.array([144,36,36])
limits = np.array([360.,90.,90.]) limits = np.array([360.,90.,90.])
rng = tuple(zip(-limits/steps*.5,limits-limits/steps*.5)) 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 = grid_filters.coordinates0_node(steps,limits)[:-1,:-1,:-1]
Eulers = np.radians(Eulers) if not degrees else Eulers 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) Eulers_r = Rotation.from_ODF(weights,Eulers.reshape(-1,3,order='F'),shape,degrees).as_Euler_angles(True)
weights_r = np.histogramdd(Eulers_r,steps,rng)[0].flatten(order='F')/N * np.sum(weights) 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 assert np.sqrt(((weights_r - weights) ** 2).mean()) < 5