From 9c7d69dc5c11d23f2db550a9d441a82c4b6a56f9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Sun, 21 Jul 2019 23:37:25 +0300 Subject: [PATCH] Get rid of dependency on Moble's quaternions. Former-commit-id: 182aeda66992aa43c2c172c28a199feba0d6878f --- TODO.md | 2 ++ qpms/qpms_c.pyx | 12 ++++++++++++ qpms/qpms_cdefs.pxd | 1 + qpms/tmatrices.py | 34 +++++++++++++++++++++++++++------- qpms/wigner.h | 7 +++++++ setup.py | 6 ++++-- 6 files changed, 53 insertions(+), 9 deletions(-) diff --git a/TODO.md b/TODO.md index d9792ce..445f787 100644 --- a/TODO.md +++ b/TODO.md @@ -7,6 +7,8 @@ TODO list before public release - Field calculations. - Complex frequencies, n's, k's. - Transforming point (meta)generators. +- Check whether moble's quaternions and my + quaternions give the same results in tmatrices.py - Ewald summations of all types of lattices (dimensionality-wise). - Split lattices.h into separate point generator and lattice vector manipulation parts. * Maybe move something from the .h to .c file. diff --git a/qpms/qpms_c.pyx b/qpms/qpms_c.pyx index fe8d5fd..eddf328 100644 --- a/qpms/qpms_c.pyx +++ b/qpms/qpms_c.pyx @@ -919,6 +919,18 @@ cdef class CQuat: if (abs(m) > l or abs(mp) > l): return 0 return qpms_wignerD_elem(self.q, l, mp, m) + + @staticmethod + def from_rotvector(vec): + if vec.shape != (3,): + raise ValueError("Single 3d vector expected") + res = CQuat() + cdef cart3_t v + v.x = vec[0] + v.y = vec[1] + v.z = vec[2] + res.q = qpms_quat_from_rotvector(v) + return res cdef class IRot3: ''' diff --git a/qpms/qpms_cdefs.pxd b/qpms/qpms_cdefs.pxd index 1bcbf69..1f06284 100644 --- a/qpms/qpms_cdefs.pxd +++ b/qpms/qpms_cdefs.pxd @@ -188,6 +188,7 @@ cdef extern from "wigner.h": qpms_m_t mp, qpms_m_t m) qpms_irot3_t qpms_irot3_mult(qpms_irot3_t p, qpms_irot3_t q) qpms_irot3_t qpms_irot3_pow(qpms_irot3_t p, int n) + qpms_quat_t qpms_quat_from_rotvector(cart3_t v) cdef extern from "groups.h": struct qpms_finite_group_irrep_t: diff --git a/qpms/tmatrices.py b/qpms/tmatrices.py index 02203c3..fac40cf 100644 --- a/qpms/tmatrices.py +++ b/qpms/tmatrices.py @@ -1,9 +1,15 @@ import numpy as np -import quaternion, spherical_functions as sf # because of the Wigner matrices. These imports are SLOW. +use_moble_quaternion = False +try: + import quaternion, spherical_functions as sf # because of the Wigner matrices. These imports are SLOW. + use_moble_quaternion = True +except ImportError: + use_moble_quaternion = False + import re from scipy import interpolate from scipy.constants import hbar, e as eV, pi, c -from qpms_c import get_mn_y, get_nelem +from qpms_c import get_mn_y, get_nelem, CQuat ň = np.newaxis from .types import NormalizationT, TMatrixSpec @@ -17,15 +23,26 @@ def WignerD_mm(l, quat): TODO doc """ - indices = np.array([ [l,i,j] for i in range(-l,l+1) for j in range(-l,l+1)]) - Delems = sf.Wigner_D_element(quat, indices).reshape(2*l+1,2*l+1) - return Delems + if use_moble_quaternion: + indices = np.array([ [l,i,j] for i in range(-l,l+1) for j in range(-l,l+1)]) + Delems = sf.Wigner_D_element(quat, indices).reshape(2*l+1,2*l+1) + return Delems + else: + Delems = np.zeros((2*l+1, 2*l+1), dtype=complex) + for i in range(-l,l+1): + for j in range(-l,l+1): + Delems[i,j] = quat.wignerDelem(l, i, j) + return Delems + def WignerD_mm_fromvector(l, vect): """ TODO doc """ - return WignerD_mm(l, quaternion.from_rotation_vector(vect)) + if use_moble_quaternion: + return WignerD_mm(l, quaternion.from_rotation_vector(vect)) + else: + return WignerD_mm(l, CQuat.from_rotvector(vect)) def WignerD_yy(lmax, quat): @@ -46,7 +63,10 @@ def WignerD_yy_fromvector(lmax, vect): """ TODO doc """ - return WignerD_yy(lmax, quaternion.from_rotation_vector(vect)) + if use_moble_quaternion: + return WignerD_yy(lmax, quaternion.from_rotation_vector(vect)) + else: + return WignerD_yy(lMax, CQuat.from_rotvector(vect)) def identity_yy(lmax): """ diff --git a/qpms/wigner.h b/qpms/wigner.h index 6fd8ad1..5328b56 100644 --- a/qpms/wigner.h +++ b/qpms/wigner.h @@ -113,6 +113,7 @@ static inline qpms_quat_t qpms_quat_pow(const qpms_quat_t q, const double expone return qpms_quat_exp(qe); } + /// Quaternion inversion. /** \f[ q^{-1} = \frac{q*}{|q|}. \f] */ static inline qpms_quat_t qpms_quat_inv(const qpms_quat_t q) { @@ -144,6 +145,12 @@ static inline cart3_t qpms_quat_rot_cart3(qpms_quat_t q, const cart3_t v) { qpms_quat_mult(vv, qc))); } +/// Versor quaternion from rotation vector (norm of the vector is the rotation angle). +static inline qpms_quat_t qpms_quat_from_rotvector(cart3_t v) { + return qpms_quat_exp(qpms_quat_rscale(0.5, + qpms_quat_from_cart3(v))); +} + /// Wigner D matrix element from a rotator quaternion for integer \a l. /** * The D matrix are calculated using formulae (3), (4), (6), (7) from diff --git a/setup.py b/setup.py index f6c78cd..8578bec 100755 --- a/setup.py +++ b/setup.py @@ -98,8 +98,10 @@ setup(name='qpms', packages=['qpms'], # libraries = [('amos', {'sources': amos_sources} )], setup_requires=['cython>=0.28',], - install_requires=['cython>=0.28','quaternion','spherical_functions','scipy>=0.18.0', 'sympy>=1.2'], - dependency_links=['https://github.com/moble/quaternion/archive/v2.0.tar.gz','https://github.com/moble/spherical_functions/archive/master.zip'], + install_requires=['cython>=0.28', + #'quaternion','spherical_functions', + 'scipy>=0.18.0', 'sympy>=1.2'], + #dependency_links=['https://github.com/moble/quaternion/archive/v2.0.tar.gz','https://github.com/moble/spherical_functions/archive/master.zip'], ext_modules=cythonize([qpms_c], include_path=['qpms', 'amos'], gdb_debug=True), cmdclass = {'build_ext': build_ext}, )