From ce905eb0a4e62df310934e1ea01019cf76478e1a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Tue, 20 Aug 2019 12:46:28 +0300 Subject: [PATCH] Expose my incomplete gamma function in Python. More sophisticated GSL error handling Former-commit-id: 4233d01d8aff64b00502c55cb50b54faa5c25ceb --- qpms/__init__.py | 2 +- qpms/ewald.h | 5 +- qpms/qpms_c.pyx | 110 ++++++++++++++++++++++++++++++++++++++++++++ qpms/qpms_cdefs.pxd | 15 ++++++ 4 files changed, 130 insertions(+), 2 deletions(-) diff --git a/qpms/__init__.py b/qpms/__init__.py index 00a1f0f..735021f 100644 --- a/qpms/__init__.py +++ b/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(): diff --git a/qpms/ewald.h b/qpms/ewald.h index 8cc7f4d..486b231 100644 --- a/qpms/ewald.h +++ b/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); diff --git a/qpms/qpms_c.pyx b/qpms/qpms_c.pyx index 2292f08..d539f2b 100644 --- a/qpms/qpms_c.pyx +++ b/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 = 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) diff --git a/qpms/qpms_cdefs.pxd b/qpms/qpms_cdefs.pxd index b331110..2efe72c 100644 --- a/qpms/qpms_cdefs.pxd +++ b/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)