From a18415b78798e53c50c08bdefe83db4a377aa461 Mon Sep 17 00:00:00 2001 From: tlatorre Date: Tue, 5 Mar 2019 13:16:52 -0600 Subject: update quad() to not abort if the matrix is singular --- src/fit.c | 9 +++++---- src/quad.c | 26 +++++++++++++++++++++++++- 2 files changed, 30 insertions(+), 5 deletions(-) diff --git a/src/fit.c b/src/fit.c index ef22d44..8fa3fad 100644 --- a/src/fit.c +++ b/src/fit.c @@ -5248,17 +5248,18 @@ int fit_event2(event *ev, double *xopt, double *fmin, int *id, size_t n, double size_t npeaks; size_t *result; size_t nvertices; + int status; /* Create the minimizer object. */ opt = nlopt_create(NLOPT_LN_BOBYQA, 4+3*n); nlopt_set_min_objective(opt,nlopt_nll2,&fpars); /* Guess the position and t0 of the event using the QUAD fitter. */ - quad(ev,pos,&t0,10000); + status = quad(ev,pos,&t0,10000); - if (t0 < 0 || t0 > GTVALID || NORM(pos) > PSUP_RADIUS) { - /* If the QUAD fitter returns something outside the PSUP or with an - * invalid time we just assume it's at the center. */ + if (status || t0 < 0 || t0 > GTVALID || NORM(pos) > PSUP_RADIUS) { + /* If the QUAD fitter fails or returns something outside the PSUP or + * with an invalid time we just assume it's at the center. */ fprintf(stderr, "quad returned pos = %.2f, %.2f, %.2f t0 = %.2f. Assuming vertex is at the center!\n", pos[0], pos[1], pos[2], t0); pos[0] = 0.0; diff --git a/src/quad.c b/src/quad.c index c869aed..848a12e 100644 --- a/src/quad.c +++ b/src/quad.c @@ -14,9 +14,16 @@ #include #include #include "misc.h" +#include char quad_err[256]; +void my_handler(const char *reason, const char *file, int line, int gsl_errno) +{ + fprintf(stderr, "gsl: %s:%d: %s: %s\n", file, line, "ERROR", reason); + return; +} + /* Returns an approximate vertex position and time using the QUAD fitter. * * The QUAD fitter was originally developed for SNO and estimates the event @@ -83,6 +90,8 @@ int quad(event *ev, double *pos, double *t0, size_t npoints) static size_t p2[MAX_PMTS]; double tmin; size_t max_tries; + gsl_error_handler_t *old_handler; + int status; wavelength0 = 400.0; n_d2o = get_index_snoman_d2o(wavelength0); @@ -156,7 +165,22 @@ int quad(event *ev, double *pos, double *t0, size_t npoints) gsl_vector_view n_view = gsl_vector_view_array(n,3); gsl_linalg_LU_decomp(&m.matrix, p, &s); - gsl_linalg_LU_invert(&m.matrix, p, &minv.matrix); + + /* Ocassionaly the matrix is singular and we can't invert it. Since we + * don't want our program to quit when this happens, we install a new + * gsl error handler, do the matrix inversion, and then restore the old + * handler. If the matrix inversion failed, then we just continue with + * the next quad point. */ + + /* save original handler, install new handler */ + old_handler = gsl_set_error_handler(&my_handler); + + status = gsl_linalg_LU_invert(&m.matrix, p, &minv.matrix); + + /* restore original handler */ + gsl_set_error_handler(old_handler); + + if (status) continue; gsl_blas_dgemv(CblasNoTrans,1.0,&minv.matrix,&k_view.vector,0.0,&g_view.vector); gsl_blas_dgemv(CblasNoTrans,1.0,&minv.matrix,&n_view.vector,0.0,&h_view.vector); -- cgit