[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: [Bug-gsl] non-linear LS fit example in documentation: bug?
From: |
Mark M. Ito |
Subject: |
Re: [Bug-gsl] non-linear LS fit example in documentation: bug? |
Date: |
Wed, 21 Jan 2009 16:12:51 -0500 |
User-agent: |
Thunderbird 2.0.0.18 (X11/20081119) |
I've enclosed a copy of the program that is returning GSL_CONTINUE from
the solver. It is basically a fit to a parabola. It is a stand-alone
program, just needs the gsl linked in.
Brian Gough wrote:
At Tue, 20 Jan 2009 15:16:26 -0500,
Mark M. Ito wrote:
Shouldn't the "if (status) break;" be an "if (status) continue;"? It is
normal for the solver to return GSL_CONTINUE in which case you want to
continue to iterate. Break exits the do loop completely, no?
Hello,
Thanks for your email. The iterate function should return GSL_SUCCESS
or an error (if it can't make progress) so I believe the example is ok
(the test functions return GSL_CONTINUE).
If you're seeing GSL_CONTINUE return values from
gsl_multifit_fdfsolver_iterate I'd be interested to see the test
function, thanks.
/***********
function f(x) = ax^2 + bx + c
three parameters: a, b, c
7 data points: (1,~1) (2,~4) (3,~9) (4,~16) (5,~25) (6,~36) (7,~49)
7 errors on y, random
so p=3 (a, b, c)
n=7
df/da = x^2/error
df/db = x/error
df/dc = 1/error
************/
#include <stdio.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_multifit_nlin.h>
size_t n = 7;
size_t p = 3;
double xdata[7] = {1.0,2.0,3.0,4.0,5.0,6.0,7.0};
double ydata[7] = {1.12,4.13,8.94,15.75,24.86,36.27,49.28};
double edata[7] = {0.5,0.2,0.3,0.4,0.5,0.6,0.7};
int myf (const gsl_vector *x, void *data, gsl_vector *f) {
size_t i;
double a = gsl_vector_get(x, 0); /* get the parameters from the input
vector */
double b = gsl_vector_get(x, 1);
double c = gsl_vector_get(x, 2);
double func, xthis;
printf("myf params: %f %f %f\n", a, b, c);
for (i = 0; i < n; i++) {
xthis = xdata[i];
func = a*xthis*xthis + b*xthis + c;
printf("function %d %f\n", i, (ydata[i] - func)/edata[i]);
gsl_vector_set(f, i, func - ydata[i]);
}
return GSL_SUCCESS;
}
int mydf (const gsl_vector *x, void *data, gsl_matrix *J) {
size_t i;
double a = gsl_vector_get(x, 0); /* get the parameters from the input
vector */
double b = gsl_vector_get(x, 1);
double c = gsl_vector_get(x, 2);
double func, xthis, ethis;
double dfda, dfdb, dfdc;
printf("mydf params: %f %f %f\n", a, b, c);
/* Jacobian matrix J(i,j) = dfi / dxj */
for (i = 0; i < n; i++) {
xthis = xdata[i];
ethis = edata[i];
dfda = xthis*xthis/ethis;
dfdb = xthis/ethis;
dfdc = 1.0/ethis;
printf("derivs %d %f %f %f\n", i, dfda, dfdb, dfdc);
gsl_matrix_set (J, i, 0, dfda);
gsl_matrix_set (J, i, 1, dfdb);
gsl_matrix_set (J, i, 2, dfdc);
}
return GSL_SUCCESS;
}
int myfdf (const gsl_vector *x, void *data, gsl_vector *f, gsl_matrix *J) {
myf (x, data, f);
mydf (x, data, J);
return GSL_SUCCESS;
}
void print_state (size_t iter, gsl_multifit_fdfsolver * s)
{
printf ("iter: %3u x = % 15.8f % 15.8f % 15.8f "
"|f(x)| = %g\n",
iter,
gsl_vector_get (s->x, 0),
gsl_vector_get (s->x, 1),
gsl_vector_get (s->x, 2),
gsl_blas_dnrm2 (s->f));
}
int main (void) {
const gsl_multifit_fdfsolver_type *T; /* pointer to solver type */
gsl_multifit_fdfsolver *s; /* pointer to solver */
gsl_multifit_function_fdf f; /* function structure */
double x_init[3] = { 0.0, 0.0, 0.0 };
gsl_vector_view x = gsl_vector_view_array(x_init, p);
unsigned int iter; /* fit iteration */
int status; /* status from fit */
gsl_matrix *covar = gsl_matrix_alloc (p, p); /* covariance matrix */
f.f = &myf;
f.df = &mydf;
f.fdf = &myfdf;
f.n = n;
f.p = p;
f.params = NULL; /* pointer to something that gets the data, data is global
here so it is not needed */
T = gsl_multifit_fdfsolver_lmsder; /* choose the solver */
s = gsl_multifit_fdfsolver_alloc (T, n, p); /* allocate work space,
get solver pointer */
gsl_multifit_fdfsolver_set (s, &f, &x.vector); /* initialize the solver */
iter = 0;
print_state(iter, s);
do {
iter++;
printf("before solver iterate\n");
print_state(iter, s);
status = gsl_multifit_fdfsolver_iterate(s); /* go to next point */
printf("after solver iterate, status = %d, %s\n", status, gsl_strerror
(status));
print_state(iter, s);
if (status) break;
/* if (status) continue; */
status = gsl_multifit_test_delta(s->dx, s->x,
1e-4, 1e-4);
printf("after test delta, status = %d, %s\n", status, gsl_strerror
(status));
}
while (status == GSL_CONTINUE && iter < 100);
gsl_multifit_covar (s->J, 0.0, covar);
gsl_multifit_fdfsolver_free (s); /* free solver space */
gsl_matrix_free (covar); /* free covariance matrix space */
}