From 1aa64248fcfd10e11f11ebbbe7ef7b8d53432aeb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Sat, 10 Aug 2019 09:55:50 +0300 Subject: [PATCH] Split away CTMatrix, TMatrixInterpolator Former-commit-id: ca92a0938f12fda858794b203e196366f2466a6f --- qpms/__init__.py | 1 + qpms/cytmatrices.pxd | 25 +++++++++ qpms/cytmatrices.pyx | 116 ++++++++++++++++++++++++++++++++++++++ qpms/qpms_c.pxd | 8 +++ qpms/qpms_c.pyx | 130 +------------------------------------------ setup.py | 7 ++- 6 files changed, 157 insertions(+), 130 deletions(-) create mode 100644 qpms/cytmatrices.pxd create mode 100644 qpms/cytmatrices.pyx create mode 100644 qpms/qpms_c.pxd diff --git a/qpms/__init__.py b/qpms/__init__.py index 1bd8c74..d978158 100644 --- a/qpms/__init__.py +++ b/qpms/__init__.py @@ -5,6 +5,7 @@ from .qpms_c import * from .qpms_p import * from .cyquaternions import CQuat, IRot3 from .cybspec import VSWFNorm, BaseSpec +from .cytmatrices import CTMatrix, TMatrixInterpolator from .cytranslations import trans_calculator from .cymaterials import MaterialInterpolator from .lattices2d import * diff --git a/qpms/cytmatrices.pxd b/qpms/cytmatrices.pxd new file mode 100644 index 0000000..b267d85 --- /dev/null +++ b/qpms/cytmatrices.pxd @@ -0,0 +1,25 @@ +cimport numpy as np +from qpms_cdefs cimport qpms_tmatrix_t, cdouble, qpms_tmatrix_interpolator_t +from cybspec cimport BaseSpec + +cdef class TMatrixInterpolator: + #cdef readonly np.ndarray m # Numpy array holding the matrix data + cdef readonly BaseSpec spec # Here we hold the base spec for the correct reference counting; TODO check if it gets copied + cdef qpms_tmatrix_t *tmatrices_array + cdef cdouble *tmdata + cdef double *freqs + cdef double *freqs_su + cdef size_t nfreqs + cdef qpms_tmatrix_interpolator_t *interp + +cdef class CTMatrix: # N.B. there is another type called TMatrix in tmatrices.py! + cdef readonly np.ndarray m # Numpy array holding the matrix data + cdef readonly BaseSpec spec # Here we hold the base spec for the correct reference counting; TODO check if it gets copied + cdef qpms_tmatrix_t t + + + cdef inline qpms_tmatrix_t *rawpointer(CTMatrix self): + '''Pointer to the qpms_tmatrix_t structure. + Don't forget to reference the BaseSpec object itself when storing the pointer anywhere!!! + ''' + return &(self.t) diff --git a/qpms/cytmatrices.pyx b/qpms/cytmatrices.pyx new file mode 100644 index 0000000..c14c088 --- /dev/null +++ b/qpms/cytmatrices.pyx @@ -0,0 +1,116 @@ +import numpy as np +from qpms_cdefs cimport * +from cybspec cimport BaseSpec +from cycommon import * +from cycommon cimport make_c_string +from qpms_c cimport FinitePointGroup +import warnings +import os +from libc.stdlib cimport free + +cdef class TMatrixInterpolator: + ''' + Wrapper over the qpms_tmatrix_interpolator_t structure. + ''' + def __cinit__(self, filename, BaseSpec bspec, *args, **kwargs): + '''Creates a T-matrix interpolator object from a scuff-tmatrix output''' + global qpms_load_scuff_tmatrix_crash_on_failure + qpms_load_scuff_tmatrix_crash_on_failure = False + self.spec = bspec + cdef char * cpath = make_c_string(filename) + retval = qpms_load_scuff_tmatrix(cpath, self.spec.rawpointer(), + &(self.nfreqs), &(self.freqs), &(self.freqs_su), + &(self.tmatrices_array), &(self.tmdata)) + if (retval != QPMS_SUCCESS): + raise IOError("Could not read T-matrix from %s: %s" % (filename, os.strerror(retval))) + if 'symmetrise' in kwargs: + sym = kwargs['symmetrise'] + if isinstance(sym, FinitePointGroup): + if QPMS_SUCCESS != qpms_symmetrise_tmdata_finite_group( + self.tmdata, self.nfreqs, self.spec.rawpointer(), + (sym).rawpointer()): + raise Exception("This should not happen.") + atol = kwargs['atol'] if 'atol' in kwargs else 1e-16 + qpms_czero_roundoff_clean(self.tmdata, self.nfreqs * len(bspec)**2, atol) + else: + warnings.warn('symmetrise argument type not supported; ignoring.') + self.interp = qpms_tmatrix_interpolator_create(self.nfreqs, + self.freqs, self.tmatrices_array, gsl_interp_cspline) + if not self.interp: raise Exception("Unexpected NULL at interpolator creation.") + def __call__(self, double freq): + '''Returns a TMatrix instance, corresponding to a given frequency.''' + if freq < self.freqs[0] or freq > self.freqs[self.nfreqs-1]:# FIXME here I assume that the input is already sorted + raise ValueError("input frequency %g is outside the interpolator domain (%g, %g)" + % (freq, self.freqs[0], self.freqs[self.nfreqs-1])) + # This is a bit stupid, I should rethink the CTMatrix constuctors + cdef qpms_tmatrix_t *t = qpms_tmatrix_interpolator_eval(self.interp, freq) + cdef CTMatrix res = CTMatrix(self.spec, (t[0].m)) + qpms_tmatrix_free(t) + return res + def __dealloc__(self): + qpms_tmatrix_interpolator_free(self.interp) + free(self.tmatrices_array) + free(self.tmdata) + free(self.freqs_su) + free(self.freqs) + property freq_interval: + def __get__(self): + return [self.freqs[0], self.freqs[self.nfreqs-1]] + +cdef class CTMatrix: # N.B. there is another type called TMatrix in tmatrices.py! + ''' + Wrapper over the C qpms_tmatrix_t stucture. + ''' + + def __cinit__(CTMatrix self, BaseSpec spec, matrix): + self.spec = spec + self.t.spec = self.spec.rawpointer(); + if (matrix is None) or not np.any(matrix): + self.m = np.zeros((len(spec),len(spec)), dtype=complex, order='C') + else: + # The following will raise an exception if shape is wrong + self.m = np.array(matrix, dtype=complex, copy=True, order='C').reshape((len(spec), len(spec))) + #self.m.setflags(write=False) # checkme + cdef cdouble[:,:] m_memview = self.m + self.t.m = &(m_memview[0,0]) + self.t.owns_m = False # Memory in self.t.m is "owned" by self.m, not by self.t... + + property rawpointer: + def __get__(self): + return &(self.t) + + # Transparent access to the T-matrix elements. + def __getitem__(self, key): + return self.m[key] + def __setitem__(self, key, value): + self.m[key] = value + + def as_ndarray(CTMatrix self): + ''' Returns a copy of the T-matrix as a numpy array.''' + # Maybe not totally needed after all, as np.array(T[...]) should be equivalent and not longer + return np.array(self.m, copy=True) + + def spherical_fill(CTMatrix self, double radius, cdouble k_int, + cdouble k_ext, cdouble mu_int = 1, cdouble mu_ext = 1): + ''' Replaces the contents of the T-matrix with those of a spherical particle.''' + qpms_tmatrix_spherical_fill(&self.t, radius, k_int, k_ext, mu_int, mu_ext) + + def spherical_perm_fill(CTMatrix self, double radius, double freq, cdouble epsilon_int, + cdouble epsilon_ext): + '''Replaces the contenst of the T-matrix with those of a spherical particle.''' + qpms_tmatrix_spherical_mu0_fill(&self.t, radius, freq, epsilon_int, epsilon_ext) + + @staticmethod + def spherical(BaseSpec spec, double radius, cdouble k_int, cdouble k_ext, + cdouble mu_int = 1, cdouble mu_ext = 1): + ''' Creates a T-matrix of a spherical nanoparticle. ''' + tm = CTMatrix(spec, 0) + tm.spherical_fill(radius, k_int, k_ext, mu_int, mu_ext) + return tm + + @staticmethod + def spherical_perm(BaseSpec spec, double radius, double freq, cdouble epsilon_int, cdouble epsilon_ext): + '''Creates a T-matrix of a spherical nanoparticle.''' + tm = CTMatrix(spec, 0) + tm.spherical_perm_fill(radius, freq, epsilon_int, epsilon_ext) + return tm diff --git a/qpms/qpms_c.pxd b/qpms/qpms_c.pxd new file mode 100644 index 0000000..19424f0 --- /dev/null +++ b/qpms/qpms_c.pxd @@ -0,0 +1,8 @@ +from qpms_cdefs cimport qpms_finite_group_t + +cdef class FinitePointGroup: + cdef readonly bint owns_data + cdef qpms_finite_group_t *G + + cdef inline qpms_finite_group_t *rawpointer(self): + return self.G diff --git a/qpms/qpms_c.pyx b/qpms/qpms_c.pyx index 7f37b2d..9da3624 100644 --- a/qpms/qpms_c.pyx +++ b/qpms/qpms_c.pyx @@ -16,136 +16,13 @@ from cybspec cimport * #from cybspec import * from cycommon import * from cycommon cimport make_c_string +from cytmatrices cimport CTMatrix cimport cython import enum import warnings import os from libc.stdlib cimport malloc, free, calloc, abort -cdef class TMatrixInterpolator: - ''' - Wrapper over the qpms_tmatrix_interpolator_t structure. - ''' - #cdef readonly np.ndarray m # Numpy array holding the matrix data - cdef readonly BaseSpec spec # Here we hold the base spec for the correct reference counting; TODO check if it gets copied - cdef qpms_tmatrix_t *tmatrices_array - cdef cdouble *tmdata - cdef double *freqs - cdef double *freqs_su - cdef size_t nfreqs - cdef qpms_tmatrix_interpolator_t *interp - - def __cinit__(self, filename, BaseSpec bspec, *args, **kwargs): - '''Creates a T-matrix interpolator object from a scuff-tmatrix output''' - global qpms_load_scuff_tmatrix_crash_on_failure - qpms_load_scuff_tmatrix_crash_on_failure = False - self.spec = bspec - cdef char * cpath = make_c_string(filename) - retval = qpms_load_scuff_tmatrix(cpath, self.spec.rawpointer(), - &(self.nfreqs), &(self.freqs), &(self.freqs_su), - &(self.tmatrices_array), &(self.tmdata)) - if (retval != QPMS_SUCCESS): - raise IOError("Could not read T-matrix from %s: %s" % (filename, os.strerror(retval))) - if 'symmetrise' in kwargs: - sym = kwargs['symmetrise'] - if isinstance(sym, FinitePointGroup): - if QPMS_SUCCESS != qpms_symmetrise_tmdata_finite_group( - self.tmdata, self.nfreqs, self.spec.rawpointer(), - (sym).rawpointer()): - raise Exception("This should not happen.") - atol = kwargs['atol'] if 'atol' in kwargs else 1e-16 - qpms_czero_roundoff_clean(self.tmdata, self.nfreqs * len(bspec)**2, atol) - else: - warnings.warn('symmetrise argument type not supported; ignoring.') - self.interp = qpms_tmatrix_interpolator_create(self.nfreqs, - self.freqs, self.tmatrices_array, gsl_interp_cspline) - if not self.interp: raise Exception("Unexpected NULL at interpolator creation.") - def __call__(self, double freq): - '''Returns a TMatrix instance, corresponding to a given frequency.''' - if freq < self.freqs[0] or freq > self.freqs[self.nfreqs-1]:# FIXME here I assume that the input is already sorted - raise ValueError("input frequency %g is outside the interpolator domain (%g, %g)" - % (freq, self.freqs[0], self.freqs[self.nfreqs-1])) - # This is a bit stupid, I should rethink the CTMatrix constuctors - cdef qpms_tmatrix_t *t = qpms_tmatrix_interpolator_eval(self.interp, freq) - cdef CTMatrix res = CTMatrix(self.spec, (t[0].m)) - qpms_tmatrix_free(t) - return res - def __dealloc__(self): - qpms_tmatrix_interpolator_free(self.interp) - free(self.tmatrices_array) - free(self.tmdata) - free(self.freqs_su) - free(self.freqs) - property freq_interval: - def __get__(self): - return [self.freqs[0], self.freqs[self.nfreqs-1]] - -cdef class CTMatrix: # N.B. there is another type called TMatrix in tmatrices.py! - ''' - Wrapper over the C qpms_tmatrix_t stucture. - ''' - cdef readonly np.ndarray m # Numpy array holding the matrix data - cdef readonly BaseSpec spec # Here we hold the base spec for the correct reference counting; TODO check if it gets copied - cdef qpms_tmatrix_t t - - def __cinit__(CTMatrix self, BaseSpec spec, matrix): - self.spec = spec - self.t.spec = self.spec.rawpointer(); - if (matrix is None) or not np.any(matrix): - self.m = np.zeros((len(spec),len(spec)), dtype=complex, order='C') - else: - # The following will raise an exception if shape is wrong - self.m = np.array(matrix, dtype=complex, copy=True, order='C').reshape((len(spec), len(spec))) - #self.m.setflags(write=False) # checkme - cdef cdouble[:,:] m_memview = self.m - self.t.m = &(m_memview[0,0]) - self.t.owns_m = False # Memory in self.t.m is "owned" by self.m, not by self.t... - - cdef qpms_tmatrix_t *rawpointer(CTMatrix self): - '''Pointer to the qpms_tmatrix_t structure. - Don't forget to reference the BaseSpec object itself when storing the pointer anywhere!!! - ''' - return &(self.t) - property rawpointer: - def __get__(self): - return &(self.t) - - # Transparent access to the T-matrix elements. - def __getitem__(self, key): - return self.m[key] - def __setitem__(self, key, value): - self.m[key] = value - - def as_ndarray(CTMatrix self): - ''' Returns a copy of the T-matrix as a numpy array.''' - # Maybe not totally needed after all, as np.array(T[...]) should be equivalent and not longer - return np.array(self.m, copy=True) - - def spherical_fill(CTMatrix self, double radius, cdouble k_int, - cdouble k_ext, cdouble mu_int = 1, cdouble mu_ext = 1): - ''' Replaces the contents of the T-matrix with those of a spherical particle.''' - qpms_tmatrix_spherical_fill(&self.t, radius, k_int, k_ext, mu_int, mu_ext) - - def spherical_perm_fill(CTMatrix self, double radius, double freq, cdouble epsilon_int, - cdouble epsilon_ext): - '''Replaces the contenst of the T-matrix with those of a spherical particle.''' - qpms_tmatrix_spherical_mu0_fill(&self.t, radius, freq, epsilon_int, epsilon_ext) - - @staticmethod - def spherical(BaseSpec spec, double radius, cdouble k_int, cdouble k_ext, - cdouble mu_int = 1, cdouble mu_ext = 1): - ''' Creates a T-matrix of a spherical nanoparticle. ''' - tm = CTMatrix(spec, 0) - tm.spherical_fill(radius, k_int, k_ext, mu_int, mu_ext) - return tm - - @staticmethod - def spherical_perm(BaseSpec spec, double radius, double freq, cdouble epsilon_int, cdouble epsilon_ext): - '''Creates a T-matrix of a spherical nanoparticle.''' - tm = CTMatrix(spec, 0) - tm.spherical_perm_fill(radius, freq, epsilon_int, epsilon_ext) - return tm - cdef class PointGroup: cdef readonly qpms_pointgroup_t G @@ -196,8 +73,6 @@ cdef class FinitePointGroup: TODO more functionality to make it better usable in Python (group element class at least) ''' - cdef readonly bint owns_data - cdef qpms_finite_group_t *G def __cinit__(self, info): '''Constructs a FinitePointGroup from PointGroupInfo''' @@ -284,9 +159,6 @@ cdef class FinitePointGroup: self.G = 0 self.owns_data = False - cdef qpms_finite_group_t *rawpointer(self): - return self.G - cdef class FinitePointGroupElement: '''TODO''' cdef readonly FinitePointGroup G diff --git a/setup.py b/setup.py index 43a1435..3fbea8e 100755 --- a/setup.py +++ b/setup.py @@ -87,6 +87,11 @@ cycommon = Extension('qpms.cycommon', extra_link_args=['qpms/libqpms.a'], libraries=['gsl', 'lapacke', 'blas', 'gslcblas', 'pthread',] ) +cytmatrices = Extension('qpms.cytmatrices', + sources = ['qpms/cytmatrices.pyx'], + extra_link_args=['qpms/libqpms.a', 'amos/libamos.a'], + libraries=['gsl', 'lapacke', 'blas', 'gslcblas', 'pthread',] + ) cytranslations = Extension('qpms.cytranslations', sources = ['qpms/cytranslations.pyx', 'qpms/translations_python.c', @@ -135,7 +140,7 @@ setup(name='qpms', #'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, cytranslations, cycommon, cyquaternions, cybspec, cymaterials], include_path=['qpms', 'amos'], gdb_debug=True), + ext_modules=cythonize([qpms_c, cytranslations, cytmatrices, cycommon, cyquaternions, cybspec, cymaterials], include_path=['qpms', 'amos'], gdb_debug=True), cmdclass = {'build_ext': build_ext}, zip_safe=False )