QPMS
/
qpms
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Start splitting qpms_c.pyx 

Former-commit-id: bb2e68dc4cb7f85769ddaf5533298ab1f0e84f5b
This commit is contained in:
Marek Nečada 2019-08-09 21:54:13 +03:00
parent 1ef0c0ad4e
commit 371a8a5f7c
8 changed files with 672 additions and 625 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
cyquat.pxd Normal file
View File

@ -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)

9
qpms/cybspec.pxd Normal file
View File

@ -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)

121
qpms/cybspec.pyx Normal file
View File

@ -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_t>(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 = <qpms_uvswfi_t[:self.s.n]> 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], <qpms_vswf_type_t*>&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 <uintptr_t> &(self.s)
property norm:
def __get__(self):
return VSWFNorm(self.s.norm)

182
qpms/cycommon.pyx Normal file
View File

@ -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(<qpms_dbgmsg_flags>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)

8
qpms/cyquaternions.pxd Normal file
View File

@ -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)

330
qpms/cyquaternions.pyx Normal file
View File

@ -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

630
qpms/qpms_c.pyx
View File

@ -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(<qpms_dbgmsg_flags>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_t>(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 = <qpms_uvswfi_t[:self.s.n]> 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], <qpms_vswf_type_t*>&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 <uintptr_t> &(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)

7
setup.py
View File

@ -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',