You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
305 lines
7.9 KiB
305 lines
7.9 KiB
=pod
|
|
|
|
=begin html
|
|
|
|
<link rel="stylesheet" href="podstyle.css" type="text/css" />
|
|
|
|
=end html
|
|
|
|
=head1 NAME
|
|
|
|
lmmin - Levenberg-Marquardt least-squares minimization
|
|
|
|
|
|
=head1 SYNOPSIS
|
|
|
|
B<#include <lmmin.h>>
|
|
|
|
B<void lmmin( const int> I<n_par>B<, double *>I<par>B<, const int> I<m_dat>B<,
|
|
constS< >void *>I<data>B<,
|
|
void *>I<evaluate>B<(
|
|
constS< >double *>I<par>B<, const int >I<m_dat>B<,
|
|
constS< >void *>I<data>B<, double *>I<fvec>B<, int *>I<userbreak>B<),
|
|
constS< >lm_control_struct *>I<control>B<,
|
|
lm_status_struct *>I<status>B< );>
|
|
|
|
B<extern const lm_control_struct lm_control_double;>
|
|
|
|
B<extern const lm_control_struct lm_control_float;>
|
|
|
|
B<extern const char *lm_infmsg[];>
|
|
|
|
B<extern const char *lm_shortmsg[];>
|
|
|
|
=head1 DESCRIPTION
|
|
|
|
B<lmmin()> determines a vector I<par> that minimizes the sum of squared elements of a vector I<fvec> that is computed by a user-supplied function I<evaluate>().
|
|
On success, I<par> represents a local minimum, not necessarily a global one; it may depend on its starting value.
|
|
|
|
For applications in curve fitting, the wrapper function B<lmcurve(3)> offers a simplified API.
|
|
|
|
The Levenberg-Marquardt minimization starts with a steepest-descent exploration of the parameter space, and achieves rapid convergence by crossing over into the Newton-Gauss method.
|
|
|
|
Function arguments:
|
|
|
|
=over
|
|
|
|
=item I<n_par>
|
|
|
|
Number of free variables.
|
|
Length of parameter vector I<par>.
|
|
|
|
=item I<par>
|
|
|
|
Parameter vector.
|
|
On input, it must contain a reasonable guess.
|
|
On output, it contains the solution found to minimize ||I<fvec>||.
|
|
|
|
=item I<m_dat>
|
|
|
|
Length of vector I<fvec>.
|
|
Must statisfy I<n_par> <= I<m_dat>.
|
|
|
|
=item I<data>
|
|
|
|
This pointer is ignored by the fit algorithm,
|
|
except for appearing as an argument in all calls to the user-supplied
|
|
routine I<evaluate>.
|
|
|
|
=item I<evaluate>
|
|
|
|
Pointer to a user-supplied function that computes I<m_dat> elements of vector I<fvec> for a given parameter vector I<par>.
|
|
If I<evaluate> return with *I<userbreak> set to a negative value, B<lmmin()> will interrupt the fitting and terminate.
|
|
|
|
=item I<control>
|
|
|
|
Parameter collection for tuning the fit procedure.
|
|
In most cases, the default &I<lm_control_double> is adequate.
|
|
If I<f> is only computed with single-precision accuracy,
|
|
I<&lm_control_float> should be used.
|
|
See also below, NOTES on initializing parameter records.
|
|
|
|
I<control> has the following members (for more details, see the source file I<lmstruct.h>):
|
|
|
|
=over
|
|
|
|
=item B<double> I<control.ftol>
|
|
|
|
Relative error desired in the sum of squares.
|
|
Recommended setting: somewhat above machine precision; less if I<fvec> is computed with reduced accuracy.
|
|
|
|
=item B<double> I<control.xtol>
|
|
|
|
Relative error between last two approximations.
|
|
Recommended setting: as I<ftol>.
|
|
|
|
=item B<double> I<control.gtol>
|
|
|
|
A measure for degeneracy.
|
|
Recommended setting: as I<ftol>.
|
|
|
|
=item B<double> I<control.epsilon>
|
|
|
|
Step used to calculate the Jacobian.
|
|
Recommended setting: as I<ftol>, but definitely less than the accuracy of I<fvec>.
|
|
|
|
=item B<double> I<control.stepbound>
|
|
|
|
Initial bound to steps in the outer loop, generally between 0.01 and 100; recommended value is 100.
|
|
|
|
=item B<int> I<control.patience>
|
|
|
|
Used to set the maximum number of function evaluations to patience*n_par.
|
|
|
|
=item B<int> I<control.scale_diag>
|
|
|
|
Logical switch (0 or 1).
|
|
If 1, then scale parameters to their initial value.
|
|
This is the recommended setting.
|
|
|
|
=item B<FILE*> I<control.msgfile>
|
|
|
|
Progress messages will be written to this file.
|
|
Typically I<stdout> or I<stderr>.
|
|
The value I<NULL> will be interpreted as I<stdout>.
|
|
|
|
=item B<int> I<control.verbosity>
|
|
|
|
If nonzero, some progress information from within the LM algorithm
|
|
is written to control.stream.
|
|
|
|
=item B<int> I<control.n_maxpri>
|
|
|
|
-1, or maximum number of parameters to print.
|
|
|
|
=item B<int> I<control.m_maxpri>
|
|
|
|
-1, or maximum number of residuals to print.
|
|
|
|
=back
|
|
|
|
=item I<status>
|
|
|
|
A record used to return information about the minimization process:
|
|
|
|
=over
|
|
|
|
=item B<double> I<status.fnorm>
|
|
|
|
Norm of the vector I<fvec>;
|
|
|
|
=item B<int> I<status.nfev>
|
|
|
|
Actual number of iterations;
|
|
|
|
=item B<int> I<status.outcome>
|
|
|
|
Status of minimization;
|
|
for the corresponding text message, print I<lm_infmsg>B<[>I<status.outcome>B<]>;
|
|
for a short code, print I<lm_shortmsg>B<[>I<status.outcome>B<]>.
|
|
|
|
=item B<int> I<status.userbreak>
|
|
|
|
Set when termination has been forced by the user-supplied routine I<evaluate>.
|
|
|
|
=back
|
|
|
|
=back
|
|
|
|
|
|
=head1 NOTES
|
|
|
|
=head2 Initializing parameter records.
|
|
|
|
The parameter record I<control> should always be initialized
|
|
from supplied default records:
|
|
|
|
lm_control_struct control = lm_control_double; /* or _float */
|
|
|
|
After this, parameters may be overwritten:
|
|
|
|
control.patience = 500; /* allow more iterations */
|
|
control.verbosity = 15; /* for verbose monitoring */
|
|
|
|
An application written this way is guaranteed to work even if new parameters
|
|
are added to I<lm_control_struct>.
|
|
|
|
Conversely, addition of parameters is not considered an API change; it may happen without increment of the major version number.
|
|
|
|
=head1 EXAMPLES
|
|
|
|
=head2 Fitting a surface
|
|
|
|
Fit a data set y(t) by a function f(t;p) where t is a two-dimensional vector:
|
|
|
|
#include "lmmin.h"
|
|
#include <stdio.h>
|
|
|
|
/* fit model: a plane p0 + p1*tx + p2*tz */
|
|
double f( double tx, double tz, const double *p )
|
|
{
|
|
return p[0] + p[1]*tx + p[2]*tz;
|
|
}
|
|
|
|
/* data structure to transmit data arays and fit model */
|
|
typedef struct {
|
|
double *tx, *tz;
|
|
double *y;
|
|
double (*f)( double tx, double tz, const double *p );
|
|
} data_struct;
|
|
|
|
/* function evaluation, determination of residues */
|
|
void evaluate_surface( const double *par, int m_dat,
|
|
const void *data, double *fvec, int *userbreak )
|
|
{
|
|
/* for readability, explicit type conversion */
|
|
data_struct *D;
|
|
D = (data_struct*)data;
|
|
|
|
int i;
|
|
for ( i = 0; i < m_dat; i++ )
|
|
fvec[i] = D->y[i] - D->f( D->tx[i], D->tz[i], par );
|
|
}
|
|
|
|
int main()
|
|
{
|
|
/* parameter vector */
|
|
int n_par = 3; /* number of parameters in model function f */
|
|
double par[3] = { -1, 0, 1 }; /* arbitrary starting value */
|
|
|
|
/* data points */
|
|
int m_dat = 4;
|
|
double tx[4] = { -1, -1, 1, 1 };
|
|
double tz[4] = { -1, 1, -1, 1 };
|
|
double y[4] = { 0, 1, 1, 2 };
|
|
|
|
data_struct data = { tx, tz, y, f };
|
|
|
|
/* auxiliary parameters */
|
|
lm_status_struct status;
|
|
lm_control_struct control = lm_control_double;
|
|
control.verbosity = 3;
|
|
|
|
/* perform the fit */
|
|
printf( "Fitting:\n" );
|
|
lmmin( n_par, par, m_dat, (const void*) &data, evaluate_surface,
|
|
&control, &status );
|
|
|
|
/* print results */
|
|
printf( "\nResults:\n" );
|
|
printf( "status after %d function evaluations:\n %s\n",
|
|
status.nfev, lm_infmsg[status.outcome] );
|
|
|
|
printf("obtained parameters:\n");
|
|
int i;
|
|
for ( i=0; i<n_par; ++i )
|
|
printf(" par[%i] = %12g\n", i, par[i]);
|
|
printf("obtained norm:\n %12g\n", status.fnorm );
|
|
|
|
printf("fitting data as follows:\n");
|
|
double ff;
|
|
for ( i=0; i<m_dat; ++i ){
|
|
ff = f(tx[i], tz[i], par);
|
|
printf( " t[%2d]=%12g,%12g y=%12g fit=%12g residue=%12g\n",
|
|
i, tx[i], tz[i], y[i], ff, y[i] - ff );
|
|
}
|
|
|
|
return 0;
|
|
}
|
|
|
|
=head2 More examples
|
|
|
|
For more examples, see the homepage and directories demo/ and test/ in the source distribution.
|
|
|
|
=head1 COPYING
|
|
|
|
Copyright (C):
|
|
1980-1999 University of Chicago
|
|
2004-2015 Joachim Wuttke, Forschungszentrum Juelich GmbH
|
|
|
|
Software: FreeBSD License
|
|
|
|
Documentation: Creative Commons Attribution Share Alike
|
|
|
|
|
|
=head1 SEE ALSO
|
|
|
|
=begin html
|
|
|
|
<a href="http://apps.jcns.fz-juelich.de/man/lmcurve.html"><b>lmcurve</b>(3)</a>
|
|
|
|
=end html
|
|
|
|
=begin man
|
|
|
|
\fBlmcurve\fR(3)
|
|
.PP
|
|
|
|
=end man
|
|
|
|
Homepage: http://apps.jcns.fz-juelich.de/lmfit
|
|
|
|
=head1 BUGS
|
|
|
|
Please send bug reports and suggestions to the author <j.wuttke@fz-juelich.de>.
|