Expose my incomplete gamma function in Python.

More sophisticated GSL error handling


Former-commit-id: 4233d01d8aff64b00502c55cb50b54faa5c25ceb

qpms/__init__.py

@@ -7,7 +7,7 @@ from sys import platform as __platform
 import warnings as __warnings
 
 try:
-    from .qpms_c import PointGroup, FinitePointGroup, FinitePointGroupElement, Particle, scatsystem_set_nthreads, ScatteringSystem, ScatteringMatrix, pitau
+    from .qpms_c import PointGroup, FinitePointGroup, FinitePointGroupElement, Particle, scatsystem_set_nthreads, ScatteringSystem, ScatteringMatrix, pitau, set_gsl_pythonic_error_handling, pgsl_ignore_error, gamma_inc
 except ImportError as ex:
     if __platform == "linux" or __platform == "linux2":
         if 'LD_LIBRARY_PATH' not in __os.environ.keys():

qpms/ewald.h

@@ -101,7 +101,7 @@ void qpms_ewald3_constants_free(qpms_ewald3_constants_t *);
 
 /// Structure for holding complex-valued result of computation and an error estimate.
 /** Similar to gsl_sf_result, but with complex val. */
-typedef struct { 
+typedef struct qpms_csf_result { 
   complex double val; ///< Calculation result.
   double err; ///< Error estimate.
 } qpms_csf_result;
@@ -136,6 +136,9 @@ static inline complex double clilgamma(complex double z) {
 /** DLMF 8.7.3 (latter expression) for complex second argument */
 int cx_gamma_inc_series_e(double a, complex z, qpms_csf_result * result);
 
+/// Incomplete Gamma function as continued fractions.
+int cx_gamma_inc_CF_e(double a, complex z, qpms_csf_result * result);
+
 /// Incomplete gamma for complex second argument.
 /** if x is (almost) real, it just uses gsl_sf_gamma_inc_e(). */
 int complex_gamma_inc_e(double a, complex double x, qpms_csf_result *result);

qpms/qpms_c.pyx

@@ -14,6 +14,98 @@ from .cycommon cimport make_c_string
 from .cycommon import string_c2py, PointGroupClass
 from .cytmatrices cimport CTMatrix
 from libc.stdlib cimport malloc, free, calloc
+import warnings
+
+
+# Set custom GSL error handler. N.B. this is obviously not thread-safe.
+cdef char *pgsl_err_reason
+cdef char *pgsl_err_file
+cdef int pgsl_err_line
+cdef int pgsl_errno = 0
+cdef int *pgsl_errno_ignorelist = NULL # list of ignored error codes, terminated by zero
+
+# This error handler only sets the variables above
+cdef void pgsl_error_handler(const char *reason, const char *_file, const int line, const int gsl_errno):
+    global pgsl_err_reason, pgsl_err_file, pgsl_err_line, pgsl_errno, pgsl_errno_ignorelist
+    cdef size_t i
+    if(pgsl_errno_ignorelist):
+        i = 0
+        while pgsl_errno_ignorelist[i] != 0:
+            if gsl_errno == pgsl_errno_ignorelist[i]:
+                return
+            i += 1
+    pgsl_err_file = _file
+    pgsl_err_reason = reason
+    pgsl_errno = gsl_errno
+    pgsl_err_line = line
+    return
+
+cdef const int* pgsl_set_ignorelist(const int *new_ignorelist):
+    global pgsl_errno_ignorelist
+    cdef const int *oldlist = pgsl_errno_ignorelist
+    pgsl_errno_ignorelist = new_ignorelist
+    return oldlist
+
+cdef class pgsl_ignore_error():
+    '''Context manager for setting a temporary list of errnos ignored by pgsl_error_handler.
+    Always sets pgsl_error_handler.
+
+    Performs pgsl_check_err() on exit unless 
+    '''
+    cdef const int *ignorelist_old
+    cdef gsl_error_handler_t *old_handler
+    cdef bint final_check
+    cdef object ignorelist_python
+
+    cdef int *ignorelist
+    def __cinit__(self, *ignorelist, **kwargs):
+        self.ignorelist = <int*>calloc((len(ignorelist)+1), sizeof(int))
+        self.ignorelist_python = ignorelist
+        for i in range(len(ignorelist)):
+            self.ignorelist[i] = ignorelist[i]
+        if "final_check" in kwargs.keys() and not kwargs["final_check"]:
+            final_check = True
+        final_check = False
+    
+    def __enter__(self):
+        global pgsl_error_handler
+        self.ignorelist_old = pgsl_set_ignorelist(self.ignorelist)
+        self.old_handler = gsl_set_error_handler(pgsl_error_handler)
+        return
+
+    def __exit__(self, type, value, traceback):
+        global pgsl_errno_ignorelist, pgsl_error_handler
+        pgsl_set_ignorelist(self.ignorelist_old)
+        gsl_set_error_handler(self.old_handler)
+        if self.final_check:
+            pgsl_check_err(retval = None, ignore = self.ignorelist_python)
+
+    def __dealloc__(self):
+        free(self.ignorelist)
+
+def pgsl_check_err(retval = None, ignorelist = None):
+    global pgsl_err_reason, pgsl_err_file, pgsl_err_line, pgsl_errno
+    '''Check for possible errors encountered by pgsl_error_handler.
+    Takes return value of a function as an optional argument, which is now ignored.
+    '''
+    cdef int errno_was
+    if (pgsl_errno != 0):
+        errno_was = pgsl_errno
+        pgsl_errno = 0
+        raise RuntimeError("Error %d in GSL calculation in %s:%d: %s" % (errno_was,
+            string_c2py(pgsl_err_file), pgsl_err_line, string_c2py(pgsl_err_reason)))
+    if (retval is not None and retval != 0 and ignorelist is not None and retval not in ignorelist):
+        warnings.warn("Got non-zero return value %d" % retval)
+    if retval is not None:
+        return retval
+    else:
+        return 0
+
+def set_gsl_pythonic_error_handling():
+    '''
+    Sets pgsl_error_handler as the GSL error handler to avoid crashing.
+    '''
+    gsl_set_error_handler(pgsl_error_handler)
 
 cdef class PointGroup:
     cdef readonly qpms_pointgroup_t G
@@ -499,3 +591,21 @@ def pitau(double theta, qpms_l_t lMax, double csphase = -1):
     cdef double[::1] tau = taua
     qpms_pitau_fill(&leg[0], &pi[0], &tau[0], theta, lMax, csphase)
     return (lega, pia, taua)
+
+def gamma_inc(double a, cdouble x):
+    cdef qpms_csf_result res
+    with pgsl_ignore_error(15): #15 is underflow
+        complex_gamma_inc_e(a, x, &res)
+    return (res.val, res.err)
+
+def gamma_inc_series(double a, cdouble x):
+    cdef qpms_csf_result res
+    with pgsl_ignore_error(15): #15 is underflow
+        cx_gamma_inc_series_e(a, x, &res)
+    return (res.val, res.err)
+
+def gamma_inc_CF(double a, cdouble x):
+    cdef qpms_csf_result res
+    with pgsl_ignore_error(15): #15 is underflow
+        cx_gamma_inc_CF_e(a, x, &res)
+    return (res.val, res.err)

qpms/qpms_cdefs.pxd

@@ -4,6 +4,12 @@ ctypedef double complex cdouble
 
 from libc.stdint cimport *
 
+cdef extern from "gsl/gsl_errno.h":
+    ctypedef void gsl_error_handler_t (const char *reason, const char *file,
+            int line, int gsl_errno)
+    gsl_error_handler_t *gsl_set_error_handler(gsl_error_handler_t *new_handler)
+    gsl_error_handler_t *gsl_set_error_handler_off();
+
 cdef extern from "qpms_types.h":
     cdef struct cart3_t:
         double x
@@ -519,4 +525,13 @@ cdef extern from "scatsystem.h":
             int *target_piv, const qpms_scatsys_t *ss, qpms_iri_t iri)
     cdouble *qpms_scatsys_scatter_solve(cdouble *target_f, const cdouble *a_inc, qpms_ss_LU ludata)
 
+cdef extern from "ewald.h":
+    struct qpms_csf_result:
+        cdouble val
+        double err
+    cdouble lilgamma(double t)
+    cdouble clilgamma(cdouble z)
+    int cx_gamma_inc_series_e(double a, cdouble x, qpms_csf_result *result)
+    int cx_gamma_inc_CF_e(double a, cdouble x, qpms_csf_result *result)
+    int complex_gamma_inc_e(double a, cdouble x, qpms_csf_result *result)