From 371a8a5f7c9faa890564b77b21269ae659d0bd0a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Fri, 9 Aug 2019 21:54:13 +0300 Subject: [PATCH] Start splitting qpms_c.pyx Former-commit-id: bb2e68dc4cb7f85769ddaf5533298ab1f0e84f5b --- cyquat.pxd | 10 + qpms/cybspec.pxd | 9 + qpms/cybspec.pyx | 121 ++++++++ qpms/cycommon.pyx | 182 ++++++++++++ qpms/cyquaternions.pxd | 8 + qpms/cyquaternions.pyx | 330 +++++++++++++++++++++ qpms/qpms_c.pyx | 630 +---------------------------------------- setup.py | 7 +- 8 files changed, 672 insertions(+), 625 deletions(-) create mode 100644 cyquat.pxd create mode 100644 qpms/cybspec.pxd create mode 100644 qpms/cybspec.pyx create mode 100644 qpms/cycommon.pyx create mode 100644 qpms/cyquaternions.pxd create mode 100644 qpms/cyquaternions.pyx diff --git a/cyquat.pxd b/cyquat.pxd new file mode 100644 index 0000000..f7de57d --- /dev/null +++ b/cyquat.pxd @@ -0,0 +1,10 @@ +from qpms_cdefs cimport * + +cdef class CQuat: + cdef readonly qpms_quat_t q + +cdef class IRot3: + cdef readonly qpms_irot3_t qd + cdef void cset(self, qpms_irot3_t qd) + + diff --git a/qpms/cybspec.pxd b/qpms/cybspec.pxd new file mode 100644 index 0000000..bc5cdf5 --- /dev/null +++ b/qpms/cybspec.pxd @@ -0,0 +1,9 @@ +from qpms_cdefs cimport * + +cimport numpy as np + +cdef class BaseSpec: + cdef qpms_vswf_set_spec_t s + cdef np.ndarray __ilist + + cdef qpms_vswf_set_spec_t *rawpointer(BaseSpec self) diff --git a/qpms/cybspec.pyx b/qpms/cybspec.pyx new file mode 100644 index 0000000..780baf2 --- /dev/null +++ b/qpms/cybspec.pyx @@ -0,0 +1,121 @@ +import numpy as np +import enum +import cycommon + +class VSWFNorm(enum.IntEnum): + # TODO try to make this an enum.IntFlag if supported + # TODO add the other flags from qpms_normalisation_t as well + UNNORM = QPMS_NORMALISATION_NORM_NONE + UNNORM_CS = QPMS_NORMALISATION_NORM_NONE | QPMS_NORMALISATION_CSPHASE + POWERNORM = QPMS_NORMALISATION_NORM_POWER + POWERNORM_CS = QPMS_NORMALISATION_NORM_POWER | QPMS_NORMALISATION_CSPHASE + SPHARMNORM = QPMS_NORMALISATION_NORM_SPHARM + SPHARMNORM_CS = QPMS_NORMALISATION_NORM_SPHARM | QPMS_NORMALISATION_CSPHASE + UNDEF = QPMS_NORMALISATION_UNDEF + +cdef class BaseSpec: + '''Cython wrapper over qpms_vswf_set_spec_t. + + It should be kept immutable. The memory is managed by numpy/cython, not directly by the C functions, therefore + whenever used in other wrapper classes that need the pointer + to qpms_vswf_set_spec_t, remember to set a (private, probably immutable) reference to qpms.basespec to ensure + correct reference counting and garbage collection. + ''' + #cdef qpms_vswf_set_spec_t s # in pxd + #cdef np.ndarray __ilist # in pxd + #cdef const qpms_uvswfi_t[:] __ilist + + def __cinit__(self, *args, **kwargs): + cdef const qpms_uvswfi_t[:] ilist_memview + if len(args) == 0: + if 'lMax' in kwargs.keys(): # if only lMax is specified, create the 'usual' definition in ('E','M') order + lMax = kwargs['lMax'] + my, ny = cycommon.get_mn_y(lMax) + nelem = len(my) + tlist = nelem * (QPMS_VSWF_ELECTRIC,) + nelem * (QPMS_VSWF_MAGNETIC,) + mlist = 2*list(my) + llist = 2*list(ny) + ilist = cycommon.tlm2uvswfi(tlist,llist,mlist) + else: + raise ValueError + else: # len(args) > 0: + ilist = args[0] + #self.__ilist = np.array(args[0], dtype=qpms_uvswfi_t, order='C', copy=True) # FIXME define the dtypes at qpms_cdef.pxd level + self.__ilist = np.array(ilist, dtype=np.ulonglong, order='C', copy=True) + self.__ilist.setflags(write=False) + ilist_memview = self.__ilist + self.s.ilist = &ilist_memview[0] + self.s.n = len(self.__ilist) + self.s.capacity = 0 # is this the best way? + if 'norm' in kwargs.keys(): + self.s.norm = kwargs['norm'] + else: + self.s.norm = (QPMS_NORMALISATION_NORM_POWER | QPMS_NORMALISATION_CSPHASE) + # set the other metadata + cdef qpms_l_t l + self.s.lMax_L = -1 + cdef qpms_m_t m + cdef qpms_vswf_type_t t + for i in range(self.s.n): + if(qpms_uvswfi2tmn(ilist_memview[i], &t, &m, &l) != QPMS_SUCCESS): + raise ValueError("Invalid uvswf index") + if (t == QPMS_VSWF_ELECTRIC): + self.s.lMax_N = max(self.s.lMax_N, l) + elif (t == QPMS_VSWF_MAGNETIC): + self.s.lMax_M = max(self.s.lMax_M, l) + elif (t == QPMS_VSWF_LONGITUDINAL): + self.s.lMax_L = max(self.s.lMax_L, l) + else: + raise ValueError # If this happens, it's probably a bug, as it should have failed already at qpms_uvswfi2tmn + self.s.lMax = max(self.s.lMax, l) + + def tlm(self): + cdef const qpms_uvswfi_t[:] ilist_memview = self.s.ilist + #cdef qpms_vswf_type_t[:] t = np.empty(shape=(self.s.n,), dtype=qpms_vswf_type_t) # does not work, workaround: + cdef size_t i + cdef np.ndarray ta = np.empty(shape=(self.s.n,), dtype=np.intc) + cdef int[:] t = ta + #cdef qpms_l_t[:] l = np.empty(shape=(self.s.n,), dtype=qpms_l_t) # FIXME explicit dtype again + cdef np.ndarray la = np.empty(shape=(self.s.n,), dtype=np.intc) + cdef qpms_l_t[:] l = la + #cdef qpms_m_t[:] m = np.empty(shape=(self.s.n,), dtype=qpms_m_t) # FIXME explicit dtype again + cdef np.ndarray ma = np.empty(shape=(self.s.n,), dtype=np.intc) + cdef qpms_m_t[:] m = ma + for i in range(self.s.n): + qpms_uvswfi2tmn(self.s.ilist[i], &t[i], &m[i], &l[i]) + return (ta, la, ma) + + def m(self): # ugly + return self.tlm()[2] + + def t(self): # ugly + return self.tlm()[0] + + def l(self): # ugly + return self.tlm()[1] + + def __len__(self): + return self.s.n + + def __getitem__(self, key): + # TODO raise correct errors (TypeError on bad type of key, IndexError on exceeding index) + return self.__ilist[key] + + property ilist: + def __get__(self): + return self.__ilist + + cdef qpms_vswf_set_spec_t *rawpointer(BaseSpec self): + '''Pointer to the qpms_vswf_set_spec_t structure. + Don't forget to reference the BaseSpec object itself when storing the pointer anywhere!!! + ''' + return &(self.s) + + property rawpointer: + def __get__(self): + return &(self.s) + + property norm: + def __get__(self): + return VSWFNorm(self.s.norm) + diff --git a/qpms/cycommon.pyx b/qpms/cycommon.pyx new file mode 100644 index 0000000..8b0d9d9 --- /dev/null +++ b/qpms/cycommon.pyx @@ -0,0 +1,182 @@ +import numpy as np +from qpms_cdefs cimport * +cimport cython +import enum + +# Here will be enum and dtype definitions; maybe move these to a separate file +class VSWFType(enum.IntEnum): + ELECTRIC = QPMS_VSWF_ELECTRIC + MAGNETIC = QPMS_VSWF_MAGNETIC + LONGITUDINAL = QPMS_VSWF_LONGITUDINAL + M = QPMS_VSWF_MAGNETIC + N = QPMS_VSWF_ELECTRIC + L = QPMS_VSWF_LONGITUDINAL + +class BesselType(enum.IntEnum): + UNDEF = QPMS_BESSEL_UNDEF + REGULAR = QPMS_BESSEL_REGULAR + SINGULAR = QPMS_BESSEL_SINGULAR + HANKEL_PLUS = QPMS_HANKEL_PLUS + HANKEL_MINUS = QPMS_HANKEL_MINUS + +class PointGroupClass(enum.IntEnum): + CN = QPMS_PGS_CN + S2N = QPMS_PGS_S2N + CNH = QPMS_PGS_CNH + CNV = QPMS_PGS_CNV + DN = QPMS_PGS_DN + DND = QPMS_PGS_DND + DNH = QPMS_PGS_DNH + T = QPMS_PGS_T + TD = QPMS_PGS_TD + TH = QPMS_PGS_TH + O = QPMS_PGS_O + OH = QPMS_PGS_OH + I = QPMS_PGS_I + IH = QPMS_PGS_IH + CINF = QPMS_PGS_CINF + CINFH = QPMS_PGS_CINFH + CINFV = QPMS_PGS_CINFV + DINF = QPMS_PGS_DINF + DINFH = QPMS_PGS_DINFH + SO3 = QPMS_PGS_SO3 + O3 = QPMS_PGS_O3 + +try: + class DebugFlags(enum.IntFlag): # Should be IntFlag if python version >= 3.6 + MISC = QPMS_DBGMSG_MISC + THREADS = QPMS_DBGMSG_THREADS + has_IntFlag = True +except AttributeError: # For old versions of enum, use IntEnum instead + class DebugFlags(enum.IntEnum): + MISC = QPMS_DBGMSG_MISC + THREADS = QPMS_DBGMSG_THREADS + has_IntFlag = False + +def dbgmsg_enable(qpms_dbgmsg_flags types): + flags = qpms_dbgmsg_enable(types) + return DebugFlags(flags) if has_IntFlag else flags +def dbgmsg_disable(qpms_dbgmsg_flags types): + flags = qpms_dbgmsg_disable(types) + return DebugFlags(flags) if has_IntFlag else flags +def dbgmsg_active(): + flags = qpms_dbgmsg_enable(0) + return DebugFlags(flags) if has_IntFlag else flags + +#import re # TODO for crep methods? + +#cimport openmp +#openmp.omp_set_dynamic(1) + +## Auxillary function for retrieving the "meshgrid-like" indices; inc. nmax +@cython.boundscheck(False) +def get_mn_y(int nmax): + """ + Auxillary function for retreiving the 'meshgrid-like' indices from the flat indexing; + inc. nmax. + ('y to mn' conversion) + + Parameters + ---------- + + nmax : int + The maximum order to which the VSWFs / Legendre functions etc. will be evaluated. + + Returns + ------- + + output : (m, n) + Tuple of two arrays of type np.array(shape=(nmax*nmax + 2*nmax), dtype=np.int), + where [(m[y],n[y]) for y in range(nmax*nmax + 2*nma)] covers all possible + integer pairs n >= 1, -n <= m <= n. + """ + cdef Py_ssize_t nelems = nmax * nmax + 2 * nmax + cdef np.ndarray[np.int_t,ndim=1] m_arr = np.empty([nelems], dtype=np.int) + cdef np.ndarray[np.int_t,ndim=1] n_arr = np.empty([nelems], dtype=np.int) + cdef Py_ssize_t i = 0 + cdef np.int_t n, m + for n in range(1,nmax+1): + for m in range(-n,n+1): + m_arr[i] = m + n_arr[i] = n + i = i + 1 + return (m_arr, n_arr) + +def get_nelem(unsigned int lMax): + return lMax * (lMax + 2) + +def get_y_mn_unsigned(int nmax): + """ + Auxillary function for mapping 'unsigned m', n indices to the flat y-indexing. + For use with functions as scipy.special.lpmn, which have to be evaluated separately + for positive and negative m. + + Parameters + ---------- + + nmax : int + The maximum order to which the VSWFs / Legendre functions etc. will be evaluated. + + output : (ymn_plus, ymn_minus) + Tuple of two arrays of shape (nmax+1,nmax+1), containing the flat y-indices corresponding + to the respective (m,n) and (-m,n). The elements for which |m| > n are set to -1. + (Therefore, the caller must not use those elements equal to -1.) + """ + cdef np.ndarray[np.intp_t, ndim=2] ymn_plus = np.full((nmax+1,nmax+1),-1, dtype=np.intp) + cdef np.ndarray[np.intp_t, ndim=2] ymn_minus = np.full((nmax+1,nmax+1),-1, dtype=np.intp) + cdef Py_ssize_t i = 0 + cdef np.int_t n, m + for n in range(1,nmax+1): + for m in range(-n,0): + ymn_minus[-m,n] = i + i = i + 1 + for m in range(0,n+1): + ymn_plus[m,n] = i + i = i + 1 + return(ymn_plus, ymn_minus) + + +def tlm2uvswfi(t, l, m): + ''' TODO doc + And TODO this should rather be an ufunc. + ''' + # Very low-priority TODO: add some types / cythonize + if isinstance(t, int) and isinstance(l, int) and isinstance(m, int): + return qpms_tmn2uvswfi(t, m, l) + elif len(t) == len(l) and len(t) == len(m): + u = list() + for i in range(len(t)): + if not (t[i] % 1 == 0 and l[i] % 1 == 0 and m[i] % 1 == 0): # maybe not the best check possible, though + raise ValueError # TODO error message + u.append(qpms_tmn2uvswfi(t[i],m[i],l[i])) + return u + else: + print(len(t), len(l), len(m)) + raise ValueError("Lengths of the t,l,m arrays must be equal, but they are %d, %d, %d." + % (len(t), len(l), len(m))) + + +def uvswfi2tlm(u): + ''' TODO doc + and TODO this should rather be an ufunc. + ''' + cdef qpms_vswf_type_t t + cdef qpms_l_t l + cdef qpms_m_t m + cdef size_t i + if isinstance(u, (int, np.ulonglong)): + if (qpms_uvswfi2tmn(u, &t, &m, &l) != QPMS_SUCCESS): + raise ValueError("Invalid uvswf index") + return (t, l, m) + else: + ta = list() + la = list() + ma = list() + for i in range(len(u)): + if (qpms_uvswfi2tmn(u[i], &t, &m, &l) != QPMS_SUCCESS): + raise ValueError("Invalid uvswf index") + ta.append(t) + la.append(l) + ma.append(m) + return (ta, la, ma) + diff --git a/qpms/cyquaternions.pxd b/qpms/cyquaternions.pxd new file mode 100644 index 0000000..4855d3d --- /dev/null +++ b/qpms/cyquaternions.pxd @@ -0,0 +1,8 @@ +from qpms_cdefs cimport * + +cdef class CQuat: + cdef readonly qpms_quat_t q + +cdef class IRot3: + cdef readonly qpms_irot3_t qd + cdef void cset(self, qpms_irot3_t qd) diff --git a/qpms/cyquaternions.pyx b/qpms/cyquaternions.pyx new file mode 100644 index 0000000..b57c025 --- /dev/null +++ b/qpms/cyquaternions.pyx @@ -0,0 +1,330 @@ +from cybspec cimport BaseSpec +import cmath +import math + + +def complex_crep(complex c, parentheses = False, shortI = True, has_Imaginary = False): + ''' + Return a C-code compatible string representation of a (python) complex number. + ''' + return ( ('(' if parentheses else '') + + repr(c.real) + + ('+' if math.copysign(1, c.imag) >= 0 else '') + + repr(c.imag) + + ('*I' if shortI else '*_Imaginary_I' if has_Imaginary else '*_Complex_I') + + (')' if parentheses else '') + ) + + +cdef class CQuat: + ''' + Wrapper of the qpms_quat_t object, with the functionality + to evaluate Wigner D-matrix elements. + ''' + # cdef readonly qpms_quat_t q # pxd + + def __cinit__(self, double w, double x, double y, double z): + cdef qpms_quat4d_t p + p.c1 = w + p.ci = x + p.cj = y + p.ck = z + self.q = qpms_quat_2c_from_4d(p) + + def copy(self): + res = CQuat(0,0,0,0) + res.q = self.q + return res + + def __repr__(self): # TODO make this look like a quaternion with i,j,k + return repr(self.r) + + def __add__(CQuat self, CQuat other): + # TODO add real numbers + res = CQuat(0,0,0,0) + res.q = qpms_quat_add(self.q, other.q) + return res + + def __mul__(self, other): + res = CQuat(0,0,0,0) + if isinstance(self, CQuat): + if isinstance(other, CQuat): + res.q = qpms_quat_mult(self.q, other.q) + elif isinstance(other, (int, float)): + res.q = qpms_quat_rscale(other, self.q) + else: return NotImplemented + elif isinstance(self, (int, float)): + if isinstance(other, CQuat): + res.q = qpms_quat_rscale(self, other.q) + else: return NotImplemented + return res + + def __neg__(CQuat self): + res = CQuat(0,0,0,0) + res.q = qpms_quat_rscale(-1, self.q) + return res + + def __sub__(CQuat self, CQuat other): + res = CQuat(0,0,0,0) + res.q = qpms_quat_add(self.q, qpms_quat_rscale(-1,other.q)) + return res + + def __abs__(self): + return qpms_quat_norm(self.q) + + def norm(self): + return qpms_quat_norm(self.q) + + def imnorm(self): + return qpms_quat_imnorm(self.q) + + def exp(self): + res = CQuat(0,0,0,0) + res.q = qpms_quat_exp(self.q) + return res + + def log(self): + res = CQuat(0,0,0,0) + res.q = qpms_quat_exp(self.q) + return res + + def __pow__(CQuat self, double other, _): + res = CQuat(0,0,0,0) + res.q = qpms_quat_pow(self.q, other) + return res + + def normalise(self): + res = CQuat(0,0,0,0) + res.q = qpms_quat_normalise(self.q) + return res + + def isclose(CQuat self, CQuat other, rtol=1e-5, atol=1e-8): + ''' + Checks whether two quaternions are "almost equal". + ''' + return abs(self - other) <= (atol + rtol * abs(other)) + + property c: + ''' + Quaternion representation as two complex numbers + ''' + def __get__(self): + return (self.q.a, self.q.b) + def __set__(self, RaRb): + self.q.a = RaRb[0] + self.q.b = RaRb[1] + + property r: + ''' + Quaternion representation as four real numbers + ''' + def __get__(self): + cdef qpms_quat4d_t p + p = qpms_quat_4d_from_2c(self.q) + return (p.c1, p.ci, p.cj, p.ck) + def __set__(self, wxyz): + cdef qpms_quat4d_t p + p.c1 = wxyz[0] + p.ci = wxyz[1] + p.cj = wxyz[2] + p.ck = wxyz[3] + self.q = qpms_quat_2c_from_4d(p) + + def crepr(self): + ''' + Returns a string that can be used in C code to initialise a qpms_irot3_t + ''' + return '{' + complex_crep(self.q.a) + ', ' + complex_crep(self.q.b) + '}' + + def wignerDelem(self, qpms_l_t l, qpms_m_t mp, qpms_m_t m): + ''' + Returns an element of a bosonic Wigner matrix. + ''' + # don't crash on bad l, m here + 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: + ''' + Wrapper over the C type qpms_irot3_t. + ''' + #cdef readonly qpms_irot3_t qd + + def __cinit__(self, *args): + ''' + TODO doc + ''' + # TODO implement a constructor with + # - tuple as argument ...? + if (len(args) == 0): # no args, return identity + self.qd.rot.a = 1 + self.qd.rot.b = 0 + self.qd.det = 1 + elif (len(args) == 2 and isinstance(args[0], CQuat) and isinstance(args[1], (int, float))): + # The original __cinit__(self, CQuat q, short det) constructor + q = args[0] + det = args[1] + if (det != 1 and det != -1): + raise ValueError("Improper rotation determinant has to be 1 or -1") + self.qd.rot = q.normalise().q + self.qd.det = det + elif (len(args) == 1 and isinstance(args[0], IRot3)): + # Copy + self.qd = args[0].qd + elif (len(args) == 1 and isinstance(args[0], CQuat)): + # proper rotation from a quaternion + q = args[0] + det = 1 + self.qd.rot = q.normalise().q + self.qd.det = det + else: + raise ValueError('Unsupported constructor arguments') + + cdef void cset(self, qpms_irot3_t qd): + self.qd = qd + + def copy(self): + res = IRot3(CQuat(1,0,0,0),1) + res.qd = self.qd + return res + + property rot: + ''' + The proper rotation part of the IRot3 type. + ''' + def __get__(self): + res = CQuat(0,0,0,0) + res.q = self.qd.rot + return res + def __set__(self, CQuat r): + # TODO check for non-zeroness and throw an exception if norm is zero + self.qd.rot = r.normalise().q + + property det: + ''' + The determinant of the improper rotation. + ''' + def __get__(self): + return self.qd.det + def __set__(self, d): + d = int(d) + if (d != 1 and d != -1): + raise ValueError("Improper rotation determinant has to be 1 or -1") + self.qd.det = d + + def __repr__(self): # TODO make this look like a quaternion with i,j,k + return '(' + repr(self.rot) + ', ' + repr(self.det) + ')' + + def crepr(self): + ''' + Returns a string that can be used in C code to initialise a qpms_irot3_t + ''' + return '{' + self.rot.crepr() + ', ' + repr(self.det) + '}' + + def __mul__(IRot3 self, IRot3 other): + res = IRot3(CQuat(1,0,0,0), 1) + res.qd = qpms_irot3_mult(self.qd, other.qd) + return res + + def __pow__(IRot3 self, n, _): + cdef int nint + if (n % 1 == 0): + nint = n + else: + raise ValueError("The exponent of an IRot3 has to have an integer value.") + res = IRot3(CQuat(1,0,0,0), 1) + res.qd = qpms_irot3_pow(self.qd, n) + return res + + def isclose(IRot3 self, IRot3 other, rtol=1e-5, atol=1e-8): + ''' + Checks whether two (improper) rotations are "almost equal". + Returns always False if the determinants are different. + ''' + if self.det != other.det: + return False + return (self.rot.isclose(other.rot, rtol=rtol, atol=atol) + # unit quaternions are a double cover of SO(3), i.e. + # minus the same quaternion represents the same rotation + or self.rot.isclose(-(other.rot), rtol=rtol, atol=atol) + ) + + # Several 'named constructors' for convenience + @staticmethod + def inversion(): + ''' + Returns an IRot3 object representing the 3D spatial inversion. + ''' + r = IRot3() + r.det = -1 + return r + + @staticmethod + def zflip(): + ''' + Returns an IRot3 object representing the 3D xy-plane mirror symmetry (z axis sign flip). + ''' + r = IRot3() + r.rot = CQuat(0,0,0,1) # π-rotation around z-axis + r.det = -1 # inversion + return r + + @staticmethod + def yflip(): + ''' + Returns an IRot3 object representing the 3D xz-plane mirror symmetry (y axis sign flip). + ''' + r = IRot3() + r.rot = CQuat(0,0,1,0) # π-rotation around y-axis + r.det = -1 # inversion + return r + + @staticmethod + def xflip(): + ''' + Returns an IRot3 object representing the 3D yz-plane mirror symmetry (x axis sign flip). + ''' + r = IRot3() + r.rot = CQuat(0,1,0,0) # π-rotation around x-axis + r.det = -1 # inversion + return r + + @staticmethod + def zrotN(int n): + ''' + Returns an IRot3 object representing a \f$ C_n $\f rotation (around the z-axis). + ''' + r = IRot3() + r.rot = CQuat(math.cos(math.pi/n),0,0,math.sin(math.pi/n)) + return r + + @staticmethod + def identity(): + ''' + An alias for the constructor without arguments; returns identity. + ''' + return IRot3() + + def as_uvswf_matrix(IRot3 self, BaseSpec bspec): + ''' + Returns the uvswf representation of the current transform as a numpy array + ''' + cdef ssize_t sz = len(bspec) + cdef np.ndarray m = np.empty((sz, sz), dtype=complex, order='C') # FIXME explicit dtype + cdef cdouble[:, ::1] view = m + qpms_irot3_uvswfi_dense(&view[0,0], bspec.rawpointer(), self.qd) + return m + diff --git a/qpms/qpms_c.pyx b/qpms/qpms_c.pyx index a8e5960..4f1e70d 100644 --- a/qpms/qpms_c.pyx +++ b/qpms/qpms_c.pyx @@ -9,157 +9,18 @@ to make them available in Python. import numpy as np import cmath -from qpms_cdefs cimport * +#from qpms_cdefs cimport * +from cyquaternions cimport * +from cyquaternions import * +from cybspec cimport * +from cybspec import * +from cycommon import * cimport cython from cython.parallel cimport parallel, prange import enum import warnings import os -# Here will be enum and dtype definitions; maybe move these to a separate file -class VSWFType(enum.IntEnum): - ELECTRIC = QPMS_VSWF_ELECTRIC - MAGNETIC = QPMS_VSWF_MAGNETIC - LONGITUDINAL = QPMS_VSWF_LONGITUDINAL - M = QPMS_VSWF_MAGNETIC - N = QPMS_VSWF_ELECTRIC - L = QPMS_VSWF_LONGITUDINAL - -class BesselType(enum.IntEnum): - UNDEF = QPMS_BESSEL_UNDEF - REGULAR = QPMS_BESSEL_REGULAR - SINGULAR = QPMS_BESSEL_SINGULAR - HANKEL_PLUS = QPMS_HANKEL_PLUS - HANKEL_MINUS = QPMS_HANKEL_MINUS - -class PointGroupClass(enum.IntEnum): - CN = QPMS_PGS_CN - S2N = QPMS_PGS_S2N - CNH = QPMS_PGS_CNH - CNV = QPMS_PGS_CNV - DN = QPMS_PGS_DN - DND = QPMS_PGS_DND - DNH = QPMS_PGS_DNH - T = QPMS_PGS_T - TD = QPMS_PGS_TD - TH = QPMS_PGS_TH - O = QPMS_PGS_O - OH = QPMS_PGS_OH - I = QPMS_PGS_I - IH = QPMS_PGS_IH - CINF = QPMS_PGS_CINF - CINFH = QPMS_PGS_CINFH - CINFV = QPMS_PGS_CINFV - DINF = QPMS_PGS_DINF - DINFH = QPMS_PGS_DINFH - SO3 = QPMS_PGS_SO3 - O3 = QPMS_PGS_O3 - -class VSWFNorm(enum.IntEnum): - # TODO try to make this an enum.IntFlag if supported - # TODO add the other flags from qpms_normalisation_t as well - UNNORM = QPMS_NORMALISATION_NORM_NONE - UNNORM_CS = QPMS_NORMALISATION_NORM_NONE | QPMS_NORMALISATION_CSPHASE - POWERNORM = QPMS_NORMALISATION_NORM_POWER - POWERNORM_CS = QPMS_NORMALISATION_NORM_POWER | QPMS_NORMALISATION_CSPHASE - SPHARMNORM = QPMS_NORMALISATION_NORM_SPHARM - SPHARMNORM_CS = QPMS_NORMALISATION_NORM_SPHARM | QPMS_NORMALISATION_CSPHASE - UNDEF = QPMS_NORMALISATION_UNDEF - -try: - class DebugFlags(enum.IntFlag): # Should be IntFlag if python version >= 3.6 - MISC = QPMS_DBGMSG_MISC - THREADS = QPMS_DBGMSG_THREADS - has_IntFlag = True -except AttributeError: # For old versions of enum, use IntEnum instead - class DebugFlags(enum.IntEnum): - MISC = QPMS_DBGMSG_MISC - THREADS = QPMS_DBGMSG_THREADS - has_IntFlag = False - -def dbgmsg_enable(qpms_dbgmsg_flags types): - flags = qpms_dbgmsg_enable(types) - return DebugFlags(flags) if has_IntFlag else flags -def dbgmsg_disable(qpms_dbgmsg_flags types): - flags = qpms_dbgmsg_disable(types) - return DebugFlags(flags) if has_IntFlag else flags -def dbgmsg_active(): - flags = qpms_dbgmsg_enable(0) - return DebugFlags(flags) if has_IntFlag else flags - -import math # for copysign in crep methods -#import re # TODO for crep methods? - -#cimport openmp -#openmp.omp_set_dynamic(1) - -## Auxillary function for retrieving the "meshgrid-like" indices; inc. nmax -@cython.boundscheck(False) -def get_mn_y(int nmax): - """ - Auxillary function for retreiving the 'meshgrid-like' indices from the flat indexing; - inc. nmax. - ('y to mn' conversion) - - Parameters - ---------- - - nmax : int - The maximum order to which the VSWFs / Legendre functions etc. will be evaluated. - - Returns - ------- - - output : (m, n) - Tuple of two arrays of type np.array(shape=(nmax*nmax + 2*nmax), dtype=np.int), - where [(m[y],n[y]) for y in range(nmax*nmax + 2*nma)] covers all possible - integer pairs n >= 1, -n <= m <= n. - """ - cdef Py_ssize_t nelems = nmax * nmax + 2 * nmax - cdef np.ndarray[np.int_t,ndim=1] m_arr = np.empty([nelems], dtype=np.int) - cdef np.ndarray[np.int_t,ndim=1] n_arr = np.empty([nelems], dtype=np.int) - cdef Py_ssize_t i = 0 - cdef np.int_t n, m - for n in range(1,nmax+1): - for m in range(-n,n+1): - m_arr[i] = m - n_arr[i] = n - i = i + 1 - return (m_arr, n_arr) - -def get_nelem(unsigned int lMax): - return lMax * (lMax + 2) - -def get_y_mn_unsigned(int nmax): - """ - Auxillary function for mapping 'unsigned m', n indices to the flat y-indexing. - For use with functions as scipy.special.lpmn, which have to be evaluated separately - for positive and negative m. - - Parameters - ---------- - - nmax : int - The maximum order to which the VSWFs / Legendre functions etc. will be evaluated. - - output : (ymn_plus, ymn_minus) - Tuple of two arrays of shape (nmax+1,nmax+1), containing the flat y-indices corresponding - to the respective (m,n) and (-m,n). The elements for which |m| > n are set to -1. - (Therefore, the caller must not use those elements equal to -1.) - """ - cdef np.ndarray[np.intp_t, ndim=2] ymn_plus = np.full((nmax+1,nmax+1),-1, dtype=np.intp) - cdef np.ndarray[np.intp_t, ndim=2] ymn_minus = np.full((nmax+1,nmax+1),-1, dtype=np.intp) - cdef Py_ssize_t i = 0 - cdef np.int_t n, m - for n in range(1,nmax+1): - for m in range(-n,0): - ymn_minus[-m,n] = i - i = i + 1 - for m in range(0,n+1): - ymn_plus[m,n] = i - i = i + 1 - return(ymn_plus, ymn_minus) - cdef int q_max(int m, int n, int mu, int nu): return min(n,nu,(n+nu-abs(m+mu)//2)) @@ -747,439 +608,6 @@ cdef class trans_calculator: # TODO make possible to access the attributes (to show normalization etc) -def complex_crep(complex c, parentheses = False, shortI = True, has_Imaginary = False): - ''' - Return a C-code compatible string representation of a (python) complex number. - ''' - return ( ('(' if parentheses else '') - + repr(c.real) - + ('+' if math.copysign(1, c.imag) >= 0 else '') - + repr(c.imag) - + ('*I' if shortI else '*_Imaginary_I' if has_Imaginary else '*_Complex_I') - + (')' if parentheses else '') - ) - -cdef class BaseSpec: - '''Cython wrapper over qpms_vswf_set_spec_t. - - It should be kept immutable. The memory is managed by numpy/cython, not directly by the C functions, therefore - whenever used in other wrapper classes that need the pointer - to qpms_vswf_set_spec_t, remember to set a (private, probably immutable) reference to qpms.basespec to ensure - correct reference counting and garbage collection. - ''' - cdef qpms_vswf_set_spec_t s - cdef np.ndarray __ilist - #cdef const qpms_uvswfi_t[:] __ilist - - def __cinit__(self, *args, **kwargs): - cdef const qpms_uvswfi_t[:] ilist_memview - if len(args) == 0: - if 'lMax' in kwargs.keys(): # if only lMax is specified, create the 'usual' definition in ('E','M') order - lMax = kwargs['lMax'] - my, ny = get_mn_y(lMax) - nelem = len(my) - tlist = nelem * (QPMS_VSWF_ELECTRIC,) + nelem * (QPMS_VSWF_MAGNETIC,) - mlist = 2*list(my) - llist = 2*list(ny) - ilist = tlm2uvswfi(tlist,llist,mlist) - else: - raise ValueError - else: # len(args) > 0: - ilist = args[0] - #self.__ilist = np.array(args[0], dtype=qpms_uvswfi_t, order='C', copy=True) # FIXME define the dtypes at qpms_cdef.pxd level - self.__ilist = np.array(ilist, dtype=np.ulonglong, order='C', copy=True) - self.__ilist.setflags(write=False) - ilist_memview = self.__ilist - self.s.ilist = &ilist_memview[0] - self.s.n = len(self.__ilist) - self.s.capacity = 0 # is this the best way? - if 'norm' in kwargs.keys(): - self.s.norm = kwargs['norm'] - else: - self.s.norm = (QPMS_NORMALISATION_NORM_POWER | QPMS_NORMALISATION_CSPHASE) - # set the other metadata - cdef qpms_l_t l - self.s.lMax_L = -1 - cdef qpms_m_t m - cdef qpms_vswf_type_t t - for i in range(self.s.n): - if(qpms_uvswfi2tmn(ilist_memview[i], &t, &m, &l) != QPMS_SUCCESS): - raise ValueError("Invalid uvswf index") - if (t == QPMS_VSWF_ELECTRIC): - self.s.lMax_N = max(self.s.lMax_N, l) - elif (t == QPMS_VSWF_MAGNETIC): - self.s.lMax_M = max(self.s.lMax_M, l) - elif (t == QPMS_VSWF_LONGITUDINAL): - self.s.lMax_L = max(self.s.lMax_L, l) - else: - raise ValueError # If this happens, it's probably a bug, as it should have failed already at qpms_uvswfi2tmn - self.s.lMax = max(self.s.lMax, l) - - def tlm(self): - cdef const qpms_uvswfi_t[:] ilist_memview = self.s.ilist - #cdef qpms_vswf_type_t[:] t = np.empty(shape=(self.s.n,), dtype=qpms_vswf_type_t) # does not work, workaround: - cdef size_t i - cdef np.ndarray ta = np.empty(shape=(self.s.n,), dtype=np.intc) - cdef int[:] t = ta - #cdef qpms_l_t[:] l = np.empty(shape=(self.s.n,), dtype=qpms_l_t) # FIXME explicit dtype again - cdef np.ndarray la = np.empty(shape=(self.s.n,), dtype=np.intc) - cdef qpms_l_t[:] l = la - #cdef qpms_m_t[:] m = np.empty(shape=(self.s.n,), dtype=qpms_m_t) # FIXME explicit dtype again - cdef np.ndarray ma = np.empty(shape=(self.s.n,), dtype=np.intc) - cdef qpms_m_t[:] m = ma - for i in range(self.s.n): - qpms_uvswfi2tmn(self.s.ilist[i], &t[i], &m[i], &l[i]) - return (ta, la, ma) - - def m(self): # ugly - return self.tlm()[2] - - def t(self): # ugly - return self.tlm()[0] - - def l(self): # ugly - return self.tlm()[1] - - def __len__(self): - return self.s.n - - def __getitem__(self, key): - # TODO raise correct errors (TypeError on bad type of key, IndexError on exceeding index) - return self.__ilist[key] - - property ilist: - def __get__(self): - return self.__ilist - - cdef qpms_vswf_set_spec_t *rawpointer(BaseSpec self): - '''Pointer to the qpms_vswf_set_spec_t structure. - Don't forget to reference the BaseSpec object itself when storing the pointer anywhere!!! - ''' - return &(self.s) - - property rawpointer: - def __get__(self): - return &(self.s) - - property norm: - def __get__(self): - return VSWFNorm(self.s.norm) - -# Quaternions from quaternions.h -# (mainly for testing; use moble's quaternions in python) - -cdef class CQuat: - ''' - Wrapper of the qpms_quat_t object, with the functionality - to evaluate Wigner D-matrix elements. - ''' - cdef readonly qpms_quat_t q - - def __cinit__(self, double w, double x, double y, double z): - cdef qpms_quat4d_t p - p.c1 = w - p.ci = x - p.cj = y - p.ck = z - self.q = qpms_quat_2c_from_4d(p) - - def copy(self): - res = CQuat(0,0,0,0) - res.q = self.q - return res - - def __repr__(self): # TODO make this look like a quaternion with i,j,k - return repr(self.r) - - def __add__(CQuat self, CQuat other): - # TODO add real numbers - res = CQuat(0,0,0,0) - res.q = qpms_quat_add(self.q, other.q) - return res - - def __mul__(self, other): - res = CQuat(0,0,0,0) - if isinstance(self, CQuat): - if isinstance(other, CQuat): - res.q = qpms_quat_mult(self.q, other.q) - elif isinstance(other, (int, float)): - res.q = qpms_quat_rscale(other, self.q) - else: return NotImplemented - elif isinstance(self, (int, float)): - if isinstance(other, CQuat): - res.q = qpms_quat_rscale(self, other.q) - else: return NotImplemented - return res - - def __neg__(CQuat self): - res = CQuat(0,0,0,0) - res.q = qpms_quat_rscale(-1, self.q) - return res - - def __sub__(CQuat self, CQuat other): - res = CQuat(0,0,0,0) - res.q = qpms_quat_add(self.q, qpms_quat_rscale(-1,other.q)) - return res - - def __abs__(self): - return qpms_quat_norm(self.q) - - def norm(self): - return qpms_quat_norm(self.q) - - def imnorm(self): - return qpms_quat_imnorm(self.q) - - def exp(self): - res = CQuat(0,0,0,0) - res.q = qpms_quat_exp(self.q) - return res - - def log(self): - res = CQuat(0,0,0,0) - res.q = qpms_quat_exp(self.q) - return res - - def __pow__(CQuat self, double other, _): - res = CQuat(0,0,0,0) - res.q = qpms_quat_pow(self.q, other) - return res - - def normalise(self): - res = CQuat(0,0,0,0) - res.q = qpms_quat_normalise(self.q) - return res - - def isclose(CQuat self, CQuat other, rtol=1e-5, atol=1e-8): - ''' - Checks whether two quaternions are "almost equal". - ''' - return abs(self - other) <= (atol + rtol * abs(other)) - - property c: - ''' - Quaternion representation as two complex numbers - ''' - def __get__(self): - return (self.q.a, self.q.b) - def __set__(self, RaRb): - self.q.a = RaRb[0] - self.q.b = RaRb[1] - - property r: - ''' - Quaternion representation as four real numbers - ''' - def __get__(self): - cdef qpms_quat4d_t p - p = qpms_quat_4d_from_2c(self.q) - return (p.c1, p.ci, p.cj, p.ck) - def __set__(self, wxyz): - cdef qpms_quat4d_t p - p.c1 = wxyz[0] - p.ci = wxyz[1] - p.cj = wxyz[2] - p.ck = wxyz[3] - self.q = qpms_quat_2c_from_4d(p) - - def crepr(self): - ''' - Returns a string that can be used in C code to initialise a qpms_irot3_t - ''' - return '{' + complex_crep(self.q.a) + ', ' + complex_crep(self.q.b) + '}' - - def wignerDelem(self, qpms_l_t l, qpms_m_t mp, qpms_m_t m): - ''' - Returns an element of a bosonic Wigner matrix. - ''' - # don't crash on bad l, m here - 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: - ''' - Wrapper over the C type qpms_irot3_t. - ''' - cdef readonly qpms_irot3_t qd - - def __cinit__(self, *args): - ''' - TODO doc - ''' - # TODO implement a constructor with - # - tuple as argument ...? - if (len(args) == 0): # no args, return identity - self.qd.rot.a = 1 - self.qd.rot.b = 0 - self.qd.det = 1 - elif (len(args) == 2 and isinstance(args[0], CQuat) and isinstance(args[1], (int, float))): - # The original __cinit__(self, CQuat q, short det) constructor - q = args[0] - det = args[1] - if (det != 1 and det != -1): - raise ValueError("Improper rotation determinant has to be 1 or -1") - self.qd.rot = q.normalise().q - self.qd.det = det - elif (len(args) == 1 and isinstance(args[0], IRot3)): - # Copy - self.qd = args[0].qd - elif (len(args) == 1 and isinstance(args[0], CQuat)): - # proper rotation from a quaternion - q = args[0] - det = 1 - self.qd.rot = q.normalise().q - self.qd.det = det - else: - raise ValueError('Unsupported constructor arguments') - - cdef void cset(self, qpms_irot3_t qd): - self.qd = qd - - def copy(self): - res = IRot3(CQuat(1,0,0,0),1) - res.qd = self.qd - return res - - property rot: - ''' - The proper rotation part of the IRot3 type. - ''' - def __get__(self): - res = CQuat(0,0,0,0) - res.q = self.qd.rot - return res - def __set__(self, CQuat r): - # TODO check for non-zeroness and throw an exception if norm is zero - self.qd.rot = r.normalise().q - - property det: - ''' - The determinant of the improper rotation. - ''' - def __get__(self): - return self.qd.det - def __set__(self, d): - d = int(d) - if (d != 1 and d != -1): - raise ValueError("Improper rotation determinant has to be 1 or -1") - self.qd.det = d - - def __repr__(self): # TODO make this look like a quaternion with i,j,k - return '(' + repr(self.rot) + ', ' + repr(self.det) + ')' - - def crepr(self): - ''' - Returns a string that can be used in C code to initialise a qpms_irot3_t - ''' - return '{' + self.rot.crepr() + ', ' + repr(self.det) + '}' - - def __mul__(IRot3 self, IRot3 other): - res = IRot3(CQuat(1,0,0,0), 1) - res.qd = qpms_irot3_mult(self.qd, other.qd) - return res - - def __pow__(IRot3 self, n, _): - cdef int nint - if (n % 1 == 0): - nint = n - else: - raise ValueError("The exponent of an IRot3 has to have an integer value.") - res = IRot3(CQuat(1,0,0,0), 1) - res.qd = qpms_irot3_pow(self.qd, n) - return res - - def isclose(IRot3 self, IRot3 other, rtol=1e-5, atol=1e-8): - ''' - Checks whether two (improper) rotations are "almost equal". - Returns always False if the determinants are different. - ''' - if self.det != other.det: - return False - return (self.rot.isclose(other.rot, rtol=rtol, atol=atol) - # unit quaternions are a double cover of SO(3), i.e. - # minus the same quaternion represents the same rotation - or self.rot.isclose(-(other.rot), rtol=rtol, atol=atol) - ) - - # Several 'named constructors' for convenience - @staticmethod - def inversion(): - ''' - Returns an IRot3 object representing the 3D spatial inversion. - ''' - r = IRot3() - r.det = -1 - return r - - @staticmethod - def zflip(): - ''' - Returns an IRot3 object representing the 3D xy-plane mirror symmetry (z axis sign flip). - ''' - r = IRot3() - r.rot = CQuat(0,0,0,1) # π-rotation around z-axis - r.det = -1 # inversion - return r - - @staticmethod - def yflip(): - ''' - Returns an IRot3 object representing the 3D xz-plane mirror symmetry (y axis sign flip). - ''' - r = IRot3() - r.rot = CQuat(0,0,1,0) # π-rotation around y-axis - r.det = -1 # inversion - return r - - @staticmethod - def xflip(): - ''' - Returns an IRot3 object representing the 3D yz-plane mirror symmetry (x axis sign flip). - ''' - r = IRot3() - r.rot = CQuat(0,1,0,0) # π-rotation around x-axis - r.det = -1 # inversion - return r - - @staticmethod - def zrotN(int n): - ''' - Returns an IRot3 object representing a \f$ C_n $\f rotation (around the z-axis). - ''' - r = IRot3() - r.rot = CQuat(math.cos(math.pi/n),0,0,math.sin(math.pi/n)) - return r - - @staticmethod - def identity(): - ''' - An alias for the constructor without arguments; returns identity. - ''' - return IRot3() - - def as_uvswf_matrix(IRot3 self, BaseSpec bspec): - ''' - Returns the uvswf representation of the current transform as a numpy array - ''' - cdef ssize_t sz = len(bspec) - cdef np.ndarray m = np.empty((sz, sz), dtype=complex, order='C') # FIXME explicit dtype - cdef cdouble[:, ::1] view = m - qpms_irot3_uvswfi_dense(&view[0,0], bspec.rawpointer(), self.qd) - return m - cdef class MaterialInterpolator: ''' Wrapper over the qpms_permittivity_interpolator_t structure. @@ -1833,49 +1261,3 @@ cdef class ScatteringMatrix: return f -def tlm2uvswfi(t, l, m): - ''' TODO doc - And TODO this should rather be an ufunc. - ''' - # Very low-priority TODO: add some types / cythonize - if isinstance(t, int) and isinstance(l, int) and isinstance(m, int): - return qpms_tmn2uvswfi(t, m, l) - elif len(t) == len(l) and len(t) == len(m): - u = list() - for i in range(len(t)): - if not (t[i] % 1 == 0 and l[i] % 1 == 0 and m[i] % 1 == 0): # maybe not the best check possible, though - raise ValueError # TODO error message - u.append(qpms_tmn2uvswfi(t[i],m[i],l[i])) - return u - else: - print(len(t), len(l), len(m)) - raise ValueError("Lengths of the t,l,m arrays must be equal, but they are %d, %d, %d." - % (len(t), len(l), len(m))) - - -def uvswfi2tlm(u): - ''' TODO doc - and TODO this should rather be an ufunc. - ''' - cdef qpms_vswf_type_t t - cdef qpms_l_t l - cdef qpms_m_t m - cdef size_t i - if isinstance(u, (int, np.ulonglong)): - if (qpms_uvswfi2tmn(u, &t, &m, &l) != QPMS_SUCCESS): - raise ValueError("Invalid uvswf index") - return (t, l, m) - else: - ta = list() - la = list() - ma = list() - for i in range(len(u)): - if (qpms_uvswfi2tmn(u[i], &t, &m, &l) != QPMS_SUCCESS): - raise ValueError("Invalid uvswf index") - ta.append(t) - la.append(l) - ma.append(m) - return (ta, la, ma) - - - diff --git a/setup.py b/setup.py index 1e3ad14..3d80b9e 100755 --- a/setup.py +++ b/setup.py @@ -62,7 +62,12 @@ amos_sources = [ qpms_c = Extension('qpms_c', - sources = ['qpms/qpms_c.pyx', #'qpms/hexpoints_c.pyx', + sources = [ + 'qpms/cycommon.pyx', + 'qpms/cyquaternions.pyx', + 'qpms/cybspec.pyx', + 'qpms/qpms_c.pyx', + #'qpms/hexpoints_c.pyx', 'qpms/gaunt.c',#'qpms/gaunt.h','qpms/vectors.h','qpms/translations.h', # FIXME http://stackoverflow.com/questions/4259170/python-setup-script-extensions-how-do-you-include-a-h-file 'qpms/translations.c',