diff --git a/PRIVATE b/PRIVATE index 459326e98..00b3eb79e 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 459326e9840c843ade72b04cf28e50889c9779f1 +Subproject commit 00b3eb79ee6f8df2ca50276d2111008e5f79b3e1 diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 5ceb38cf2..2963925c4 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -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) #################################################################################################### diff --git a/python/tests/test_Orientation.py b/python/tests/test_Orientation.py index bb4a14ed7..4a7b549fd 100644 --- a/python/tests/test_Orientation.py +++ b/python/tests/test_Orientation.py @@ -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',[ diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py index 57f078ae5..a9b8b36f0 100644 --- a/python/tests/test_Rotation.py +++ b/python/tests/test_Rotation.py @@ -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