From 5f3f87cd68b26081b10ad8a96c655b3a56348830 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 3 May 2020 18:51:30 +0200 Subject: [PATCH] vectorizing cubochoric conversions --- python/damask/_rotation.py | 104 ++++++++++++++++++++++------------ python/tests/test_Rotation.py | 91 +++++++++++++++++------------ 2 files changed, 124 insertions(+), 71 deletions(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 878face97..a13a9af0b 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -241,7 +241,7 @@ class Rotation: """Homochoric vector: (h_1, h_2, h_3).""" return Rotation.qu2ho(self.quaternion) - def asCubochoric(self): + def as_cubochoric(self): """Cubochoric vector: (c_1, c_2, c_3).""" return Rotation.qu2cu(self.quaternion) @@ -265,6 +265,7 @@ class Rotation: asMatrix = as_matrix asRodrigues = as_Rodrigues asHomochoric = as_homochoric + asCubochoric = as_cubochoric ################################################################################################ # Static constructors. The input data needs to follow the conventions, options allow to @@ -386,7 +387,7 @@ class Rotation: return Rotation(Rotation.ho2qu(ho)) @staticmethod - def fromCubochoric(cubochoric, + def from_cubochoric(cubochoric, P = -1): cu = np.array(cubochoric,dtype=float) @@ -461,6 +462,7 @@ class Rotation: fromMatrix = from_matrix fromRodrigues = from_Rodrigues fromHomochoric = from_homochoric + fromCubochoric = from_cubochoric fromRandom = from_random #################################################################################################### @@ -1066,7 +1068,7 @@ class Rotation: if np.allclose(ho,0.0,rtol=0.0,atol=1.0e-16): cu = np.zeros(3) else: - xyz3 = ho[Rotation._get_order(ho,'forward')] + xyz3 = ho[Rotation._get_pyramid_order(ho,'forward')] # inverse M_3 xyz2 = xyz3[0:2] * np.sqrt( 2.0*rs/(rs+np.abs(xyz3[2])) ) @@ -1087,20 +1089,35 @@ class Rotation: # inverse M_1 cu = np.array([ Tinv[0], Tinv[1], (-1.0 if xyz3[2] < 0.0 else 1.0) * rs / np.sqrt(6.0/np.pi) ]) /_sc - # reverse the coordinates back to the regular order according to the original pyramid number - cu = cu[Rotation._get_order(ho,'backward')] - - return cu + cu = cu[Rotation._get_pyramid_order(ho,'backward')] else: - q2 = qxy + np.max(np.abs(xyz2))**2 - sq2 = np.sqrt(q2) - q = (_beta/np.sqrt(2.0)/_R1) * np.sqrt(q2*qxy/(q2-np.max(np.abs(xyz2))*sq2)) - tt = np.clip((np.min(np.abs(xyz2))**2+np.max(np.abs(xyz2))*sq2)/np.sqrt(2.0)/qxy,-1.0,1.0) - Tinv = np.array([1.0,np.arccos(tt)/np.pi*12.0]) if np.abs(xyz2[1]) <= np.abs(xyz2[0]) else \ - np.array([np.arccos(tt)/np.pi*12.0,1.0]) - Tinv = q * np.where(xyz2<0.0,-Tinv,Tinv) - raise NotImplementedError('Support for multiple rotations missing') + rs = np.linalg.norm(ho,axis=-1,keepdims=True) + xyz3 = np.take_along_axis(ho,Rotation._get_pyramid_order(ho,'forward'),-1) + + with np.errstate(invalid='ignore',divide='ignore'): + # inverse M_3 + xyz2 = xyz3[...,0:2] * np.sqrt( 2.0*rs/(rs+np.abs(xyz3[...,2:3])) ) + qxy = np.sum(xyz2**2,axis=-1,keepdims=True) + + q2 = qxy + np.max(np.abs(xyz2),axis=-1,keepdims=True)**2 + sq2 = np.sqrt(q2) + q = (_beta/np.sqrt(2.0)/_R1) * np.sqrt(q2*qxy/(q2-np.max(np.abs(xyz2),axis=-1,keepdims=True)*sq2)) + tt = np.clip((np.min(np.abs(xyz2),axis=-1,keepdims=True)**2\ + +np.max(np.abs(xyz2),axis=-1,keepdims=True)*sq2)/np.sqrt(2.0)/qxy,-1.0,1.0) + T_inv = np.where(np.abs(xyz2[...,1:2]) <= np.abs(xyz2[...,0:1]), + np.block([np.ones_like(tt),np.arccos(tt)/np.pi*12.0]), + np.block([np.arccos(tt)/np.pi*12.0,np.ones_like(tt)]))*q + T_inv[xyz2<0.0] *= -1.0 + T_inv[np.broadcast_to(np.isclose(qxy,0.0,rtol=0.0,atol=1.0e-12),T_inv.shape)] = 0.0 + cu = np.block([T_inv, np.where(xyz3[...,2:3]<0.0,-np.ones_like(xyz3[...,2:3]),np.ones_like(xyz3[...,2:3])) \ + * rs/np.sqrt(6.0/np.pi), + ])/ _sc + + cu[np.isclose(np.sum(np.abs(ho),axis=-1),0.0,rtol=0.0,atol=1.0e-16)] = 0.0 + cu = np.take_along_axis(cu,Rotation._get_pyramid_order(ho,'backward'),-1) + + return cu #---------- Cubochoric ---------- @staticmethod @@ -1145,7 +1162,7 @@ class Rotation: ho = np.zeros(3) else: # get pyramide and scale by grid parameter ratio - XYZ = cu[Rotation._get_order(cu,'forward')] * _sc + XYZ = cu[Rotation._get_pyramid_order(cu,'forward')] * _sc # intercept all the points along the z-axis if np.allclose(XYZ[0:2],0.0,rtol=0.0,atol=1.0e-16): @@ -1163,31 +1180,46 @@ class Rotation: c = np.sum(T**2) s = c * np.pi/24.0 /XYZ[2]**2 c = c * np.sqrt(np.pi/24.0)/XYZ[2] + q = np.sqrt( 1.0 - s ) ho = np.array([ T[order[1]] * q, T[order[0]] * q, np.sqrt(6.0/np.pi) * XYZ[2] - c ]) - # reverse the coordinates back to the regular order according to the original pyramid number - ho = ho[Rotation._get_order(cu,'backward')] - - return ho + ho = ho[Rotation._get_pyramid_order(cu,'backward')] else: - # get pyramide and scale by grid parameter ratio - XYZ = cu[Rotation._get_order(cu,'forward')] * _sc - order = np.where(np.abs(XYZ[...,1]) <= np.abs(XYZ[...,0]), - np.broadcast_to([1,0],XYZ.shape[:-1]+(2,)), - np.broadcast_to([0,1],XYZ.shape[:-1]+(2,))) - q = np.pi/12.0 * XYZ[order[...,0]]/XYZ[order[...,1]] - c = np.cos(q) - s = np.sin(q) - q = _R1*2.0**0.25/_beta * XYZ[order[...,1]] / np.sqrt(np.sqrt(2.0)-c) - T = np.block([ (np.sqrt(2.0)*c - 1.0), np.sqrt(2.0) * s]) * q - s = c * np.pi/24.0 /XYZ[...,2]**2 - c = c * np.sqrt(np.pi/24.0)/XYZ[...,2] - q = np.sqrt( 1.0 - s ) - raise NotImplementedError('Support for multiple rotations missing') + with np.errstate(invalid='ignore',divide='ignore'): + # get pyramide and scale by grid parameter ratio + XYZ = np.take_along_axis(cu,Rotation._get_pyramid_order(cu,'forward'),-1) * _sc + order = np.abs(XYZ[...,1:2]) <= np.abs(XYZ[...,0:1]) + q = np.pi/12.0 * np.where(order,XYZ[...,1:2],XYZ[...,0:1]) \ + / np.where(order,XYZ[...,0:1],XYZ[...,1:2]) + c = np.cos(q) + s = np.sin(q) + q = _R1*2.0**0.25/_beta/ np.sqrt(np.sqrt(2.0)-c) \ + * np.where(order,XYZ[...,0:1],XYZ[...,1:2]) + + T = np.block([ (np.sqrt(2.0)*c - 1.0), np.sqrt(2.0) * s]) * q + + # transform to sphere grid (inverse Lambert) + c = np.sum(T**2,axis=-1,keepdims=True) + s = c * np.pi/24.0 /XYZ[...,2:3]**2 + c = c * np.sqrt(np.pi/24.0)/XYZ[...,2:3] + q = np.sqrt( 1.0 - s) + + ho = np.where(np.isclose(np.sum(np.abs(XYZ[...,0:2]),axis=-1,keepdims=True),0.0,rtol=0.0,atol=1.0e-16), + np.block([np.zeros_like(XYZ[...,0:2]),np.sqrt(6.0/np.pi) *XYZ[...,2:3]]), + np.block([np.where(order,T[...,0:1],T[...,1:2])*q, + np.where(order,T[...,1:2],T[...,0:1])*q, + np.sqrt(6.0/np.pi) * XYZ[...,2:3] - c]) + ) + + ho[np.isclose(np.sum(np.abs(cu),axis=-1),0.0,rtol=0.0,atol=1.0e-16)] = 0.0 + ho = np.take_along_axis(ho,Rotation._get_pyramid_order(cu,'backward'),-1) + + return ho + @staticmethod - def _get_order(xyz,direction=None): + def _get_pyramid_order(xyz,direction=None): """ Get order of the coordinates. @@ -1206,7 +1238,7 @@ class Rotation: """ order = {'forward':np.array([[0,1,2],[1,2,0],[2,0,1]]), - 'backward':np.array([[0,1,2],[2,0,1],[1,2,0]])} + 'backward':np.array([[0,1,2],[2,0,1],[1,2,0]])} if len(xyz.shape) == 1: if np.maximum(abs(xyz[0]),abs(xyz[1])) <= xyz[2] or \ np.maximum(abs(xyz[0]),abs(xyz[1])) <=-xyz[2]: diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py index a46cf1a23..b7442035f 100644 --- a/python/tests/test_Rotation.py +++ b/python/tests/test_Rotation.py @@ -78,9 +78,9 @@ def default(): specials_scatter /= np.linalg.norm(specials_scatter,axis=1).reshape(-1,1) specials_scatter[specials_scatter[:,0]<0]*=-1 - return [Rotation.fromQuaternion(s) for s in specials] + \ - [Rotation.fromQuaternion(s) for s in specials_scatter] + \ - [Rotation.fromRandom() for _ in range(n-len(specials)-len(specials_scatter))] + return [Rotation.from_quaternion(s) for s in specials] + \ + [Rotation.from_quaternion(s) for s in specials_scatter] + \ + [Rotation.from_random() for _ in range(n-len(specials)-len(specials_scatter))] @pytest.fixture def reference_dir(reference_dir_base): @@ -92,41 +92,41 @@ class TestRotation: def test_Eulers(self,default): for rot in default: - m = rot.asQuaternion() - o = Rotation.fromEulers(rot.asEulers()).asQuaternion() + m = rot.as_quaternion() + o = Rotation.from_Eulers(rot.as_Eulers()).as_quaternion() ok = np.allclose(m,o,atol=atol) - if np.isclose(rot.asQuaternion()[0],0.0,atol=atol): + if np.isclose(rot.as_quaternion()[0],0.0,atol=atol): ok = ok or np.allclose(m*-1.,o,atol=atol) - print(m,o,rot.asQuaternion()) + print(m,o,rot.as_quaternion()) assert ok and np.isclose(np.linalg.norm(o),1.0) def test_AxisAngle(self,default): for rot in default: - m = rot.asEulers() - o = Rotation.fromAxisAngle(rot.asAxisAngle()).asEulers() + m = rot.as_Eulers() + o = Rotation.from_axis_angle(rot.as_axis_angle()).as_Eulers() u = np.array([np.pi*2,np.pi,np.pi*2]) 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) 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]]) ok = ok or np.isclose(sum_phi[0],sum_phi[1],atol=atol) - print(m,o,rot.asQuaternion()) + print(m,o,rot.as_quaternion()) assert ok and (np.zeros(3)-1.e-9 <= o).all() and (o <= np.array([np.pi*2.,np.pi,np.pi*2.])+1.e-9).all() def test_Matrix(self,default): for rot in default: - m = rot.asAxisAngle() - o = Rotation.fromAxisAngle(rot.asAxisAngle()).asAxisAngle() + m = rot.as_axis_angle() + o = Rotation.from_axis_angle(rot.as_axis_angle()).as_axis_angle() ok = np.allclose(m,o,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) - print(m,o,rot.asQuaternion()) + print(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 def test_Rodrigues(self,default): for rot in default: - m = rot.asMatrix() - o = Rotation.fromRodrigues(rot.asRodrigues()).asMatrix() + m = rot.as_matrix() + o = Rotation.from_Rodrigues(rot.as_Rodrigues()).as_matrix() ok = np.allclose(m,o,atol=atol) print(m,o) assert ok and np.isclose(np.linalg.det(o),1.0) @@ -134,27 +134,27 @@ class TestRotation: def test_Homochoric(self,default): cutoff = np.tan(np.pi*.5*(1.-1e-4)) for rot in default: - m = rot.asRodrigues() - o = Rotation.fromHomochoric(rot.asHomochoric()).asRodrigues() + m = rot.as_Rodrigues() + o = Rotation.from_homochoric(rot.as_homochoric()).as_Rodrigues() 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) - print(m,o,rot.asQuaternion()) + print(m,o,rot.as_quaternion()) assert ok and np.isclose(np.linalg.norm(o[:3]),1.0) def test_Cubochoric(self,default): for rot in default: - m = rot.asHomochoric() - o = Rotation.fromCubochoric(rot.asCubochoric()).asHomochoric() + m = rot.as_homochoric() + o = Rotation.from_cubochoric(rot.as_cubochoric()).as_homochoric() ok = np.allclose(m,o,atol=atol) - print(m,o,rot.asQuaternion()) + print(m,o,rot.as_quaternion()) assert ok and np.linalg.norm(o) < (3.*np.pi/4.)**(1./3.) + 1.e-9 def test_Quaternion(self,default): for rot in default: - m = rot.asCubochoric() - o = Rotation.fromQuaternion(rot.asQuaternion()).asCubochoric() + m = rot.as_cubochoric() + o = Rotation.from_quaternion(rot.as_quaternion()).as_cubochoric() ok = np.allclose(m,o,atol=atol) - print(m,o,rot.asQuaternion()) + print(m,o,rot.as_quaternion()) assert ok and o.max() < np.pi**(2./3.)*0.5+1.e-9 @pytest.mark.parametrize('function',[Rotation.from_quaternion, @@ -185,9 +185,11 @@ class TestRotation: Rotation.qu2eu, Rotation.qu2ax, Rotation.qu2ro, - Rotation.qu2ho]) + Rotation.qu2ho, + Rotation.qu2cu + ]) def test_quaternion_vectorization(self,default,conversion): - qu = np.array([rot.asQuaternion() for rot in default]) + qu = np.array([rot.as_quaternion() for rot in default]) conversion(qu.reshape(qu.shape[0]//2,-1,4)) co = conversion(qu) for q,c in zip(qu,co): @@ -199,9 +201,10 @@ class TestRotation: Rotation.om2ax, Rotation.om2ro, Rotation.om2ho, + Rotation.om2cu ]) def test_matrix_vectorization(self,default,conversion): - om = np.array([rot.asMatrix() for rot in default]) + om = np.array([rot.as_matrix() for rot in default]) conversion(om.reshape(om.shape[0]//2,-1,3,3)) co = conversion(om) for o,c in zip(om,co): @@ -213,9 +216,10 @@ class TestRotation: Rotation.eu2ax, Rotation.eu2ro, Rotation.eu2ho, + Rotation.eu2cu ]) def test_Euler_vectorization(self,default,conversion): - eu = np.array([rot.asEulers() for rot in default]) + eu = np.array([rot.as_Eulers() for rot in default]) conversion(eu.reshape(eu.shape[0]//2,-1,3)) co = conversion(eu) for e,c in zip(eu,co): @@ -227,9 +231,10 @@ class TestRotation: Rotation.ax2eu, Rotation.ax2ro, Rotation.ax2ho, + Rotation.ax2cu ]) def test_axisAngle_vectorization(self,default,conversion): - ax = np.array([rot.asAxisAngle() for rot in default]) + ax = np.array([rot.as_axis_angle() for rot in default]) conversion(ax.reshape(ax.shape[0]//2,-1,4)) co = conversion(ax) for a,c in zip(ax,co): @@ -242,9 +247,10 @@ class TestRotation: Rotation.ro2eu, Rotation.ro2ax, Rotation.ro2ho, + Rotation.ro2cu ]) def test_Rodrigues_vectorization(self,default,conversion): - ro = np.array([rot.asRodrigues() for rot in default]) + ro = np.array([rot.as_Rodrigues() for rot in default]) conversion(ro.reshape(ro.shape[0]//2,-1,4)) co = conversion(ro) for r,c in zip(ro,co): @@ -256,26 +262,41 @@ class TestRotation: Rotation.ho2eu, Rotation.ho2ax, Rotation.ho2ro, + Rotation.ho2cu ]) def test_homochoric_vectorization(self,default,conversion): - ho = np.array([rot.asHomochoric() for rot in default]) + ho = np.array([rot.as_homochoric() for rot in default]) conversion(ho.reshape(ho.shape[0]//2,-1,3)) co = conversion(ho) for h,c in zip(ho,co): print(h,c) assert np.allclose(conversion(h),c) + @pytest.mark.parametrize('conversion',[Rotation.cu2qu, + Rotation.cu2om, + Rotation.cu2eu, + Rotation.cu2ax, + Rotation.cu2ro, + Rotation.cu2ho + ]) + def test_cubochoric_vectorization(self,default,conversion): + cu = np.array([rot.as_cubochoric() for rot in default]) + conversion(cu.reshape(cu.shape[0]//2,-1,3)) + co = conversion(cu) + for u,c in zip(cu,co): + print(u,c) + assert np.allclose(conversion(u),c) @pytest.mark.parametrize('direction',['forward', 'backward']) def test_pyramid_vectorization(self,direction): p = np.random.rand(n,3) - o = Rotation._get_order(p,direction) + o = Rotation._get_pyramid_order(p,direction) for i,o_i in enumerate(o): - assert np.allclose(o_i,Rotation._get_order(p[i],direction)) + assert np.all(o_i==Rotation._get_pyramid_order(p[i],direction)) def test_pyramid_invariant(self): a = np.random.rand(n,3) - f = damask.Rotation._get_order(a,'forward') - b = damask.Rotation._get_order(a,'backward') + f = Rotation._get_pyramid_order(a,'forward') + b = Rotation._get_pyramid_order(a,'backward') assert np.all(np.take_along_axis(np.take_along_axis(a,f,-1),b,-1) == a)