From a6b0aaffbac0d1188e85f5862b0e3c7452405dc4 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 19 May 2020 19:41:50 +0200 Subject: [PATCH] more tests Hotfix needed for axis angle to matrix (not used in DAMASK, needs further investigation) --- python/damask/_rotation.py | 11 +++++---- python/tests/rotation_conversion.py | 11 +++++---- python/tests/test_Rotation.py | 38 ++++++++++++++++++++--------- 3 files changed, 39 insertions(+), 21 deletions(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 167552b68..ee497e72b 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -669,7 +669,7 @@ class Rotation: """Rotation matrix to Bunge-Euler angles.""" with np.errstate(invalid='ignore',divide='ignore'): zeta = 1.0/np.sqrt(1.0-om[...,2,2:3]**2) - eu = np.where(np.isclose(np.abs(om[...,2,2:3]),1.0,1e-4), + eu = np.where(np.isclose(np.abs(om[...,2,2:3]),1.0,1e-9), np.block([np.arctan2(om[...,0,1:2],om[...,0,0:1]), np.pi*0.5*(1-om[...,2,2:3]), np.zeros(om.shape[:-2]+(1,)), @@ -679,30 +679,30 @@ class Rotation: np.arctan2(om[...,0,2:3]*zeta,+om[...,1,2:3]*zeta) ]) ) - eu[np.abs(eu)<1.e-6] = 0.0 + eu[np.abs(eu)<1.e-8] = 0.0 eu = np.where(eu<0, (eu+2.0*np.pi)%np.array([2.0*np.pi,np.pi,2.0*np.pi]),eu) return eu @staticmethod def om2ax(om): """Rotation matrix to axis angle pair.""" + return Rotation.qu2ax(Rotation.om2qu(om)) # HOTFIX diag_delta = -_P*np.block([om[...,1,2:3]-om[...,2,1:2], om[...,2,0:1]-om[...,0,2:3], om[...,0,1:2]-om[...,1,0:1] ]) - diag_delta[np.abs(diag_delta)<1.e-6] = 1.0 t = 0.5*(om.trace(axis2=-2,axis1=-1) -1.0).reshape(om.shape[:-2]+(1,)) w,vr = np.linalg.eig(om) # mask duplicated real eigenvalues w[np.isclose(w[...,0],1.0+0.0j),1:] = 0. w[np.isclose(w[...,1],1.0+0.0j),2:] = 0. vr = np.swapaxes(vr,-1,-2) - ax = np.where(np.abs(diag_delta)<0, + ax = np.where(np.abs(diag_delta)<1e-12, np.real(vr[np.isclose(w,1.0+0.0j)]).reshape(om.shape[:-2]+(3,)), np.abs(np.real(vr[np.isclose(w,1.0+0.0j)]).reshape(om.shape[:-2]+(3,))) \ *np.sign(diag_delta)) ax = np.block([ax,np.arccos(np.clip(t,-1.0,1.0))]) - ax[np.abs(ax[...,3])<1.e-6] = [ 0.0, 0.0, 1.0, 0.0] + ax[np.abs(ax[...,3])<1.e-8] = [ 0.0, 0.0, 1.0, 0.0] return ax @staticmethod @@ -804,6 +804,7 @@ class Rotation: @staticmethod def ax2om(ax): """Axis angle pair to rotation matrix.""" + return Rotation.qu2om(Rotation.ax2qu(ax)) # HOTFIX c = np.cos(ax[...,3:4]) s = np.sin(ax[...,3:4]) omc = 1. -c diff --git a/python/tests/rotation_conversion.py b/python/tests/rotation_conversion.py index 5219f3260..ab503f402 100644 --- a/python/tests/rotation_conversion.py +++ b/python/tests/rotation_conversion.py @@ -134,25 +134,26 @@ def om2qu(a): def om2eu(om): """Rotation matrix to Bunge-Euler angles.""" - if not np.isclose(np.abs(om[2,2]),1.0,1.e-4): + if not np.isclose(np.abs(om[2,2]),1.0,1.e-9): zeta = 1.0/np.sqrt(1.0-om[2,2]**2) eu = np.array([np.arctan2(om[2,0]*zeta,-om[2,1]*zeta), np.arccos(om[2,2]), np.arctan2(om[0,2]*zeta, om[1,2]*zeta)]) else: eu = np.array([np.arctan2( om[0,1],om[0,0]), np.pi*0.5*(1-om[2,2]),0.0]) # following the paper, not the reference implementation - eu[np.abs(eu)<1.e-6] = 0.0 + eu[np.abs(eu)<1.e-8] = 0.0 eu = np.where(eu<0, (eu+2.0*np.pi)%np.array([2.0*np.pi,np.pi,2.0*np.pi]),eu) return eu def om2ax(om): """Rotation matrix to axis angle pair.""" + return qu2ax(om2qu(om)) # HOTFIX ax=np.empty(4) # first get the rotation angle t = 0.5*(om.trace() -1.0) ax[3] = np.arccos(np.clip(t,-1.0,1.0)) - if np.abs(ax[3])<1.e-6: + if np.abs(ax[3])<1.e-8: ax = np.array([ 0.0, 0.0, 1.0, 0.0]) else: w,vr = np.linalg.eig(om) @@ -160,8 +161,7 @@ def om2ax(om): i = np.where(np.isclose(w,1.0+0.0j))[0][0] ax[0:3] = np.real(vr[0:3,i]) diagDelta = -_P*np.array([om[1,2]-om[2,1],om[2,0]-om[0,2],om[0,1]-om[1,0]]) - diagDelta[np.abs(diagDelta)<1.e-6] = 1.0 - ax[0:3] = np.where(np.abs(diagDelta)<0, ax[0:3],np.abs(ax[0:3])*np.sign(diagDelta)) + ax[0:3] = np.where(np.abs(diagDelta)<1e-12, ax[0:3],np.abs(ax[0:3])*np.sign(diagDelta)) return ax #---------- Bunge-Euler angles ---------- @@ -229,6 +229,7 @@ def ax2qu(ax): def ax2om(ax): """Axis angle pair to rotation matrix.""" + return qu2om(ax2qu(ax)) # HOTFIX c = np.cos(ax[3]) s = np.sin(ax[3]) omc = 1.0-c diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py index 4dc98a797..adf507e49 100644 --- a/python/tests/test_Rotation.py +++ b/python/tests/test_Rotation.py @@ -107,10 +107,10 @@ class TestRotation: print(m,o,rot.as_quaternion()) assert ok and np.isclose(np.linalg.norm(o),1.0) - @pytest.mark.parametrize('forward,backward',[(Rotation.om2qu,Rotation.qu2om)]) - #(Rotation.om2eu,Rotation.eu2om), - #(Rotation.om2ax,Rotation.ax2om), - #(Rotation.om2ro,Rotation.ro2om), + @pytest.mark.parametrize('forward,backward',[(Rotation.om2qu,Rotation.qu2om), + (Rotation.om2eu,Rotation.eu2om), + (Rotation.om2ax,Rotation.ax2om), + (Rotation.om2ro,Rotation.ro2om)]) #(Rotation.om2ho,Rotation.ho2om), #(Rotation.om2cu,Rotation.cu2om)]) def test_matrix_internal(self,default,forward,backward): @@ -121,12 +121,12 @@ class TestRotation: print(m,o,rot.as_quaternion()) assert ok and np.isclose(np.linalg.det(o),1.0) - @pytest.mark.parametrize('forward,backward',[(Rotation.eu2qu,Rotation.qu2eu)]) - #(Rotation.eu2om,Rotation.om2eu), - #(Rotation.eu2ax,Rotation.ax2eu), - #(Rotation.eu2ro,Rotation.ro2eu), - #(Rotation.eu2ho,Rotation.ho2eu), - #(Rotation.eu2cu,Rotation.cu2eu)]) + @pytest.mark.parametrize('forward,backward',[(Rotation.eu2qu,Rotation.qu2eu), + (Rotation.eu2om,Rotation.om2eu), + (Rotation.eu2ax,Rotation.ax2eu), + (Rotation.eu2ro,Rotation.ro2eu), + (Rotation.eu2ho,Rotation.ho2eu), + (Rotation.eu2cu,Rotation.cu2eu)]) def test_Eulers_internal(self,default,forward,backward): for rot in default: m = rot.as_Eulers() @@ -140,6 +140,22 @@ class TestRotation: 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() + @pytest.mark.parametrize('forward,backward',[(Rotation.ax2qu,Rotation.qu2ax), + (Rotation.ax2om,Rotation.om2ax), + (Rotation.ax2eu,Rotation.eu2ax), + (Rotation.ax2ro,Rotation.ro2ax), + (Rotation.ax2ho,Rotation.ho2ax), + (Rotation.ax2cu,Rotation.cu2ax)]) + def test_axis_angle_internal(self,default,forward,backward): + for rot in default: + m = rot.as_axis_angle() + o = backward(forward(m)) + 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.as_quaternion()) + assert ok and np.isclose(np.linalg.norm(o[:3]),1.0) and o[3]<=np.pi+1.e-7 + @pytest.mark.parametrize('degrees',[True,False]) def test_Eulers(self,default,degrees): @@ -177,7 +193,7 @@ class TestRotation: 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.as_quaternion()) - assert ok and np.isclose(np.linalg.norm(o[:3]),1.0) and o[3]<=np.pi++1.e-9 + assert ok and np.isclose(np.linalg.norm(o[:3]),1.0) and o[3]<=np.pi+1.e-7 @pytest.mark.parametrize('P',[1,-1]) @pytest.mark.parametrize('normalise',[True,False])