Split away CTMatrix, TMatrixInterpolator

Former-commit-id: ca92a0938f12fda858794b203e196366f2466a6f

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 *

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)

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(),
+                        (<FinitePointGroup?>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, <cdouble[:len(self.spec),:len(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 <uintptr_t> &(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

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

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(),
-                        (<FinitePointGroup?>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, <cdouble[:len(self.spec),:len(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 <uintptr_t> &(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 = <qpms_finite_group_t *>0
             self.owns_data = False
 
-    cdef qpms_finite_group_t *rawpointer(self):
-        return self.G
-
 cdef class FinitePointGroupElement:
     '''TODO'''
     cdef readonly FinitePointGroup G

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
         )