|
Packit |
67cb25 |
#include <stdlib.h>
|
|
Packit |
67cb25 |
#include <stdio.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_rng.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_randist.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_matrix.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_vector.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_blas.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_multifit_nlinear.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#define N 100 /* number of data points to fit */
|
|
Packit |
67cb25 |
#define TMAX (40.0) /* time variable in [0,TMAX] */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
struct data {
|
|
Packit |
67cb25 |
size_t n;
|
|
Packit |
67cb25 |
double * t;
|
|
Packit |
67cb25 |
double * y;
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
expb_f (const gsl_vector * x, void *data,
|
|
Packit |
67cb25 |
gsl_vector * f)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t n = ((struct data *)data)->n;
|
|
Packit |
67cb25 |
double *t = ((struct data *)data)->t;
|
|
Packit |
67cb25 |
double *y = ((struct data *)data)->y;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double A = gsl_vector_get (x, 0);
|
|
Packit |
67cb25 |
double lambda = gsl_vector_get (x, 1);
|
|
Packit |
67cb25 |
double b = gsl_vector_get (x, 2);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
size_t i;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < n; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* Model Yi = A * exp(-lambda * t_i) + b */
|
|
Packit |
67cb25 |
double Yi = A * exp (-lambda * t[i]) + b;
|
|
Packit |
67cb25 |
gsl_vector_set (f, i, Yi - y[i]);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
expb_df (const gsl_vector * x, void *data,
|
|
Packit |
67cb25 |
gsl_matrix * J)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t n = ((struct data *)data)->n;
|
|
Packit |
67cb25 |
double *t = ((struct data *)data)->t;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double A = gsl_vector_get (x, 0);
|
|
Packit |
67cb25 |
double lambda = gsl_vector_get (x, 1);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
size_t i;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < n; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* Jacobian matrix J(i,j) = dfi / dxj, */
|
|
Packit |
67cb25 |
/* where fi = (Yi - yi)/sigma[i], */
|
|
Packit |
67cb25 |
/* Yi = A * exp(-lambda * t_i) + b */
|
|
Packit |
67cb25 |
/* and the xj are the parameters (A,lambda,b) */
|
|
Packit |
67cb25 |
double e = exp(-lambda * t[i]);
|
|
Packit |
67cb25 |
gsl_matrix_set (J, i, 0, e);
|
|
Packit |
67cb25 |
gsl_matrix_set (J, i, 1, -t[i] * A * e);
|
|
Packit |
67cb25 |
gsl_matrix_set (J, i, 2, 1.0);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
void
|
|
Packit |
67cb25 |
callback(const size_t iter, void *params,
|
|
Packit |
67cb25 |
const gsl_multifit_nlinear_workspace *w)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_vector *f = gsl_multifit_nlinear_residual(w);
|
|
Packit |
67cb25 |
gsl_vector *x = gsl_multifit_nlinear_position(w);
|
|
Packit |
67cb25 |
double rcond;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute reciprocal condition number of J(x) */
|
|
Packit |
67cb25 |
gsl_multifit_nlinear_rcond(&rcond, w);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
fprintf(stderr, "iter %2zu: A = %.4f, lambda = %.4f, b = %.4f, cond(J) = %8.4f, |f(x)| = %.4f\n",
|
|
Packit |
67cb25 |
iter,
|
|
Packit |
67cb25 |
gsl_vector_get(x, 0),
|
|
Packit |
67cb25 |
gsl_vector_get(x, 1),
|
|
Packit |
67cb25 |
gsl_vector_get(x, 2),
|
|
Packit |
67cb25 |
1.0 / rcond,
|
|
Packit |
67cb25 |
gsl_blas_dnrm2(f));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
main (void)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const gsl_multifit_nlinear_type *T = gsl_multifit_nlinear_trust;
|
|
Packit |
67cb25 |
gsl_multifit_nlinear_workspace *w;
|
|
Packit |
67cb25 |
gsl_multifit_nlinear_fdf fdf;
|
|
Packit |
67cb25 |
gsl_multifit_nlinear_parameters fdf_params =
|
|
Packit |
67cb25 |
gsl_multifit_nlinear_default_parameters();
|
|
Packit |
67cb25 |
const size_t n = N;
|
|
Packit |
67cb25 |
const size_t p = 3;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector *f;
|
|
Packit |
67cb25 |
gsl_matrix *J;
|
|
Packit |
67cb25 |
gsl_matrix *covar = gsl_matrix_alloc (p, p);
|
|
Packit |
67cb25 |
double t[N], y[N], weights[N];
|
|
Packit |
67cb25 |
struct data d = { n, t, y };
|
|
Packit |
67cb25 |
double x_init[3] = { 1.0, 1.0, 0.0 }; /* starting values */
|
|
Packit |
67cb25 |
gsl_vector_view x = gsl_vector_view_array (x_init, p);
|
|
Packit |
67cb25 |
gsl_vector_view wts = gsl_vector_view_array(weights, n);
|
|
Packit |
67cb25 |
gsl_rng * r;
|
|
Packit |
67cb25 |
double chisq, chisq0;
|
|
Packit |
67cb25 |
int status, info;
|
|
Packit |
67cb25 |
size_t i;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
const double xtol = 1e-8;
|
|
Packit |
67cb25 |
const double gtol = 1e-8;
|
|
Packit |
67cb25 |
const double ftol = 0.0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_rng_env_setup();
|
|
Packit |
67cb25 |
r = gsl_rng_alloc(gsl_rng_default);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* define the function to be minimized */
|
|
Packit |
67cb25 |
fdf.f = expb_f;
|
|
Packit |
67cb25 |
fdf.df = expb_df; /* set to NULL for finite-difference Jacobian */
|
|
Packit |
67cb25 |
fdf.fvv = NULL; /* not using geodesic acceleration */
|
|
Packit |
67cb25 |
fdf.n = n;
|
|
Packit |
67cb25 |
fdf.p = p;
|
|
Packit |
67cb25 |
fdf.params = &d;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* this is the data to be fitted */
|
|
Packit |
67cb25 |
for (i = 0; i < n; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double ti = i * TMAX / (n - 1.0);
|
|
Packit |
67cb25 |
double yi = 1.0 + 5 * exp (-0.1 * ti);
|
|
Packit |
67cb25 |
double si = 0.1 * yi;
|
|
Packit |
67cb25 |
double dy = gsl_ran_gaussian(r, si);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
t[i] = ti;
|
|
Packit |
67cb25 |
y[i] = yi + dy;
|
|
Packit |
67cb25 |
weights[i] = 1.0 / (si * si);
|
|
Packit |
67cb25 |
printf ("data: %g %g %g\n", ti, y[i], si);
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* allocate workspace with default parameters */
|
|
Packit |
67cb25 |
w = gsl_multifit_nlinear_alloc (T, &fdf_params, n, p);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* initialize solver with starting point and weights */
|
|
Packit |
67cb25 |
gsl_multifit_nlinear_winit (&x.vector, &wts.vector, &fdf, w);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute initial cost function */
|
|
Packit |
67cb25 |
f = gsl_multifit_nlinear_residual(w);
|
|
Packit |
67cb25 |
gsl_blas_ddot(f, f, &chisq0);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* solve the system with a maximum of 100 iterations */
|
|
Packit |
67cb25 |
status = gsl_multifit_nlinear_driver(100, xtol, gtol, ftol,
|
|
Packit |
67cb25 |
callback, NULL, &info, w);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute covariance of best fit parameters */
|
|
Packit |
67cb25 |
J = gsl_multifit_nlinear_jac(w);
|
|
Packit |
67cb25 |
gsl_multifit_nlinear_covar (J, 0.0, covar);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute final cost */
|
|
Packit |
67cb25 |
gsl_blas_ddot(f, f, &chisq);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#define FIT(i) gsl_vector_get(w->x, i)
|
|
Packit |
67cb25 |
#define ERR(i) sqrt(gsl_matrix_get(covar,i,i))
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
fprintf(stderr, "summary from method '%s/%s'\n",
|
|
Packit |
67cb25 |
gsl_multifit_nlinear_name(w),
|
|
Packit |
67cb25 |
gsl_multifit_nlinear_trs_name(w));
|
|
Packit |
67cb25 |
fprintf(stderr, "number of iterations: %zu\n",
|
|
Packit |
67cb25 |
gsl_multifit_nlinear_niter(w));
|
|
Packit |
67cb25 |
fprintf(stderr, "function evaluations: %zu\n", fdf.nevalf);
|
|
Packit |
67cb25 |
fprintf(stderr, "Jacobian evaluations: %zu\n", fdf.nevaldf);
|
|
Packit |
67cb25 |
fprintf(stderr, "reason for stopping: %s\n",
|
|
Packit |
67cb25 |
(info == 1) ? "small step size" : "small gradient");
|
|
Packit |
67cb25 |
fprintf(stderr, "initial |f(x)| = %f\n", sqrt(chisq0));
|
|
Packit |
67cb25 |
fprintf(stderr, "final |f(x)| = %f\n", sqrt(chisq));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double dof = n - p;
|
|
Packit |
67cb25 |
double c = GSL_MAX_DBL(1, sqrt(chisq / dof));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
fprintf(stderr, "chisq/dof = %g\n", chisq / dof);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
fprintf (stderr, "A = %.5f +/- %.5f\n", FIT(0), c*ERR(0));
|
|
Packit |
67cb25 |
fprintf (stderr, "lambda = %.5f +/- %.5f\n", FIT(1), c*ERR(1));
|
|
Packit |
67cb25 |
fprintf (stderr, "b = %.5f +/- %.5f\n", FIT(2), c*ERR(2));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
fprintf (stderr, "status = %s\n", gsl_strerror (status));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_multifit_nlinear_free (w);
|
|
Packit |
67cb25 |
gsl_matrix_free (covar);
|
|
Packit |
67cb25 |
gsl_rng_free (r);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return 0;
|
|
Packit |
67cb25 |
}
|