From 99655c9f61060547ff76cfde2121e324f5587cb9 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 11 Apr 2020 13:57:05 +0200 Subject: [PATCH] more vectorized functions --- python/damask/_rotation.py | 45 ++++++++++++++++++++++------------- python/tests/test_Rotation.py | 32 ++++++++++++++++++++----- 2 files changed, 54 insertions(+), 23 deletions(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index e58de982f..56c5d0593 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -767,8 +767,8 @@ class Rotation: def eu2ro(eu): """Bunge-Euler angles to Rodrigues-Frank vector.""" if len(eu.shape) == 1: - ro = Rotation.eu2ax(eu) # convert to axis angle pair representation - if ro[3] >= np.pi: # Differs from original implementation. check convention 5 + ro = Rotation.eu2ax(eu) # convert to axis angle pair representation + if ro[3] >= np.pi: # Differs from original implementation. check convention 5 ro[3] = np.inf elif iszero(ro[3]): ro = np.array([ 0.0, 0.0, P, 0.0 ]) @@ -902,27 +902,38 @@ class Rotation: @staticmethod def ro2ax(ro): """Rodrigues-Frank vector to axis angle pair.""" - ta = ro[3] - - if iszero(ta): - ax = [ 0.0, 0.0, 1.0, 0.0 ] - elif not np.isfinite(ta): - ax = [ ro[0], ro[1], ro[2], np.pi ] + if len(ro.shape) == 1: + if np.abs(ro[3]) < 1.e-6: + ax = np.array([ 0.0, 0.0, 1.0, 0.0 ]) + elif not np.isfinite(ro[3]): + ax = np.array([ ro[0], ro[1], ro[2], np.pi ]) + else: + angle = 2.0*np.arctan(ro[3]) + ta = np.linalg.norm(ro[0:3]) + ax = np.array([ ro[0]*ta, ro[1]*ta, ro[2]*ta, angle ]) else: - angle = 2.0*np.arctan(ta) - ta = 1.0/np.linalg.norm(ro[0:3]) - ax = [ ro[0]/ta, ro[1]/ta, ro[2]/ta, angle ] - return np.array(ax) + with np.errstate(invalid='ignore',divide='ignore'): + ax = np.where(np.isfinite(ro[...,3:4]), + np.block([ro[...,0:3]*np.linalg.norm(ro[...,0:3],axis=-1,keepdims=True),2.*np.arctan(ro[...,3:4])]), + np.block([ro[...,0:3],np.broadcast_to(np.pi,ro[...,3:4].shape)])) + ax[np.abs(ro[...,3]) < 1.e-6] = np.array([ 0.0, 0.0, 1.0, 0.0 ]) + return ax @staticmethod def ro2ho(ro): """Rodrigues-Frank vector to homochoric vector.""" - if iszero(np.sum(ro[0:3]**2.0)): - ho = [ 0.0, 0.0, 0.0 ] + if len(ro.shape) == 1: + if np.sum(ro[0:3]**2.0) < 1.e-6: + ho = np.zeros(3) + else: + f = 2.0*np.arctan(ro[3]) -np.sin(2.0*np.arctan(ro[3])) if np.isfinite(ro[3]) else np.pi + ho = ro[0:3] * (0.75*f)**(1.0/3.0) else: - f = 2.0*np.arctan(ro[3]) -np.sin(2.0*np.arctan(ro[3])) if np.isfinite(ro[3]) else np.pi - ho = ro[0:3] * (0.75*f)**(1.0/3.0) - return np.array(ho) + f = np.where(np.isfinite(ro[...,3:4]),2.0*np.arctan(ro[...,3:4]) -np.sin(2.0*np.arctan(ro[...,3:4])),np.pi) + ho = np.where(np.broadcast_to(np.sum(ro[...,0:3]**2.0,axis=-1,keepdims=True) < 1.e-6,ro[...,0:3].shape), + np.zeros(3), ro[...,0:3]* (0.75*f)**(1.0/3.0)) + + return ho @staticmethod def ro2cu(ro): diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py index b3d3d881d..544e088c7 100644 --- a/python/tests/test_Rotation.py +++ b/python/tests/test_Rotation.py @@ -123,7 +123,7 @@ class TestRotation: print(m,o,rot.asQuaternion()) assert ok and np.isclose(np.linalg.norm(o[:3]),1.0) and o[3]<=np.pi++1.e-9 - def test_Rodriques(self,default): + def test_Rodrigues(self,default): for rot in default: m = rot.asMatrix() o = Rotation.fromRodrigues(rot.asRodrigues()).asMatrix() @@ -164,18 +164,21 @@ class TestRotation: Rotation.qu2ho]) def test_quaternion_vectorization(self,default,conversion): qu = np.array([rot.asQuaternion() for rot in default]) - dev_null = conversion(qu.reshape(qu.shape[0]//2,-1,4)) + conversion(qu.reshape(qu.shape[0]//2,-1,4)) co = conversion(qu) for q,c in zip(qu,co): print(q,c) assert np.allclose(conversion(q),c) - @pytest.mark.parametrize('conversion',[Rotation.om2eu, + @pytest.mark.parametrize('conversion',[Rotation.om2qu, + Rotation.om2eu, Rotation.om2ax, + Rotation.om2ro, + Rotation.om2ho, ]) def test_matrix_vectorization(self,default,conversion): om = np.array([rot.asMatrix() for rot in default]) - dev_null = conversion(om.reshape(om.shape[0]//2,-1,3,3)) + conversion(om.reshape(om.shape[0]//2,-1,3,3)) co = conversion(om) for o,c in zip(om,co): print(o,c) @@ -185,10 +188,11 @@ class TestRotation: Rotation.eu2om, Rotation.eu2ax, Rotation.eu2ro, + Rotation.eu2ho, ]) def test_Euler_vectorization(self,default,conversion): eu = np.array([rot.asEulers() for rot in default]) - dev_null = conversion(eu.reshape(eu.shape[0]//2,-1,3)) + conversion(eu.reshape(eu.shape[0]//2,-1,3)) co = conversion(eu) for e,c in zip(eu,co): print(e,c) @@ -196,13 +200,29 @@ class TestRotation: @pytest.mark.parametrize('conversion',[Rotation.ax2qu, Rotation.ax2om, + Rotation.ax2eu, Rotation.ax2ro, Rotation.ax2ho, ]) def test_axisAngle_vectorization(self,default,conversion): ax = np.array([rot.asAxisAngle() for rot in default]) - dev_null = conversion(ax.reshape(ax.shape[0]//2,-1,4)) + conversion(ax.reshape(ax.shape[0]//2,-1,4)) co = conversion(ax) for a,c in zip(ax,co): print(a,c) assert np.allclose(conversion(a),c) + + + @pytest.mark.parametrize('conversion',[Rotation.ro2qu, + Rotation.ro2om, + Rotation.ro2eu, + Rotation.ro2ax, + Rotation.ro2ho, + ]) + def test_Rodrigues_vectorization(self,default,conversion): + ro = np.array([rot.asRodrigues() for rot in default]) + conversion(ro.reshape(ro.shape[0]//2,-1,4)) + co = conversion(ro) + for r,c in zip(ro,co): + print(r,c) + assert np.allclose(conversion(r),c)