|
Packit |
67cb25 |
/* multifit/fdfsolver.c
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Brian Gough
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is free software; you can redistribute it and/or modify
|
|
Packit |
67cb25 |
* it under the terms of the GNU General Public License as published by
|
|
Packit |
67cb25 |
* the Free Software Foundation; either version 3 of the License, or (at
|
|
Packit |
67cb25 |
* your option) any later version.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is distributed in the hope that it will be useful, but
|
|
Packit |
67cb25 |
* WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
Packit |
67cb25 |
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
Packit |
67cb25 |
* General Public License for more details.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* You should have received a copy of the GNU General Public License
|
|
Packit |
67cb25 |
* along with this program; if not, write to the Free Software
|
|
Packit |
67cb25 |
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include <config.h>
|
|
Packit |
67cb25 |
#include <stdlib.h>
|
|
Packit |
67cb25 |
#include <string.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_math.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_errno.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_multifit_nlin.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_multifit_fdfsolver *
|
|
Packit |
67cb25 |
gsl_multifit_fdfsolver_alloc (const gsl_multifit_fdfsolver_type * T,
|
|
Packit |
67cb25 |
size_t n, size_t p)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int status;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_multifit_fdfsolver * s;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (n < p)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR_VAL ("insufficient data points, n < p", GSL_EINVAL, 0);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
s = (gsl_multifit_fdfsolver *) calloc (1, sizeof (gsl_multifit_fdfsolver));
|
|
Packit |
67cb25 |
if (s == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR_VAL ("failed to allocate space for multifit solver struct",
|
|
Packit |
67cb25 |
GSL_ENOMEM, 0);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
s->x = gsl_vector_calloc (p);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (s->x == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_multifit_fdfsolver_free (s);
|
|
Packit |
67cb25 |
GSL_ERROR_VAL ("failed to allocate space for x", GSL_ENOMEM, 0);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
s->f = gsl_vector_calloc (n);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (s->f == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_multifit_fdfsolver_free (s);
|
|
Packit |
67cb25 |
GSL_ERROR_VAL ("failed to allocate space for f", GSL_ENOMEM, 0);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
s->dx = gsl_vector_calloc (p);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (s->dx == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_multifit_fdfsolver_free (s);
|
|
Packit |
67cb25 |
GSL_ERROR_VAL ("failed to allocate space for dx", GSL_ENOMEM, 0);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
s->g = gsl_vector_alloc (p);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (s->g == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_multifit_fdfsolver_free (s);
|
|
Packit |
67cb25 |
GSL_ERROR_VAL ("failed to allocate space for g", GSL_ENOMEM, 0);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
s->sqrt_wts = gsl_vector_calloc (n);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (s->sqrt_wts == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_multifit_fdfsolver_free (s);
|
|
Packit |
67cb25 |
GSL_ERROR_VAL ("failed to allocate space for sqrt_wts", GSL_ENOMEM, 0);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
s->state = calloc (1, T->size);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (s->state == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_multifit_fdfsolver_free (s);
|
|
Packit |
67cb25 |
GSL_ERROR_VAL ("failed to allocate space for multifit solver state",
|
|
Packit |
67cb25 |
GSL_ENOMEM, 0);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
s->type = T ;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
status = (s->type->alloc)(s->state, n, p);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (status != GSL_SUCCESS)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_multifit_fdfsolver_free (s);
|
|
Packit |
67cb25 |
GSL_ERROR_VAL ("failed to set solver", status, 0);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
s->fdf = NULL;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
s->niter = 0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return s;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_multifit_fdfsolver_set (gsl_multifit_fdfsolver * s,
|
|
Packit |
67cb25 |
gsl_multifit_function_fdf * f,
|
|
Packit |
67cb25 |
const gsl_vector * x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return gsl_multifit_fdfsolver_wset(s, f, x, NULL);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_multifit_fdfsolver_wset (gsl_multifit_fdfsolver * s,
|
|
Packit |
67cb25 |
gsl_multifit_function_fdf * f,
|
|
Packit |
67cb25 |
const gsl_vector * x,
|
|
Packit |
67cb25 |
const gsl_vector * wts)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t n = s->f->size;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (n != f->n)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("function size does not match solver", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (s->x->size != x->size)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("vector length does not match solver", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (wts != NULL && n != wts->size)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("weight vector length does not match solver", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t i;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
s->fdf = f;
|
|
Packit |
67cb25 |
gsl_vector_memcpy(s->x, x);
|
|
Packit |
67cb25 |
s->niter = 0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (wts)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
for (i = 0; i < n; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double wi = gsl_vector_get(wts, i);
|
|
Packit |
67cb25 |
gsl_vector_set(s->sqrt_wts, i, sqrt(wi));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
gsl_vector_set_all(s->sqrt_wts, 1.0);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return (s->type->set) (s->state, s->sqrt_wts, s->fdf, s->x, s->f, s->dx);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_multifit_fdfsolver_iterate (gsl_multifit_fdfsolver * s)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int status =
|
|
Packit |
67cb25 |
(s->type->iterate) (s->state, s->sqrt_wts, s->fdf, s->x, s->f, s->dx);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
s->niter++;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return status;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
gsl_multifit_fdfsolver_driver()
|
|
Packit |
67cb25 |
Iterate the nonlinear least squares solver until completion
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Inputs: s - fdfsolver
|
|
Packit |
67cb25 |
maxiter - maximum iterations to allow
|
|
Packit |
67cb25 |
xtol - tolerance in step x
|
|
Packit |
67cb25 |
gtol - tolerance in gradient
|
|
Packit |
67cb25 |
ftol - tolerance in ||f||
|
|
Packit |
67cb25 |
info - (output) info flag on why iteration terminated
|
|
Packit |
67cb25 |
1 = stopped due to small step size ||dx|
|
|
Packit |
67cb25 |
2 = stopped due to small gradient
|
|
Packit |
67cb25 |
3 = stopped due to small change in f
|
|
Packit |
67cb25 |
GSL_ETOLX = ||dx|| has converged to within machine
|
|
Packit |
67cb25 |
precision (and xtol is too small)
|
|
Packit |
67cb25 |
GSL_ETOLG = ||g||_inf is smaller than machine
|
|
Packit |
67cb25 |
precision (gtol is too small)
|
|
Packit |
67cb25 |
GSL_ETOLF = change in ||f|| is smaller than machine
|
|
Packit |
67cb25 |
precision (ftol is too small)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Return: GSL_SUCCESS if converged, GSL_MAXITER if maxiter exceeded without
|
|
Packit |
67cb25 |
converging
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_multifit_fdfsolver_driver (gsl_multifit_fdfsolver * s,
|
|
Packit |
67cb25 |
const size_t maxiter,
|
|
Packit |
67cb25 |
const double xtol,
|
|
Packit |
67cb25 |
const double gtol,
|
|
Packit |
67cb25 |
const double ftol,
|
|
Packit |
67cb25 |
int *info)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int status;
|
|
Packit |
67cb25 |
size_t iter = 0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
do
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
status = gsl_multifit_fdfsolver_iterate (s);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* if status is GSL_ENOPROG or GSL_SUCCESS, continue iterating,
|
|
Packit |
67cb25 |
* otherwise the method has converged with a GSL_ETOLx flag
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
if (status != GSL_SUCCESS && status != GSL_ENOPROG)
|
|
Packit |
67cb25 |
break;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* test for convergence */
|
|
Packit |
67cb25 |
status = gsl_multifit_fdfsolver_test(s, xtol, gtol, ftol, info);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
while (status == GSL_CONTINUE && ++iter < maxiter);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* the following error codes mean that the solution has converged
|
|
Packit |
67cb25 |
* to within machine precision, so record the error code in info
|
|
Packit |
67cb25 |
* and return success
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
if (status == GSL_ETOLF || status == GSL_ETOLX || status == GSL_ETOLG)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
*info = status;
|
|
Packit |
67cb25 |
status = GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* check if max iterations reached */
|
|
Packit |
67cb25 |
if (iter >= maxiter && status != GSL_SUCCESS)
|
|
Packit |
67cb25 |
status = GSL_EMAXITER;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return status;
|
|
Packit |
67cb25 |
} /* gsl_multifit_fdfsolver_driver() */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_multifit_fdfsolver_jac (gsl_multifit_fdfsolver * s, gsl_matrix * J)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t n = s->f->size;
|
|
Packit |
67cb25 |
const size_t p = s->x->size;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (n != J->size1 || p != J->size2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("Jacobian dimensions do not match solver", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return (s->type->jac) (s->state, J);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
} /* gsl_multifit_fdfsolver_jac() */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
void
|
|
Packit |
67cb25 |
gsl_multifit_fdfsolver_free (gsl_multifit_fdfsolver * s)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
RETURN_IF_NULL (s);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (s->state)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
(s->type->free) (s->state);
|
|
Packit |
67cb25 |
free (s->state);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (s->dx)
|
|
Packit |
67cb25 |
gsl_vector_free (s->dx);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (s->x)
|
|
Packit |
67cb25 |
gsl_vector_free (s->x);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (s->f)
|
|
Packit |
67cb25 |
gsl_vector_free (s->f);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (s->sqrt_wts)
|
|
Packit |
67cb25 |
gsl_vector_free (s->sqrt_wts);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (s->g)
|
|
Packit |
67cb25 |
gsl_vector_free (s->g);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
free (s);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
const char *
|
|
Packit |
67cb25 |
gsl_multifit_fdfsolver_name (const gsl_multifit_fdfsolver * s)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return s->type->name;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector *
|
|
Packit |
67cb25 |
gsl_multifit_fdfsolver_position (const gsl_multifit_fdfsolver * s)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return s->x;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector *
|
|
Packit |
67cb25 |
gsl_multifit_fdfsolver_residual (const gsl_multifit_fdfsolver * s)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return s->f;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
size_t
|
|
Packit |
67cb25 |
gsl_multifit_fdfsolver_niter (const gsl_multifit_fdfsolver * s)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return s->niter;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
gsl_multifit_eval_wf()
|
|
Packit |
67cb25 |
Compute residual vector y with user callback function, and apply
|
|
Packit |
67cb25 |
weighting transform if given:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
y~ = sqrt(W) y
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Inputs: fdf - callback function
|
|
Packit |
67cb25 |
x - model parameters
|
|
Packit |
67cb25 |
swts - weight matrix sqrt(W) = sqrt(diag(w1,w2,...,wn))
|
|
Packit |
67cb25 |
set to NULL for unweighted fit
|
|
Packit |
67cb25 |
y - (output) (weighted) residual vector
|
|
Packit |
67cb25 |
y_i = sqrt(w_i) f_i where f_i is unweighted residual
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_multifit_eval_wf(gsl_multifit_function_fdf *fdf, const gsl_vector *x,
|
|
Packit |
67cb25 |
const gsl_vector *swts, gsl_vector *y)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int s = ((*((fdf)->f)) (x, fdf->params, y));
|
|
Packit |
67cb25 |
++(fdf->nevalf);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* y <- sqrt(W) y */
|
|
Packit |
67cb25 |
if (swts)
|
|
Packit |
67cb25 |
gsl_vector_mul(y, swts);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return s;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
gsl_multifit_eval_wdf()
|
|
Packit |
67cb25 |
Compute Jacobian matrix J with user callback function, and apply
|
|
Packit |
67cb25 |
weighting transform if given:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
J~ = sqrt(W) J
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Inputs: fdf - callback function
|
|
Packit |
67cb25 |
x - model parameters
|
|
Packit |
67cb25 |
swts - weight matrix W = diag(w1,w2,...,wn)
|
|
Packit |
67cb25 |
set to NULL for unweighted fit
|
|
Packit |
67cb25 |
dy - (output) (weighted) Jacobian matrix
|
|
Packit |
67cb25 |
dy = sqrt(W) dy where dy is unweighted Jacobian
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_multifit_eval_wdf(gsl_multifit_function_fdf *fdf, const gsl_vector *x,
|
|
Packit |
67cb25 |
const gsl_vector *swts, gsl_matrix *dy)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int status = ((*((fdf)->df)) (x, fdf->params, dy));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
++(fdf->nevaldf);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* J <- sqrt(W) J */
|
|
Packit |
67cb25 |
if (swts)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t n = swts->size;
|
|
Packit |
67cb25 |
size_t i;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < n; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double swi = gsl_vector_get(swts, i);
|
|
Packit |
67cb25 |
gsl_vector_view v = gsl_matrix_row(dy, i);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_scale(&v.vector, swi);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return status;
|
|
Packit |
67cb25 |
}
|