Running the attached program nlfit_bug.c triggers undefined behavior,
calculating unreasonable
results. The program "nlfit_debug.c" is adapted from the example program
provided
by GSL's doc/examples/nlfit.c. Their difference lies in only one parameter
40.0
versus 60.0:
1d0
<
12c11,12
< #define TMAX (60.0) /* time variable in [0,TMAX] */ //40.0 instead of
60.0 in nlfit.c
---
> //#define TMAX (40.0) /* time variable in [0,TMAX] */
> #define TMAX (60.0) /* time variable in [0,TMAX] */
The undefined behavior can be better visualized if GCC's Undefined
Behavior Sanitizer is enabled (-fsanitize=undefined). Below is the
execution results:
data: 59.3939 1.03153 0.101317
data: 60 1.14181 0.101239
iter 0: A = 1.0000, lambda = 1.0000, b = 0.0000, cond(J) = inf,
|f(x)| = 101.0200
iter 1: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = 92.8216,
|f(x)| = -nan
iter 2: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 3: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 4: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
*nielsen.c:98:7: runtime error: left shift of 4611686018427387904 by 1
places cannot be represented in type 'long int'*iter 5: A = 3.5110, lambda
= -12.8820, b = 1.2364, cond(J) = nan, |f(x)| = -nan
iter 6: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 7: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 8: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 9: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 10: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 11: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 12: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 13: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 14: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 15: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 16: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 17: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 18: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 19: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 20: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 21: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 22: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 23: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 24: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 25: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 26: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 27: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 28: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 29: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 30: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 31: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 32: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 33: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 34: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 35: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 36: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 37: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 38: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 39: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 40: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 41: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 42: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 43: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 44: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 45: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 46: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 47: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 48: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 49: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 50: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 51: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 52: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 53: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 54: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 55: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 56: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 57: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 58: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 59: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 60: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 61: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 62: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 63: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 64: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 65: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 66: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 67: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 68: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 69: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 70: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 71: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 72: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 73: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 74: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 75: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 76: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 77: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 78: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 79: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 80: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 81: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 82: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 83: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 84: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 85: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 86: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 87: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 88: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 89: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 90: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 91: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 92: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 93: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 94: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 95: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 96: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 97: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 98: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 99: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
iter 100: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan,
|f(x)| = -nan
summary from method 'trust-region/levenberg-marquardt'
number of iterations: 100
function evaluations: 1586
Jacobian evaluations: 2
reason for stopping: small gradient
initial |f(x)| = 101.020000
final |f(x)| = inf
chisq/dof = inf
A = 3.51103 +/- nan
lambda = -12.88201 +/- nan
b = 1.23642 +/- nan
status = exceeded max number of iterations
Preliminary debugging: The undefined behavior occurs in
multifit_nlinear/nielsen.c, line 98, triggered by overflow of a left shift
of a long int.
*nu <<= 1;
BR,
Zhoulai
#include <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_multifit_nlinear.h>
#define N 100 /* number of data points to fit */
#define TMAX (60.0) /* time variable in [0,TMAX] */ //40.0 instead of 60.0 in nlfit.c
struct data {
size_t n;
double * t;
double * y;
};
int
expb_f (const gsl_vector * x, void *data,
gsl_vector * f)
{
size_t n = ((struct data *)data)->n;
double *t = ((struct data *)data)->t;
double *y = ((struct data *)data)->y;
double A = gsl_vector_get (x, 0);
double lambda = gsl_vector_get (x, 1);
double b = gsl_vector_get (x, 2);
size_t i;
for (i = 0; i < n; i++)
{
/* Model Yi = A * exp(-lambda * t_i) + b */
double Yi = A * exp (-lambda * t[i]) + b;
gsl_vector_set (f, i, Yi - y[i]);
}
return GSL_SUCCESS;
}
int
expb_df (const gsl_vector * x, void *data,
gsl_matrix * J)
{
size_t n = ((struct data *)data)->n;
double *t = ((struct data *)data)->t;
double A = gsl_vector_get (x, 0);
double lambda = gsl_vector_get (x, 1);
size_t i;
for (i = 0; i < n; i++)
{
/* Jacobian matrix J(i,j) = dfi / dxj, */
/* where fi = (Yi - yi)/sigma[i], */
/* Yi = A * exp(-lambda * t_i) + b */
/* and the xj are the parameters (A,lambda,b) */
double e = exp(-lambda * t[i]);
gsl_matrix_set (J, i, 0, e);
gsl_matrix_set (J, i, 1, -t[i] * A * e);
gsl_matrix_set (J, i, 2, 1.0);
}
return GSL_SUCCESS;
}
void
callback(const size_t iter, void *params,
const gsl_multifit_nlinear_workspace *w)
{
gsl_vector *f = gsl_multifit_nlinear_residual(w);
gsl_vector *x = gsl_multifit_nlinear_position(w);
double rcond;
/* compute reciprocal condition number of J(x) */
gsl_multifit_nlinear_rcond(&rcond, w);
fprintf(stderr, "iter %2zu: A = %.4f, lambda = %.4f, b = %.4f, cond(J) = %8.4f, |f(x)| = %.4f\n",
iter,
gsl_vector_get(x, 0),
gsl_vector_get(x, 1),
gsl_vector_get(x, 2),
1.0 / rcond,
gsl_blas_dnrm2(f));
}
int
main (void)
{
const gsl_multifit_nlinear_type *T = gsl_multifit_nlinear_trust;
gsl_multifit_nlinear_workspace *w;
gsl_multifit_nlinear_fdf fdf;
gsl_multifit_nlinear_parameters fdf_params =
gsl_multifit_nlinear_default_parameters();
const size_t n = N;
const size_t p = 3;
gsl_vector *f;
gsl_matrix *J;
gsl_matrix *covar = gsl_matrix_alloc (p, p);
double t[N], y[N], weights[N];
struct data d = { n, t, y };
double x_init[3] = { 1.0, 1.0, 0.0 }; /* starting values */
gsl_vector_view x = gsl_vector_view_array (x_init, p);
gsl_vector_view wts = gsl_vector_view_array(weights, n);
gsl_rng * r;
double chisq, chisq0;
int status, info;
size_t i;
const double xtol = 1e-8;
const double gtol = 1e-8;
const double ftol = 0.0;
gsl_rng_env_setup();
r = gsl_rng_alloc(gsl_rng_default);
/* define the function to be minimized */
fdf.f = expb_f;
fdf.df = expb_df; /* set to NULL for finite-difference Jacobian */
fdf.fvv = NULL; /* not using geodesic acceleration */
fdf.n = n;
fdf.p = p;
fdf.params = &d;
/* this is the data to be fitted */
for (i = 0; i < n; i++)
{
double ti = i * TMAX / (n - 1.0);
double yi = 1.0 + 5 * exp (-0.1 * ti);
double si = 0.1 * yi;
double dy = gsl_ran_gaussian(r, si);
t[i] = ti;
y[i] = yi + dy;
weights[i] = 1.0 / (si * si);
printf ("data: %g %g %g\n", ti, y[i], si);
};
/* allocate workspace with default parameters */
w = gsl_multifit_nlinear_alloc (T, &fdf_params, n, p);
/* initialize solver with starting point and weights */
gsl_multifit_nlinear_winit (&x.vector, &wts.vector, &fdf, w);
/* compute initial cost function */
f = gsl_multifit_nlinear_residual(w);
gsl_blas_ddot(f, f, &chisq0);
/* solve the system with a maximum of 100 iterations */
status = gsl_multifit_nlinear_driver(100, xtol, gtol, ftol,
callback, NULL, &info, w);
/* compute covariance of best fit parameters */
J = gsl_multifit_nlinear_jac(w);
gsl_multifit_nlinear_covar (J, 0.0, covar);
/* compute final cost */
gsl_blas_ddot(f, f, &chisq);
#define FIT(i) gsl_vector_get(w->x, i)
#define ERR(i) sqrt(gsl_matrix_get(covar,i,i))
fprintf(stderr, "summary from method '%s/%s'\n",
gsl_multifit_nlinear_name(w),
gsl_multifit_nlinear_trs_name(w));
fprintf(stderr, "number of iterations: %zu\n",
gsl_multifit_nlinear_niter(w));
fprintf(stderr, "function evaluations: %zu\n", fdf.nevalf);
fprintf(stderr, "Jacobian evaluations: %zu\n", fdf.nevaldf);
fprintf(stderr, "reason for stopping: %s\n",
(info == 1) ? "small step size" : "small gradient");
fprintf(stderr, "initial |f(x)| = %f\n", sqrt(chisq0));
fprintf(stderr, "final |f(x)| = %f\n", sqrt(chisq));
{
double dof = n - p;
double c = GSL_MAX_DBL(1, sqrt(chisq / dof));
fprintf(stderr, "chisq/dof = %g\n", chisq / dof);
fprintf (stderr, "A = %.5f +/- %.5f\n", FIT(0), c*ERR(0));
fprintf (stderr, "lambda = %.5f +/- %.5f\n", FIT(1), c*ERR(1));
fprintf (stderr, "b = %.5f +/- %.5f\n", FIT(2), c*ERR(2));
}
fprintf (stderr, "status = %s\n", gsl_strerror (status));
gsl_multifit_nlinear_free (w);
gsl_matrix_free (covar);
gsl_rng_free (r);
return 0;
}
