polishing

This commit is contained in:
Martin Diehl 2020-04-09 14:20:43 +02:00
parent b025c1838e
commit e502573e05
2 changed files with 32 additions and 14 deletions

View File

@ -640,15 +640,25 @@ class Rotation:
i = np.where(np.isclose(w,1.0+0.0j))[0][0] i = np.where(np.isclose(w,1.0+0.0j))[0][0]
ax[0:3] = np.real(vr[0:3,i]) ax[0:3] = np.real(vr[0:3,i])
diagDelta = np.array([om[1,2]-om[2,1],om[2,0]-om[0,2],om[0,1]-om[1,0]]) diagDelta = np.array([om[1,2]-om[2,1],om[2,0]-om[0,2],om[0,1]-om[1,0]])
ax[0:3] = np.where(np.abs(diagDelta)<1.e-6, ax[0:3],np.abs(ax[0:3])*np.sign(-P*diagDelta)) ax[0:3] = np.where(np.abs(diagDelta)<0, ax[0:3],np.abs(ax[0:3])*np.sign(-P*diagDelta))
return ax
else: else:
diag_delta = np.block([om[...,1,2:3]-om[...,2,1:2], diag_delta = np.block([om[...,1,2:3]-om[...,2,1:2],
om[...,2,0:1]-om[...,0,2:3], om[...,2,0:1]-om[...,0,2:3],
om[...,0,1:2]-om[...,1,0:1] om[...,0,1:2]-om[...,1,0:1]
]) ])
t = 0.5*(om.trace(axis2=-2,axis1=-1) -1.0).reshape(om.shape[:-2]+(1,))
w,vr = np.linalg.eig(om) w,vr = np.linalg.eig(om)
# TODO ------------------ # 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.
ax = np.where(np.abs(diag_delta)<0,
np.real(vr[np.isclose(w,1.0+0.0j)]).reshape(om.shape[:-2]+(3,)),
np.real(vr[np.isclose(w,1.0+0.0j)]).reshape(om.shape[:-2]+(3,)) \
* np.abs(diag_delta)*np.sign(-P*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]
return ax
@staticmethod @staticmethod
def om2ro(om): def om2ro(om):

View File

@ -80,7 +80,7 @@ def default():
return [Rotation.fromQuaternion(s) for s in specials] + \ return [Rotation.fromQuaternion(s) for s in specials] + \
[Rotation.fromQuaternion(s) for s in specials_scatter] + \ [Rotation.fromQuaternion(s) for s in specials_scatter] + \
[Rotation.fromRandom() for _ in range(n-len(specials))] [Rotation.fromRandom() for _ in range(n-len(specials)-len(specials_scatter))]
@pytest.fixture @pytest.fixture
def reference_dir(reference_dir_base): def reference_dir(reference_dir_base):
@ -159,17 +159,22 @@ class TestRotation:
Rotation.qu2ho]) Rotation.qu2ho])
def test_quaternion_vectorization(self,default,conversion): def test_quaternion_vectorization(self,default,conversion):
qu = np.array([rot.asQuaternion() for rot in default]) qu = np.array([rot.asQuaternion() for rot in default])
#dev_null = conversion(qu.reshape(qu.shape[0]//2,-1,4))
co = conversion(qu) co = conversion(qu)
for q,c in zip(qu,co): for q,c in zip(qu,co):
assert np.allclose(conversion(q),c) print(q,c)
assert np.allclose(conversion(q),c)
@pytest.mark.parametrize('conversion',[Rotation.om2eu, @pytest.mark.parametrize('conversion',[Rotation.om2eu,
#Rotation.om2ax,
]) ])
def test_matrix_vectorization(self,default,conversion): def test_matrix_vectorization(self,default,conversion):
qu = np.array([rot.asMatrix() for rot in default]) om = np.array([rot.asMatrix() for rot in default])
co = conversion(qu) #dev_null = conversion(om.reshape(om.shape[0]//2,-1,3,3))
for q,c in zip(qu,co): co = conversion(om)
assert np.allclose(conversion(q),c) for o,c in zip(om,co):
print(o,c)
assert np.allclose(conversion(o),c)
@pytest.mark.parametrize('conversion',[Rotation.eu2qu, @pytest.mark.parametrize('conversion',[Rotation.eu2qu,
Rotation.eu2om, Rotation.eu2om,
@ -177,10 +182,12 @@ class TestRotation:
Rotation.eu2ro, Rotation.eu2ro,
]) ])
def test_Euler_vectorization(self,default,conversion): def test_Euler_vectorization(self,default,conversion):
qu = np.array([rot.asEulers() for rot in default]) eu = np.array([rot.asEulers() for rot in default])
co = conversion(qu) #dev_null = conversion(eu.reshape(eu.shape[0]//2,-1,3))
for q,c in zip(qu,co): co = conversion(eu)
assert np.allclose(conversion(q),c) for e,c in zip(eu,co):
print(e,c)
assert np.allclose(conversion(e),c)
@pytest.mark.parametrize('conversion',[Rotation.ax2qu, @pytest.mark.parametrize('conversion',[Rotation.ax2qu,
Rotation.ax2om, Rotation.ax2om,
@ -192,4 +199,5 @@ class TestRotation:
dev_null = conversion(ax.reshape(ax.shape[0]//2,-1,4)) dev_null = conversion(ax.reshape(ax.shape[0]//2,-1,4))
co = conversion(ax) co = conversion(ax)
for a,c in zip(ax,co): for a,c in zip(ax,co):
assert np.allclose(conversion(a),c) print(a,c)
assert np.allclose(conversion(a),c)