From ae3eca5f98a21c745c3dc3ba780517a5848c8bdc Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 21 Apr 2020 11:59:42 +0200 Subject: [PATCH] fix for vectorized from_random --- python/damask/_rotation.py | 45 ++++++++++++++++++++------------------ 1 file changed, 24 insertions(+), 21 deletions(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 8a2441012..69a725e99 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -265,9 +265,9 @@ class Rotation: qu[qu[...,0] < 0.0] *= -1 else: if np.any(qu[...,0] < 0.0): - raise ValueError('Quaternions need to have positive first(real) component.') + raise ValueError('Quaternion with negative first (real) component.') if not np.all(np.isclose(np.linalg.norm(qu,axis=-1), 1.0)): - raise ValueError('Quaternions need to have unit length.') + raise ValueError('Quaternion is not of unit length.') return Rotation(qu) @@ -279,7 +279,7 @@ class Rotation: eu = np.radians(eu) if degrees else eu if np.any(eu < 0.0) or np.any(eu > 2.0*np.pi) or np.any(eu[...,1] > np.pi): # ToDo: No separate check for PHI - raise ValueError('Euler angles need to be in [0..2π],[0..π],[0..2π].') + raise ValueError('Euler angles outside of [0..2π],[0..π],[0..2π].') return Rotation(Rotation.eu2qu(eu)) @@ -315,11 +315,11 @@ class Rotation: (U,S,Vh) = np.linalg.svd(om) # singular value decomposition om = np.einsum('...ij,...jl->...il',U,Vh) if not np.all(np.isclose(np.linalg.det(om),1.0)): - raise ValueError('matrix is not a proper rotation: {}.'.format(om)) + raise ValueError('Orientation matrix has determinant ≠ 1.') if not np.all(np.isclose(np.einsum('...i,...i',om[...,0],om[...,1]), 0.0)) \ or not np.all(np.isclose(np.einsum('...i,...i',om[...,1],om[...,2]), 0.0)) \ or not np.all(np.isclose(np.einsum('...i,...i',om[...,2],om[...,0]), 0.0)): - raise ValueError('matrix is not orthogonal.') + raise ValueError('Orientation matrix is not orthogonal.') return Rotation(Rotation.om2qu(om)) @@ -338,9 +338,9 @@ class Rotation: if P > 0: ro[...,0:3] *= -1 # convert from P=1 to P=-1 if normalise: ro[...,0:3] /= np.linalg.norm(ro[...,0:3],axis=-1) if np.any(ro[...,3] < 0.0): - raise ValueError('Rodrigues rotation angle not positive.') + raise ValueError('Rodrigues vector rotation angle not positive.') if not np.all(np.isclose(np.linalg.norm(ro[...,0:3],axis=-1), 1.0)): - raise ValueError('Rodrigues rotation axis is not of unit length.') + raise ValueError('Rodrigues vector rotation axis is not of unit length.') return Rotation(Rotation.ro2qu(ro)) @@ -353,7 +353,7 @@ class Rotation: if P > 0: ho *= -1 # convert from P=1 to P=-1 if np.any(np.linalg.norm(ho,axis=-1) > (3.*np.pi/4.)**(1./3.)+1e-9): - raise ValueError('Coordinate outside of the sphere.') + raise ValueError('Homochoric coordinate outside of the sphere.') return Rotation(Rotation.ho2qu(ho)) @@ -364,7 +364,7 @@ class Rotation: cu = np.array(cubochoric,dtype=float) if np.abs(np.max(cu))>np.pi**(2./3.) * 0.5+1e-9: - raise ValueError('Coordinate outside of the cube: {} {} {}.'.format(*cu)) + raise ValueError('Cubochoric coordinate outside of the cube: {} {} {}.'.format(*cu)) ho = Rotation.cu2ho(cu) if P > 0: ho *= -1 # convert from P=1 to P=-1 @@ -408,14 +408,20 @@ class Rotation: def from_random(shape=None): if shape is None: r = np.random.random(3) - else: + elif hasattr(shape, '__iter__'): r = np.random.random(tuple(shape)+(3,)) + else: + r = np.random.random((shape,3)) + A = np.sqrt(r[...,2]) B = np.sqrt(1.0-r[...,2]) - return Rotation(np.block([np.cos(2.0*np.pi*r[...,0])*A, - np.sin(2.0*np.pi*r[...,1])*B, - np.cos(2.0*np.pi*r[...,1])*B, - np.sin(2.0*np.pi*r[...,0])*A])).standardize() + q = np.stack([np.cos(2.0*np.pi*r[...,0])*A, + np.sin(2.0*np.pi*r[...,1])*B, + np.cos(2.0*np.pi*r[...,1])*B, + np.sin(2.0*np.pi*r[...,0])*A],axis=-1) + + return Rotation(q.reshape(r.shape[:-1]+(4,)) if shape is not None else q).standardize() + # for compatibility (old names do not follow convention) fromQuaternion = from_quaternion @@ -820,12 +826,11 @@ class Rotation: c = np.cos(ax[3]*0.5) s = np.sin(ax[3]*0.5) qu = np.array([ c, ax[0]*s, ax[1]*s, ax[2]*s ]) - return qu else: c = np.cos(ax[...,3:4]*.5) s = np.sin(ax[...,3:4]*.5) qu = np.where(np.abs(ax[...,3:4])<1.e-6,[1.0, 0.0, 0.0, 0.0],np.block([c, ax[...,:3]*s])) - return qu + return qu @staticmethod def ax2om(ax): @@ -871,7 +876,7 @@ class Rotation: # 180 degree case ro += [np.inf] if np.isclose(ax[3],np.pi,atol=1.0e-15,rtol=0.0) else \ [np.tan(ax[3]*0.5)] - return np.array(ro) + ro = np.array(ro) else: ro = np.block([ax[...,:3], np.where(np.isclose(ax[...,3:4],np.pi,atol=1.e-15,rtol=.0), @@ -879,7 +884,7 @@ class Rotation: np.tan(ax[...,3:4]*0.5)) ]) ro[np.abs(ax[...,3])<1.e-6] = [.0,.0,_P,.0] - return ro + return ro @staticmethod def ax2ho(ax): @@ -887,11 +892,10 @@ class Rotation: if len(ax.shape) == 1: f = (0.75 * ( ax[3] - np.sin(ax[3]) ))**(1.0/3.0) ho = ax[0:3] * f - return ho else: f = (0.75 * ( ax[...,3:4] - np.sin(ax[...,3:4]) ))**(1.0/3.0) ho = ax[...,:3] * f - return ho + return ho @staticmethod def ax2cu(ax): @@ -948,7 +952,6 @@ class Rotation: 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