From 0a34e342e4e1f6815eae75922325e3ec1e1736ca Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Tue, 15 Sep 2020 18:14:15 -0400 Subject: [PATCH] fixed fiber component sampling and testing --- python/damask/_rotation.py | 40 +++++++++++++++-------------------- python/tests/test_Rotation.py | 21 +++++++++--------- 2 files changed, 27 insertions(+), 34 deletions(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 8950ee38c..ab22c8a6f 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -687,55 +687,49 @@ class Rotation: @staticmethod - def from_fiber_component(alpha,beta,FWHM=0.0,N=500,degrees=True,seed=None): + def from_fiber_component(alpha,beta,sigma=0.0,N=500,degrees=True,seed=None): """ Calculate set of rotations with Gaussian distribution around direction. - References - ---------- - K. Helming, Texturapproximation durch Modellkomponenten - Cuvillier Verlag, 1996 - Parameters ---------- alpha : numpy.ndarray of size 2 - tbd. + Polar coordinates (phi from x,theta from z) of fiber direction in crystal frame. beta : numpy.ndarray of size 2 - tbd. - FWHM : float, optional - Full width at half maximum of the Gaussian distribution. + Polar coordinates (phi from x,theta from z) of fiber direction in sample frame. + sigma : float, optional + Standard deviation of (Gaussian) misorientation distribution. Defaults to 0. - N_samples : int, optional + N : int, optional Number of samples, defaults to 500. degrees : boolean, optional - FWHM, alpha, and beta are given in degrees. + sigma, alpha, and beta are given in degrees. seed : {None, int, array_like[ints], SeedSequence, BitGenerator, Generator}, optional - A seed to initialize the BitGenerator. Defaults to None. - If None, then fresh, unpredictable entropy will be pulled from the OS. + A seed to initialize the BitGenerator. Defaults to None, i.e. unpredictable entropy + will be pulled from the OS. """ rng = np.random.default_rng(seed) - FWHM_,alpha_,beta_ = map(np.radians,(FWHM,alpha,beta)) if degrees else (FWHM,alpha,beta) + sigma_,alpha_,beta_ = map(np.radians,(sigma,alpha,beta)) if degrees else (sigma,alpha,beta) d_cr = np.array([np.sin(alpha_[0])*np.cos(alpha_[1]), np.sin(alpha_[0])*np.sin(alpha_[1]), np.cos(alpha_[0])]) - d_lab = np.array([np.sin(beta_[0] )*np.cos(beta_[1] ), np.sin(beta_[0] )*np.sin(beta_[1] ), np.cos(beta_[0] )]) + d_lab = np.array([np.sin( beta_[0])*np.cos( beta_[1]), np.sin( beta_[0])*np.sin( beta_[1]), np.cos( beta_[0])]) ax_align = np.append(np.cross(d_lab,d_cr), np.arccos(np.dot(d_lab,d_cr))) 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) # rotation to align fiber axis in crystal and sample system u,v,b = (np.random.random((N,3)) * 2 * np.array([1,np.pi,np.pi]) - np.array([1,0,np.pi])).T - a = abs(np.random.normal(scale=FWHM_,size=N)) + a = abs(np.random.normal(scale=sigma_,size=N)) p = np.vstack((np.sqrt(1-u**2)*np.cos(v), np.sqrt(1-u**2)*np.sin(v), u, a)).T - - f = np.hstack((np.broadcast_to(d_cr,(N,3)),b.reshape(N,1))) + p[:,:3] = np.einsum('ij,...j->...i',np.eye(3)-np.outer(d_lab,d_lab),p[:,:3]) + f = np.hstack((np.broadcast_to(d_lab,(N,3)),b.reshape(N,1))) f[f[:,3]<0] *= -1. - return Rotation.from_axis_angle(p) \ - * Rotation.from_axis_angle(f) \ - * R_align - + return R_align.broadcast_to(N) \ + @ Rotation.from_axis_angle(p,normalize=True) \ + @ Rotation.from_axis_angle(f) #################################################################################################### diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py index 92566b42c..1ddc03ef6 100644 --- a/python/tests/test_Rotation.py +++ b/python/tests/test_Rotation.py @@ -925,14 +925,14 @@ class TestRotation: dist = angles * (np.random.randint(0,2,N_samples)*2-1) p = stats.normaltest(dist)[1] - FWHM_out = np.std(dist) * (2*np.sqrt(2*np.log(2))) + FWHM_out = np.std(dist) # * (2*np.sqrt(2*np.log(2))) print(f'\np: {p}, FWHM ratio {FWHM/FWHM_out}') assert (.9 < FWHM/FWHM_out < 1.1) and p > 0.001 - @pytest.mark.parametrize('FWHM',[10,15,20]) - @pytest.mark.parametrize('N_samples',[500,1000,2000]) - def test_from_fiber_component(self,N_samples,FWHM): + @pytest.mark.parametrize('sigma',[5,10,15,20]) + @pytest.mark.parametrize('N',[1000,10000,100000]) + def test_from_fiber_component(self,N,sigma): """https://en.wikipedia.org/wiki/Full_width_at_half_maximum.""" alpha = np.random.random(2)*np.pi beta = np.random.random(2)*np.pi @@ -942,12 +942,11 @@ 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(FWHM),N_samples,False) - angles=[np.arccos(np.clip(np.dot(n@f_in_S,o[i]@f_in_S),-1,1)) for i in range(N_samples)] - dist = np.array(angles) * (np.random.randint(0,2,N_samples)*2-1) + 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) p = stats.normaltest(dist)[1] - FWHM_out = np.degrees(np.std(dist)) * (2*np.sqrt(2*np.log(2))) - print(f'\np: {p}, FWHM ratio {FWHM/FWHM_out}') - assert (.85 < FWHM/FWHM_out < 1.15) and p > 0.001 - + sigma_out = np.degrees(np.std(dist)) + print(f'\np: {p}, sigma ratio {sigma/sigma_out}') + assert (.85 < sigma/sigma_out < 1.15) and p > 0.001