From 3ebc1af946119240aace3f956011ea8b0278d9b1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Mon, 23 Mar 2020 11:55:26 +0200 Subject: [PATCH] Finer error handling. Former-commit-id: f30b79d2cc321fc06030374176f5013ff179ffc8 --- qpms/bessel.c | 17 +++++++++++++---- qpms/qpms_error.h | 17 +++++++++++++++-- qpms/qpms_types.h | 7 ++++--- qpms/scatsystem.c | 2 +- qpms/vswf.c | 5 +++-- qpms/vswf.h | 5 +++++ tests/test_vswf_transop_consistency.py | 5 ++++- 7 files changed, 45 insertions(+), 13 deletions(-) diff --git a/qpms/bessel.c b/qpms/bessel.c index 1d11805..4b3edb1 100644 --- a/qpms/bessel.c +++ b/qpms/bessel.c @@ -59,6 +59,11 @@ qpms_errno_t qpms_sph_bessel_realx_fill(qpms_bessel_t typ, qpms_l_t lmax, double return retval; break; case QPMS_BESSEL_SINGULAR: //FIXME: is this precise enough? Would it be better to do it one-by-one? + if(QPMS_UNLIKELY(x == 0)) { + for (int l = 0; l <= lmax; ++l) + result_array[l] = NAN; + return QPMS_ESING; // GSL would have returned GSL_EDOM without setting NANs. + } retval = gsl_sf_bessel_yl_array(lmax,x,tmparr); for (int l = 0; l <= lmax; ++l) result_array[l] = tmparr[l]; return retval; @@ -69,6 +74,11 @@ qpms_errno_t qpms_sph_bessel_realx_fill(qpms_bessel_t typ, qpms_l_t lmax, double for (int l = 0; l <= lmax; ++l) result_array[l] = tmparr[l]; if(retval) return retval; retval = gsl_sf_bessel_yl_array(lmax, x, tmparr); + if(QPMS_UNLIKELY(x == 0)) { + for (int l = 0; l <= lmax; ++l) + result_array[l] += I * NAN; + return QPMS_ESING; // GSL would have returned GSL_EDOM without setting NANs. + } if (typ==QPMS_HANKEL_PLUS) for (int l = 0; l <= lmax; ++l) result_array[l] += I * tmparr[l]; else @@ -76,10 +86,9 @@ qpms_errno_t qpms_sph_bessel_realx_fill(qpms_bessel_t typ, qpms_l_t lmax, double return retval; break; default: - abort(); - //return GSL_EDOM; + QPMS_INVALID_ENUM(typ); } - assert(0); + QPMS_WTF; } // TODO DOC @@ -119,7 +128,7 @@ qpms_errno_t qpms_sph_bessel_fill(qpms_bessel_t typ, qpms_l_t lmax, complex doub } break; default: - QPMS_WTF; + QPMS_INVALID_ENUM(typ); } // TODO check for underflows? (nz != 0) if (ierr == 0 || ierr == 3) { diff --git a/qpms/qpms_error.h b/qpms/qpms_error.h index e40e105..14b68b9 100644 --- a/qpms/qpms_error.h +++ b/qpms/qpms_error.h @@ -83,6 +83,8 @@ qpms_dbgmsg_flags qpms_dbgmsg_enable(qpms_dbgmsg_flags types); #define QPMS_WTF qpms_pr_error_at_flf(__FILE__,__LINE__,__func__,"Unexpected error.") +#define QPMS_INVALID_ENUM(x) qpms_pr_error_at_flf(__FILE__,__LINE__,__func__,"Invalid enum value (" #x " == %d)", (int) (x)) + #define QPMS_UNTESTED {\ static _Bool already_bitched = 0; \ if (QPMS_UNLIKELY(!already_bitched)) {\ @@ -93,9 +95,9 @@ qpms_dbgmsg_flags qpms_dbgmsg_enable(qpms_dbgmsg_flags types); #define QPMS_PR_ERROR(msg, ...) qpms_pr_error_at_flf(__FILE__,__LINE__,__func__,msg, ##__VA_ARGS__) -// TODO replace this macro with an inline function? +/// Raises an error if \a x is not zero. #define QPMS_ENSURE_SUCCESS(x) { \ - int errorcode = (x); \ + int errorcode = (x); /* evaluate x only once */ \ if(QPMS_UNLIKELY(errorcode)) \ qpms_pr_error_at_flf(__FILE__,__LINE__,__func__,"Unexpected error (%d)", errorcode); \ } @@ -107,6 +109,17 @@ qpms_dbgmsg_flags qpms_dbgmsg_enable(qpms_dbgmsg_flags types); qpms_pr_error_at_flf(__FILE__,__LINE__,__func__,msg, ##__VA_ARGS__); \ } +/// Raises an error if \a x is not zero or one of the values listed in arguments. +#define QPMS_ENSURE_SUCCESS_OR(x, ...) { \ + int errorcode = (x); /* evaluate x only once */ \ + static const int allowed_errorcodes[] = {0, ##__VA_ARGS__};\ + static const int n_allowed = sizeof(allowed_errorcodes) / sizeof(int); \ + int i; \ + for(i = 0; i < n_allowed; ++i) \ + if (errorcode == allowed_errorcodes[i]) break; \ + if (QPMS_UNLIKELY(i >= n_allowed)) \ + qpms_pr_error_at_flf(__FILE__,__LINE__,__func__,"Unexpected error (%d)", errorcode); \ +} #define QPMS_ENSURE(x, msg, ...) {if(QPMS_UNLIKELY(!(x))) qpms_pr_error_at_flf(__FILE__,__LINE__,__func__,msg, ##__VA_ARGS__); } diff --git a/qpms/qpms_types.h b/qpms/qpms_types.h index 8b169b2..1d8f7dd 100644 --- a/qpms/qpms_types.h +++ b/qpms/qpms_types.h @@ -78,13 +78,14 @@ typedef enum { typedef unsigned long long qpms_uvswfi_t; /// Error codes / return values for certain numerical functions. -/** These are de facto a subset of the GSL error codes. */ +/** Those with values between -2 and 32 are subset of the GSL error codes. */ typedef enum { QPMS_SUCCESS = 0, ///< Success. - QPMS_ERROR = 1, ///< Unspecified error. + QPMS_ERROR = -1, ///< Unspecified error. + QPMS_FAILURE = QPMS_ERROR, QPMS_ENOMEM = 8, ///< Out of memory. + QPMS_ESING = 21, ///< Apparent singularity detected. QPMS_NAN_ENCOUNTERED = 1024 ///< NaN value encountered in data processing. - } qpms_errno_t; diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c index 17a5e47..bc6d60b 100644 --- a/qpms/scatsystem.c +++ b/qpms/scatsystem.c @@ -2039,7 +2039,7 @@ ccart3_t qpms_scatsys_scattered_E__alt(const qpms_scatsys_t *ss, qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, 0, 1), qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, +1, 1), }, - .lMax=1, .lMax_M=1, .lMax_N=1, .lMax_L=-1, + .lMax=1, .lMax_M=0, .lMax_N=1, .lMax_L=-1, .capacity=0, .norm = ss->c->normalisation, }; diff --git a/qpms/vswf.c b/qpms/vswf.c index 19a87d0..f827364 100644 --- a/qpms/vswf.c +++ b/qpms/vswf.c @@ -231,7 +231,8 @@ qpms_errno_t qpms_vswf_fill_csph(csphvec_t *const longtarget, csph_t kr, qpms_bessel_t btyp, const qpms_normalisation_t norm) { assert(lMax >= 1); complex double *bessel = malloc((lMax+1)*sizeof(complex double)); - QPMS_ENSURE_SUCCESS(qpms_sph_bessel_fill(btyp, lMax, kr.r, bessel)); + qpms_errno_t bessel_retval = qpms_sph_bessel_fill(btyp, lMax, kr.r, bessel); + QPMS_ENSURE_SUCCESS_OR(bessel_retval, QPMS_ESING); qpms_pitau_t pt = qpms_pitau_get(kr.theta, lMax, qpms_normalisation_t_csphase(norm)); complex double const *pbes = bessel + 1; // starting from l = 1 double const *pleg = pt.leg; @@ -282,7 +283,7 @@ qpms_errno_t qpms_vswf_fill_csph(csphvec_t *const longtarget, } free(bessel); qpms_pitau_free(pt); - return QPMS_SUCCESS; + return bessel_retval; } qpms_errno_t qpms_vswf_fill(csphvec_t *const longtarget, diff --git a/qpms/vswf.h b/qpms/vswf.h index 33dae0e..d2d01ac 100644 --- a/qpms/vswf.h +++ b/qpms/vswf.h @@ -177,6 +177,11 @@ csphvec_t qpms_vswf_L00( * * Does not evaluate the zeroth-order wave \f$ \mathbf{L}_0^0 \f$. * If you need that, use qpms_vswf_L00(). + * + * If \a kdrj.r == 0 and \a btyp != QPMS_BESSEL_REGULAR, returns QPMS_ESING and fills + * targets with NaNs. + * + * \return 0 (QPMS_SUCCESS) on success, otherwise. */ qpms_errno_t qpms_vswf_fill( csphvec_t *resultL, //< Target array for longitudinal VSWFs. diff --git a/tests/test_vswf_transop_consistency.py b/tests/test_vswf_transop_consistency.py index d6ea7ba..2674109 100755 --- a/tests/test_vswf_transop_consistency.py +++ b/tests/test_vswf_transop_consistency.py @@ -6,6 +6,7 @@ waves at origin """ from qpms import Particle, CTMatrix, lorentz_drude, EpsMuGenerator, TMatrixGenerator, BaseSpec, FinitePointGroup, ScatteringSystem, TMatrixInterpolator, EpsMu, dbgmsg_enable, dbgmsg_disable, dbgmsg_active, BesselType,eV, hbar +from qpms.qpms_c import set_gsl_pythonic_error_handling from qpms.symmetries import point_group_info import numpy as np eh = eV/hbar @@ -15,6 +16,8 @@ dbgmsg_enable(2) part_radius = 80e-9 p = 1580e-9 +set_gsl_pythonic_error_handling() + sym = FinitePointGroup(point_group_info['D4h']) bspec1 = BaseSpec(lMax=3) medium=EpsMuGenerator(EpsMu(1.52**2)) @@ -36,7 +39,7 @@ for i in range(ss.fecv_size): E = ssw.scattered_E(fvc, points) E_alt = ssw.scattered_E(fvc, points,alt=True) diff = abs(E-E_alt) - reldiffavg = np.average(diff/(abs(E)+abs(E_alt))) + reldiffavg = np.nanmean(diff/(abs(E)+abs(E_alt))) fail = reldiffavg > 1e-3 fails += fail