Finer error handling.

Former-commit-id: f30b79d2cc321fc06030374176f5013ff179ffc8
This commit is contained in:
Marek Nečada 2020-03-23 11:55:26 +02:00
parent 2798fcce49
commit 3ebc1af946
7 changed files with 45 additions and 13 deletions

View File

@ -59,6 +59,11 @@ qpms_errno_t qpms_sph_bessel_realx_fill(qpms_bessel_t typ, qpms_l_t lmax, double
return retval; return retval;
break; break;
case QPMS_BESSEL_SINGULAR: //FIXME: is this precise enough? Would it be better to do it one-by-one? 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); retval = gsl_sf_bessel_yl_array(lmax,x,tmparr);
for (int l = 0; l <= lmax; ++l) result_array[l] = tmparr[l]; for (int l = 0; l <= lmax; ++l) result_array[l] = tmparr[l];
return retval; 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]; for (int l = 0; l <= lmax; ++l) result_array[l] = tmparr[l];
if(retval) return retval; if(retval) return retval;
retval = gsl_sf_bessel_yl_array(lmax, x, tmparr); 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) if (typ==QPMS_HANKEL_PLUS)
for (int l = 0; l <= lmax; ++l) result_array[l] += I * tmparr[l]; for (int l = 0; l <= lmax; ++l) result_array[l] += I * tmparr[l];
else else
@ -76,10 +86,9 @@ qpms_errno_t qpms_sph_bessel_realx_fill(qpms_bessel_t typ, qpms_l_t lmax, double
return retval; return retval;
break; break;
default: default:
abort(); QPMS_INVALID_ENUM(typ);
//return GSL_EDOM;
} }
assert(0); QPMS_WTF;
} }
// TODO DOC // TODO DOC
@ -119,7 +128,7 @@ qpms_errno_t qpms_sph_bessel_fill(qpms_bessel_t typ, qpms_l_t lmax, complex doub
} }
break; break;
default: default:
QPMS_WTF; QPMS_INVALID_ENUM(typ);
} }
// TODO check for underflows? (nz != 0) // TODO check for underflows? (nz != 0)
if (ierr == 0 || ierr == 3) { if (ierr == 0 || ierr == 3) {

View File

@ -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_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 {\ #define QPMS_UNTESTED {\
static _Bool already_bitched = 0; \ static _Bool already_bitched = 0; \
if (QPMS_UNLIKELY(!already_bitched)) {\ 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__) #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) { \ #define QPMS_ENSURE_SUCCESS(x) { \
int errorcode = (x); \ int errorcode = (x); /* evaluate x only once */ \
if(QPMS_UNLIKELY(errorcode)) \ if(QPMS_UNLIKELY(errorcode)) \
qpms_pr_error_at_flf(__FILE__,__LINE__,__func__,"Unexpected error (%d)", 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__); \ 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__); } #define QPMS_ENSURE(x, msg, ...) {if(QPMS_UNLIKELY(!(x))) qpms_pr_error_at_flf(__FILE__,__LINE__,__func__,msg, ##__VA_ARGS__); }

View File

@ -78,13 +78,14 @@ typedef enum {
typedef unsigned long long qpms_uvswfi_t; typedef unsigned long long qpms_uvswfi_t;
/// Error codes / return values for certain numerical functions. /// 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 { typedef enum {
QPMS_SUCCESS = 0, ///< Success. 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_ENOMEM = 8, ///< Out of memory.
QPMS_ESING = 21, ///< Apparent singularity detected.
QPMS_NAN_ENCOUNTERED = 1024 ///< NaN value encountered in data processing. QPMS_NAN_ENCOUNTERED = 1024 ///< NaN value encountered in data processing.
} qpms_errno_t; } qpms_errno_t;

View File

@ -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, 0, 1),
qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, +1, 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, .capacity=0,
.norm = ss->c->normalisation, .norm = ss->c->normalisation,
}; };

View File

@ -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) { csph_t kr, qpms_bessel_t btyp, const qpms_normalisation_t norm) {
assert(lMax >= 1); assert(lMax >= 1);
complex double *bessel = malloc((lMax+1)*sizeof(complex double)); 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)); 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 complex double const *pbes = bessel + 1; // starting from l = 1
double const *pleg = pt.leg; double const *pleg = pt.leg;
@ -282,7 +283,7 @@ qpms_errno_t qpms_vswf_fill_csph(csphvec_t *const longtarget,
} }
free(bessel); free(bessel);
qpms_pitau_free(pt); qpms_pitau_free(pt);
return QPMS_SUCCESS; return bessel_retval;
} }
qpms_errno_t qpms_vswf_fill(csphvec_t *const longtarget, qpms_errno_t qpms_vswf_fill(csphvec_t *const longtarget,

View File

@ -177,6 +177,11 @@ csphvec_t qpms_vswf_L00(
* *
* Does not evaluate the zeroth-order wave \f$ \mathbf{L}_0^0 \f$. * Does not evaluate the zeroth-order wave \f$ \mathbf{L}_0^0 \f$.
* If you need that, use qpms_vswf_L00(). * 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( qpms_errno_t qpms_vswf_fill(
csphvec_t *resultL, //< Target array for longitudinal VSWFs. csphvec_t *resultL, //< Target array for longitudinal VSWFs.

View File

@ -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 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 from qpms.symmetries import point_group_info
import numpy as np import numpy as np
eh = eV/hbar eh = eV/hbar
@ -15,6 +16,8 @@ dbgmsg_enable(2)
part_radius = 80e-9 part_radius = 80e-9
p = 1580e-9 p = 1580e-9
set_gsl_pythonic_error_handling()
sym = FinitePointGroup(point_group_info['D4h']) sym = FinitePointGroup(point_group_info['D4h'])
bspec1 = BaseSpec(lMax=3) bspec1 = BaseSpec(lMax=3)
medium=EpsMuGenerator(EpsMu(1.52**2)) medium=EpsMuGenerator(EpsMu(1.52**2))
@ -36,7 +39,7 @@ for i in range(ss.fecv_size):
E = ssw.scattered_E(fvc, points) E = ssw.scattered_E(fvc, points)
E_alt = ssw.scattered_E(fvc, points,alt=True) E_alt = ssw.scattered_E(fvc, points,alt=True)
diff = abs(E-E_alt) 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 fail = reldiffavg > 1e-3
fails += fail fails += fail