Default gsl error handler invoked

I am writing a python code that needs to calculate a lot of integrals really fast, so instead of using scipy.integrate using a c function, I am using ctypes to calculate all of my integrals in c.

I am importing the c function with this code:

ctypes.CDLL('/usr/lib/i386-linux-gnu/libgslcblas.so', mode=ctypes.RTLD_GLOBAL)
ctypes.CDLL('/usr/lib/i386-linux-gnu/libgsl.so', mode=ctypes.RTLD_GLOBAL)
lib = ctypes.CDLL('/home/aurelien/Desktop/Project/within_opt_model_1/integral_test.so')

In c, I am calculating the integral using:

gsl_integration_workspace *w = gsl_integration_workspace_alloc(200);
double result, err;

gsl_function F;
F.function = &integral;
F.params = &params;

gsl_integration_qags(&F, z, 1, 1e-6, 1e-6, 200, w, &result, &err);
gsl_integration_workspace_free(w);

The problem is that sometimes (~5% of the time), since my integral is slowly convergent, it will kill my program and print:

gsl: qags.c:563: ERROR: interval is divergent, or slowly convergent
Default GSL error handler invoked
Aborted (core dumped)

This problem does not occur if I put 1e-10 instead of 1e-6, but it makes the integral calculation go a lot slower.

I would like a way to do something like a try, catch so that most of the time it uses 1e-6 and goes fast, and when it fails it uses 1e-10.

I have tried doing:

int status = gsl_integration_qags(&F, z, 1, 1e-6, 1e-6, 200, w, &result, &err);
printf("%in", status);

But this only prints 0, since the error aborts before returning any value.

My guess is that I need create my own error handler method, but I don’t know how to do that.

Thank you very much!
(I can show you more code if that’s helpful, I tried to keep everything concise).

@cindex error handling This chapter describes the way that GSL functions report and handle errors. By examining the status information returned by every function you can determine whether it succeeded or failed, and if it failed you can find out what the precise cause of failure was. You can also define your own error handling functions to modify the default behavior of the library. The functions described in this section are declared in the header file @file{gsl_errno.h}. @menu * Error Reporting:: * Error Codes:: * Error Handlers:: * Using GSL error reporting in your own functions:: * Error Reporting Examples:: @end menu @node Error Reporting @section Error Reporting The library follows the thread-safe error reporting conventions of the @sc{posix} Threads library. Functions return a non-zero error code to indicate an error and @code{0} to indicate success. @example int status = gsl_function (…) if (status) @{ /* an error occurred */ ….. /* status value specifies the type of error */ @} @end example The routines report an error whenever they cannot perform the task requested of them. For example, a root-finding function would return a non-zero error code if could not converge to the requested accuracy, or exceeded a limit on the number of iterations. Situations like this are a normal occurrence when using any mathematical library and you should check the return status of the functions that you call. Whenever a routine reports an error the return value specifies the type of error. The return value is analogous to the value of the variable @code{errno} in the C library. The caller can examine the return code and decide what action to take, including ignoring the error if it is not considered serious. In addition to reporting errors by return codes the library also has an error handler function @code{gsl_error}. This function is called by other library functions when they report an error, just before they return to the caller. The default behavior of the error handler is to print a message and abort the program, @example gsl: file.c:67: ERROR: invalid argument supplied by user Default GSL error handler invoked. Aborted @end example The purpose of the @code{gsl_error} handler is to provide a function where a breakpoint can be set that will catch library errors when running under the debugger. It is not intended for use in production programs, which should handle any errors using the return codes. @node Error Codes @section Error Codes @cindex error codes, reserved The error code numbers returned by library functions are defined in the file @file{gsl_errno.h}. They all have the prefix @code{GSL_} and expand to non-zero constant integer values. Error codes above 1024 are reserved for applications, and are not used by the library. Many of the error codes use the same base name as the corresponding error code in the C library. Here are some of the most common error codes, @cindex error codes @deftypevr {Macro} int GSL_EDOM Domain error; used by mathematical functions when an argument value does not fall into the domain over which the function is defined (like EDOM in the C library) @end deftypevr @deftypevr {Macro} int GSL_ERANGE Range error; used by mathematical functions when the result value is not representable because of overflow or underflow (like ERANGE in the C library) @end deftypevr @deftypevr {Macro} int GSL_ENOMEM No memory available. The system cannot allocate more virtual memory because its capacity is full (like ENOMEM in the C library). This error is reported when a GSL routine encounters problems when trying to allocate memory with @code{malloc}. @end deftypevr @deftypevr {Macro} int GSL_EINVAL Invalid argument. This is used to indicate various kinds of problems with passing the wrong argument to a library function (like EINVAL in the C library). @end deftypevr The error codes can be converted into an error message using the function @code{gsl_strerror}. @deftypefun {const char *} gsl_strerror (const int @var{gsl_errno}) This function returns a pointer to a string describing the error code @var{gsl_errno}. For example, @example printf («error: %sn», gsl_strerror (status)); @end example @noindent would print an error message like @code{error: output range error} for a status value of @code{GSL_ERANGE}. @end deftypefun @node Error Handlers @section Error Handlers @cindex Error handlers The default behavior of the GSL error handler is to print a short message and call @code{abort}. When this default is in use programs will stop with a core-dump whenever a library routine reports an error. This is intended as a fail-safe default for programs which do not check the return status of library routines (we don’t encourage you to write programs this way). If you turn off the default error handler it is your responsibility to check the return values of routines and handle them yourself. You can also customize the error behavior by providing a new error handler. For example, an alternative error handler could log all errors to a file, ignore certain error conditions (such as underflows), or start the debugger and attach it to the current process when an error occurs. All GSL error handlers have the type @code{gsl_error_handler_t}, which is defined in @file{gsl_errno.h}, @deftp {Data Type} gsl_error_handler_t This is the type of GSL error handler functions. An error handler will be passed four arguments which specify the reason for the error (a string), the name of the source file in which it occurred (also a string), the line number in that file (an integer) and the error number (an integer). The source file and line number are set at compile time using the @code{__FILE__} and @code{__LINE__} directives in the preprocessor. An error handler function returns type @code{void}. Error handler functions should be defined like this, @example void handler (const char * reason, const char * file, int line, int gsl_errno) @end example @end deftp @comment @noindent To request the use of your own error handler you need to call the function @code{gsl_set_error_handler} which is also declared in @file{gsl_errno.h}, @deftypefun {gsl_error_handler_t *} gsl_set_error_handler (gsl_error_handler_t * @var{new_handler}) This function sets a new error handler, @var{new_handler}, for the GSL library routines. The previous handler is returned (so that you can restore it later). Note that the pointer to a user defined error handler function is stored in a static variable, so there can be only one error handler per program. This function should be not be used in multi-threaded programs except to set up a program-wide error handler from a master thread. The following example shows how to set and restore a new error handler, @example /* save original handler, install new handler */ old_handler = gsl_set_error_handler (&my_handler); /* code uses new handler */ ….. /* restore original handler */ gsl_set_error_handler (old_handler); @end example @noindent To use the default behavior (@code{abort} on error) set the error handler to @code{NULL}, @example old_handler = gsl_set_error_handler (NULL); @end example @end deftypefun @deftypefun {gsl_error_handler_t *} gsl_set_error_handler_off () This function turns off the error handler by defining an error handler which does nothing. This will cause the program to continue after any error, so the return values from any library routines must be checked. This is the recommended behavior for production programs. The previous handler is returned (so that you can restore it later). @end deftypefun The error behavior can be changed for specific applications by recompiling the library with a customized definition of the @code{GSL_ERROR} macro in the file @file{gsl_errno.h}. @node Using GSL error reporting in your own functions @section Using GSL error reporting in your own functions @cindex error handling macros If you are writing numerical functions in a program which also uses GSL code you may find it convenient to adopt the same error reporting conventions as in the library. To report an error you need to call the function @code{gsl_error} with a string describing the error and then return an appropriate error code from @code{gsl_errno.h}, or a special value, such as @code{NaN}. For convenience the file @file{gsl_errno.h} defines two macros which carry out these steps: @deffn {Macro} GSL_ERROR (@var{reason}, @var{gsl_errno}) This macro reports an error using the GSL conventions and returns a status value of @code{gsl_errno}. It expands to the following code fragment, @example gsl_error (reason, __FILE__, __LINE__, gsl_errno); return gsl_errno; @end example @noindent The macro definition in @file{gsl_errno.h} actually wraps the code in a @code{do @{ … @} while (0)} block to prevent possible parsing problems. @end deffn Here is an example of how the macro could be used to report that a routine did not achieve a requested tolerance. To report the error the routine needs to return the error code @code{GSL_ETOL}. @example if (residual > tolerance) @{ GSL_ERROR(«residual exceeds tolerance», GSL_ETOL); @} @end example @deffn {Macro} GSL_ERROR_VAL (@var{reason}, @var{gsl_errno}, @var{value}) This macro is the same as @code{GSL_ERROR} but returns a user-defined value of @var{value} instead of an error code. It can be used for mathematical functions that return a floating point value. @end deffn The following example shows how to return a @code{NaN} at a mathematical singularity using the @code{GSL_ERROR_VAL} macro, @example if (x == 0) @{ GSL_ERROR_VAL(«argument lies on singularity», GSL_ERANGE, GSL_NAN); @} @end example @node Error Reporting Examples @section Examples Here is an example of some code which checks the return value of a function where an error might be reported, @example #include <stdio.h> #include <gsl/gsl_errno.h> #include <gsl/gsl_fft_complex.h> … int status; size_t n = 37; gsl_set_error_handler_off(); status = gsl_fft_complex_radix2_forward (data, stride, n); if (status) @{ if (status == GSL_EINVAL) @{ fprintf (stderr, «invalid argument, n=%dn», n); @} else @{ fprintf (stderr, «failed, gsl_errno=%dn», status); @} exit (-1); @} … @end example @comment @noindent The function @code{gsl_fft_complex_radix2} only accepts integer lengths which are a power of two. If the variable @code{n} is not a power of two then the call to the library function will return @code{GSL_EINVAL}, indicating that the length argument is invalid. The function call to @code{gsl_set_error_handler_off} stops the default error handler from aborting the program. The @code{else} clause catches any other possible errors.

Hi,
I met the same problem when I read someother’s codes. But unfortunately, I can’t solve it.
The function is:
CvMat* ransac_xform( struct feature* features, int n, int mtype,
ransac_xform_fn xform_fn, int m, double p_badxform,
ransac_err_fn err_fn, double err_tol,
struct feature*** inliers, int* n_in )
{
struct feature** matched, ** sample, ** consensus, ** consensus_max = NULL;
struct ransac_data* rdata;
CvPoint2D64f* pts, * mpts;
CvMat* M = NULL;
gsl_rng* rng;
double p, in_frac = RANSAC_INLIER_FRAC_EST;
int i, nm, in, in_min, in_max = 0, k = 0;

nm = get_matched_features( features, n, mtype, &matched );
if( nm < m )
{
    fprintf( stderr, "Warning: not enough matches to compute xform, %s" 
        " line %dn", __FILE__, __LINE__ );
    goto end;
}

/* initialize random number generator */
rng = gsl_rng_alloc( gsl_rng_mt19937 );
gsl_rng_set( rng, time(NULL) );

in_min = calc_min_inliers( nm, m, RANSAC_PROB_BAD_SUPP, p_badxform );
p = pow( 1.0 - pow( in_frac, m ), k );
i = 0;
while( p > p_badxform )
{
    sample = draw_ransac_sample( matched, nm, m, rng );
    extract_corresp_pts( sample, m, mtype, &pts, &mpts );
    M = xform_fn( pts, mpts, m );
    if( ! M )
        goto iteration_end;
    in = find_consensus( matched, nm, mtype, M, err_fn, err_tol, &consensus);
    if( in > in_max )
    {
        if( consensus_max )
            free( consensus_max );
        consensus_max = consensus;
        in_max = in;
        in_frac = (double)in_max / nm;
    }
    else
        free( consensus );
    cvReleaseMat( &M );

iteration_end:
release_mem( pts, mpts, sample );
p = pow( 1.0 — pow( in_frac, m ), ++k );
}

/* calculate final transform based on best consensus set */
if( in_max >= in_min )
{
    extract_corresp_pts( consensus_max, in_max, mtype, &pts, &mpts );
    M = xform_fn( pts, mpts, in_max );
    in = find_consensus( matched, nm, mtype, M, err_fn, err_tol, &consensus);
    cvReleaseMat( &M );
    release_mem( pts, mpts, consensus_max );
    extract_corresp_pts( consensus, in, mtype, &pts, &mpts );
    M = xform_fn( pts, mpts, in );
    if( inliers )
    {
        *inliers = consensus;
        consensus = NULL;
    }
    if( n_in )
        *n_in = in;
    release_mem( pts, mpts, consensus );
}
else if( consensus_max )
{
    if( inliers )
        *inliers = NULL;
    if( n_in )
        *n_in = 0;
    free( consensus_max );
}

gsl_rng_free( rng );

end:
for( i = 0; i < nm; i++ )
{
rdata = feat_ransac_data( matched[i] );
matched[i]->feature_data = rdata->orig_feat_data;
free( rdata );
}
free( matched );
return M;
}

I didn’t find the function gsl_rng_alloc that was repeatlly called, and the function ransac_xform is only called by the main function:

int main( int argc, char** argv )
{
IplImage* img1, * img2, * stacked;
struct feature* feat1, * feat2, * feat;
struct feature** nbrs;
struct kd_node* kd_root;
CvPoint pt1, pt2;
double d0, d1;
int n1, n2, k, i, m = 0;
img1 = cvLoadImage( img1_file, 1 );
if( ! img1 )
fatal_error( «unable to load image from %s», img1_file );
img2 = cvLoadImage( img2_file, 1 );
if( ! img2 )
fatal_error( «unable to load image from %s», img2_file );
stacked = stack_imgs( img1, img2 );

fprintf( stderr, "Finding features in %s...n", img1_file );
n1 = sift_features( img1, &feat1 );
fprintf( stderr, "Finding features in %s...n", img2_file );
n2 = sift_features( img2, &feat2 );
kd_root = kdtree_build( feat2, n2 );
for( i = 0; i < n1; i++ )
{
    feat = feat1 + i;
    k = kdtree_bbf_knn( kd_root, feat, 2, &nbrs, KDTREE_BBF_MAX_NN_CHKS );
    if( k == 2 )
    {
        d0 = descr_dist_sq( feat, nbrs[0] );
        d1 = descr_dist_sq( feat, nbrs[1] );
        if( d0 < d1 * NN_SQ_DIST_RATIO_THR )
        {
            pt1 = cvPoint( cvRound( feat->x ), cvRound( feat->y ) );
            pt2 = cvPoint( cvRound( nbrs[0]->x ), cvRound( nbrs[0]->y ) );
            pt2.y += img1->height;
            cvLine( stacked, pt1, pt2, CV_RGB(255,0,255), 1, 8, 0 );
            m++;
            feat1[i].fwd_match = nbrs[0];
        }
    }
    free( nbrs );
}

fprintf( stderr, "Found %d total matchesn", m );
cvNamedWindow( "Matches", 1 );
cvShowImage( "Matches", stacked );
cvWaitKey( 0 );


/* 
UNCOMMENT BELOW TO SEE HOW RANSAC FUNCTION WORKS

Note that this line above:

feat1[i].fwd_match = nbrs[0];

is important for the RANSAC function to work.
*/

{
    CvMat* H;
    H = ransac_xform( feat1, n1, FEATURE_FWD_MATCH, lsq_homog, 4, 0.01,
        homog_xfer_err, 3.0, NULL, NULL );

    if( H )
    {
        IplImage* xformed;
        xformed = cvCreateImage( cvGetSize( img2 ), IPL_DEPTH_8U, 3 );
        cvWarpPerspective( img1, xformed, H, 
            CV_INTER_LINEAR + CV_WARP_FILL_OUTLIERS,
            cvScalarAll( 0 ) );
        cvNamedWindow( "Xformed", 1 );
        cvShowImage( "Xformed", xformed );
        cvWaitKey( 0 );
        cvReleaseImage( &xformed );
        cvReleaseMat( &H );
    }
}


cvReleaseImage( &stacked );
cvReleaseImage( &img1 );
cvReleaseImage( &img2 );
kdtree_release( kd_root );
free( feat1 );



free( feat2 );
return 0;

}

Any help very much appreciated.
Thanks!

Login before adding your answer.

I hope someone could say me where can i find the meaning of:

gsl: sinint.c:359: ERROR: domain error
Default GSL error handler invoked.
Aborted

or where i can find it.
I don’t understand the meaning. I tried to find the error but it seems
inside an «if ()», and it’s not a GSL instruction. The program was ok
since i addicted at the function:

if (fabs (j — j_esterno) == 1)
{
distanza_collinear = dist_adiacenti_collinear ;
}
else
{
distanza_collinear = (dist_adiacenti_collinear / fabs (j —
j_esterno)) + (lunghezza_dipolo / fabs (j_esterno — j — 1)) ;
}

Now the function is:

void onda_vert_piano_gp (gsl_matrix_complex *matrice_acc_dipoli,
gsl_vector *elementi_guasti_x, gsl_vector *elementi_guasti_y, int
guasti, int N, int M, int dist_adiacenti_side_by_side, int
dist_adiacenti_collinear, int lunghezza_dipolo, char *parassita)
{
int i, j, i_esterno, j_esterno, guasto = 0 ;
double distanza_collinear, distanza_side_by_side ;
gsl_complex accoppiamento, collinear_x, side_y ;

for (i_esterno = 0 ; i_esterno < M ; i_esterno++)
{
for (j_esterno = 0 ; j_esterno < N ; j_esterno++)
{
for (i = 0 ; i < M ; i++)
{
/* Lungo x sono collinear. */
// distanza_collinear = dist_adiacenti_collinear * fabs
(i_esterno — i) ;
distanza_side_by_side =
dist_adiacenti_side_by_side / fabs (i_esterno — i) ;
for (j = 0 ; j < N ; j++)
{
if ((i_esterno == i) && (j_esterno == j))
{
guasto = ricerca_guasto_piano
(elementi_guasti_x, elementi_guasti_y, guasti, i, j, i_esterno,
j_esterno) ;
if (guasto == 1)
{
if (strcmp (parassita, «ca») == 0)
{
GSL_SET_COMPLEX (&accoppiamento, 10e30,
0.0) ; /* Caso di ca. */
}
else
{
GSL_SET_COMPLEX (&accoppiamento, 0.0,
0.0) ; /* Caso di cc. */
}
guasto = 0 ;
}
else
{
auto_impedenza (&accoppiamento,
lunghezza_dipolo) ;
}
}
else
{
/* Lungo y sono side-by-side. */
if (fabs (j — j_esterno) == 1.0)
{
distanza_collinear = dist_adiacenti_collinear ;
}
else
{
distanza_collinear = (dist_adiacenti_collinear / fabs (j —
j_esterno)) + (lunghezza_dipolo / fabs (j_esterno — j — 1)) ;
}
guasto_piano_d_gp (distanza_collinear,
distanza_side_by_side, lunghezza_dipolo, &side_y, &collinear_x,
&accoppiamento) ;
}
gsl_matrix_complex_set (matrice_acc_dipoli, N *
i_esterno + j_esterno, N * i + j, accoppiamento) ;
}
}
}
}
}

Help Please!!!