cython wrapper for C T-matrix interpolator.

Former-commit-id: 3c2c27b504d354513c0dd10b767f265fa6bd19c1

qpms/qpms_c.pyx

@@ -1067,13 +1067,39 @@ cdef class TMatrixInterpolator:
     '''
     Wrapper over the qpms_tmatrix_interpolator_t structure.
     '''
-    def __cinit__(self, filename, *args, **kwargs):
+    #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'''
-        pass
+        self.spec = bspec
+        cdef char * cpath = make_c_string(filename)
+        if QPMS_SUCCESS != qpms_load_scuff_tmatrix(cpath, self.spec.rawpointer(),
+                &(self.nfreqs), &(self.freqs), &(self.freqs_su),
+                &(self.tmatrices_array), &(self.tmdata)):
+            raise IOError("Could not read T-matrix from %s" % filename)
+        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, freq):
         '''Returns a TMatrix instance, corresponding to a given frequency.'''
-        pass
+        # This is a bit stupid, I should rethink the CTMatrix constuctors
-    pass
+        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)
 
 cdef class CTMatrix: # N.B. there is another type called TMatrix in tmatrices.py!
     '''

qpms/qpms_cdefs.pxd

@@ -164,6 +164,8 @@ cdef extern from "symmetries.h":
     cdouble *qpms_zrot_uvswi_dense(cdouble *target, const qpms_vswf_set_spec_t *bspec, double phi)
     cdouble *qpms_zrot_rational_uvswi_dense(cdouble *target, const qpms_vswf_set_spec_t *bspec, int N, int w)
     cdouble *qpms_irot3_uvswfi_dense(cdouble *target, const qpms_vswf_set_spec_t *bspec, qpms_irot3_t transf)
+    size_t qpms_zero_roundoff_clean(double *arr, size_t nmemb, double atol)
+    size_t qpms_czero_roundoff_clean(cdouble *arr, size_t nmemb, double atol)
 
 #cdef extern from "numpy/arrayobject.h":
 #    cdef enum NPY_TYPES:
@@ -212,6 +214,13 @@ cdef extern from "translations.h":
             char *phi_data, np.npy_intp *phi_shape, np.npy_intp *phi_strides,
             char *r_ge_d_data, np.npy_intp *phi_shape, np.npy_intp *phi_strides) nogil
 
+cdef extern from "gsl/gsl_interp.h":
+    struct gsl_interp_type:
+        pass
+    gsl_interp_type gsl_interp_linear
+    gsl_interp_type gsl_interp_cspline
+    # ^^^ These are probably the only relevant ones.
+
 cdef extern from "scatsystem.h":
     struct qpms_tmatrix_t:
         qpms_vswf_set_spec_t *spec
@@ -225,6 +234,11 @@ cdef extern from "scatsystem.h":
         qpms_ss_tmi_t tmatrix_id
     struct qpms_tmatrix_interpolator_t:
         const qpms_vswf_set_spec_t *bspec
+    void qpms_tmatrix_interpolator_free(qpms_tmatrix_interpolator_t *interp)
+    qpms_tmatrix_t *qpms_tmatrix_interpolator_eval(const qpms_tmatrix_interpolator_t *interp, double freq)
+    qpms_tmatrix_interpolator_t *qpms_tmatrix_interpolator_create(size_t n, double *freqs, 
+            const qpms_tmatrix_t *tmatrices_array, const gsl_interp_type *iptype)
+    void qpms_tmatrix_free(qpms_tmatrix_t *tmatrix)
     struct qpms_scatsys_t:
         qpms_tmatrix_t **tm
         qpms_ss_tmi_t tm_count
@@ -237,8 +251,27 @@ cdef extern from "scatsystem.h":
     qpms_scatsys_t *qpms_scatsys_load(char *path) #NI
     qpms_tmatrix_isclose(const qpms_tmatrix_t *A, const qpms_tmatrix_t *B,
                 const double rtol, const double atol)
-        
+    qpms_errno_t qpms_symmetrise_tmdata_irot3arr(
-        
+            cdouble *tmdata, const size_t tmcount,
-
+            const qpms_vswf_set_spec_t *bspec,
-
+            size_t n_symops,
+            const qpms_irot3_t *symops
+            )
+    qpms_errno_t qpms_symmetrise_tmdata_finite_group(
+            cdouble *tmdata, const size_t tmcount,
+            const qpms_vswf_set_spec_t *bspec,
+            const qpms_finite_group_t *pointgroup
+            )
+    qpms_tmatrix_t *qpms_tmatrix_symmetrise_irot3arr_inplace(
+            qpms_tmatrix_t *T,
+            size_t n_symops,
+            const qpms_irot3_t *symops
+            )
+    qpms_tmatrix_t *qpms_tmatrix_symmetrise_finite_group_inplace(
+            qpms_tmatrix_t *T,
+            const qpms_finite_group_t *pointgroup
+            )
+    qpms_errno_t qpms_load_scuff_tmatrix(const char *path, const qpms_vswf_set_spec_t *bspec,
+            size_t *n, double **freqs, double **freqs_su, qpms_tmatrix_t **tmatrices_array,
+            cdouble **tmdata)
 

qpms/scatsystem.h

@@ -269,6 +269,7 @@ typedef struct qpms_tmatrix_interpolator_t {
 void qpms_tmatrix_interpolator_free(qpms_tmatrix_interpolator_t *interp);
 
 /// Evaluate a T-matrix interpolated value.
+/** The result is to be freed using qpms_tmatrix_free().*/
 qpms_tmatrix_t *qpms_tmatrix_interpolator_eval(const qpms_tmatrix_interpolator_t *interp, double freq);
 
 /// Create a T-matrix interpolator from frequency and T-matrix arrays.