mpfit.pro

非线性最小二乘-LM法 是用IDL语言编写的

;
;             Default value:  all parameters are free and unconstrained.
;
;   PERROR - The formal 1-sigma errors in each parameter, computed
;            from the covariance matrix.  If a parameter is held
;            fixed, or if it touches a boundary, then the error is
;            reported as zero.
;
;            If the fit is unweighted (i.e. no errors were given, or
;            the weights were uniformly set to unity), then PERROR
;            will probably not represent the true parameter
;            uncertainties.  
;
;            *If* you can assume that the true reduced chi-squared
;            value is unity -- meaning that the fit is implicitly
;            assumed to be of good quality -- then the estimated
;            parameter uncertainties can be computed by scaling PERROR
;            by the measured chi-squared value.
;
;              DOF     = N_ELEMENTS(X) - N_ELEMENTS(PARMS) ; deg of freedom
;              PCERROR = PERROR * SQRT(BESTNORM / DOF)   ; scaled uncertainties
;
;   QUIET - set this keyword when no textual output should be printed
;           by MPFIT
;
;   RESDAMP - a scalar number, indicating the cut-off value of
;             residuals where "damping" will occur.  Residuals with
;             magnitudes greater than this number will be replaced by
;             their logarithm.  This partially mitigates the so-called
;             large residual problem inherent in least-squares solvers
;             (as for the test problem CURVI, http://www.maxthis.com/-
;             curviex.htm).  A value of 0 indicates no damping.
;             Default: 0
;
;             Note: RESDAMP doesn't work with AUTODERIV=0
;
;   STATUS - an integer status code is returned.  All values greater
;            than zero can represent success (however STATUS EQ 5 may
;            indicate failure to converge).  It can have one of the
;            following values:
;
;        -18  a fatal execution error has occurred.  More information
;             may be available in the ERRMSG string.
;
;        -16  a parameter or function value has become infinite or an
;             undefined number.  This is usually a consequence of
;             numerical overflow in the user's model function, which
;             must be avoided.
;
;        -15 to -1 
;             these are error codes that either MYFUNCT or ITERPROC
;             may return to terminate the fitting process (see
;             description of MPFIT_ERROR common below).  If either
;             MYFUNCT or ITERPROC set ERROR_CODE to a negative number,
;             then that number is returned in STATUS.  Values from -15
;             to -1 are reserved for the user functions and will not
;             clash with MPFIT.
;
;	   0  improper input parameters.
;         
;	   1  both actual and predicted relative reductions
;	      in the sum of squares are at most FTOL.
;         
;	   2  relative error between two consecutive iterates
;	      is at most XTOL
;         
;	   3  conditions for STATUS = 1 and STATUS = 2 both hold.
;         
;	   4  the cosine of the angle between fvec and any
;	      column of the jacobian is at most GTOL in
;	      absolute value.
;         
;	   5  the maximum number of iterations has been reached
;         
;	   6  FTOL is too small. no further reduction in
;	      the sum of squares is possible.
;         
;	   7  XTOL is too small. no further improvement in
;	      the approximate solution x is possible.
;         
;	   8  GTOL is too small. fvec is orthogonal to the
;	      columns of the jacobian to machine precision.
;
;          9  A successful single iteration has been completed, and
;             the user must supply another "EXTERNAL" evaluation of
;             the function and its derivatives.  This status indicator
;             is neither an error nor a convergence indicator.
;
;   XTOL - a nonnegative input variable. Termination occurs when the
;          relative error between two consecutive iterates is at most
;          XTOL (and STATUS is accordingly set to 2 or 3).  Therefore,
;          XTOL measures the relative error desired in the approximate
;          solution.  Default: 1D-10
;
;
; EXAMPLE:
;
;   p0 = [5.7D, 2.2, 500., 1.5, 2000.]
;   fa = {X:x, Y:y, ERR:err}
;   p = mpfit('MYFUNCT', p0, functargs=fa)
;
;   Minimizes sum of squares of MYFUNCT.  MYFUNCT is called with the X,
;   Y, and ERR keyword parameters that are given by FUNCTARGS.  The
;   resulting parameter values are returned in p.
;
;
; COMMON BLOCKS:
;
;   COMMON MPFIT_ERROR, ERROR_CODE
;
;     User routines may stop the fitting process at any time by
;     setting an error condition.  This condition may be set in either
;     the user's model computation routine (MYFUNCT), or in the
;     iteration procedure (ITERPROC).
;
;     To stop the fitting, the above common block must be declared,
;     and ERROR_CODE must be set to a negative number.  After the user
;     procedure or function returns, MPFIT checks the value of this
;     common block variable and exits immediately if the error
;     condition has been set.  This value is also returned in the
;     STATUS keyword: values of -1 through -15 are reserved error
;     codes for the user routines.  By default the value of ERROR_CODE
;     is zero, indicating a successful function/procedure call.
;
;   COMMON MPFIT_PROFILE
;   COMMON MPFIT_MACHAR
;   COMMON MPFIT_CONFIG
;
;     These are undocumented common blocks are used internally by
;     MPFIT and may change in future implementations.
;
; THEORY OF OPERATION:
;
;   There are many specific strategies for function minimization.  One
;   very popular technique is to use function gradient information to
;   realize the local structure of the function.  Near a local minimum
;   the function value can be taylor expanded about x0 as follows:
;
;      f(x) = f(x0) + f'(x0) . (x-x0) + (1/2) (x-x0) . f''(x0) . (x-x0)
;             -----   ---------------   -------------------------------  (1)
;     Order    0th          1st                      2nd
;
;   Here f'(x) is the gradient vector of f at x, and f''(x) is the
;   Hessian matrix of second derivatives of f at x.  The vector x is
;   the set of function parameters, not the measured data vector.  One
;   can find the minimum of f, f(xm) using Newton's method, and
;   arrives at the following linear equation:
;
;      f''(x0) . (xm-x0) = - f'(x0)                            (2)
;
;   If an inverse can be found for f''(x0) then one can solve for
;   (xm-x0), the step vector from the current position x0 to the new
;   projected minimum.  Here the problem has been linearized (ie, the
;   gradient information is known to first order).  f''(x0) is
;   symmetric n x n matrix, and should be positive definite.
;
;   The Levenberg - Marquardt technique is a variation on this theme.
;   It adds an additional diagonal term to the equation which may aid the
;   convergence properties:
;
;      (f''(x0) + nu I) . (xm-x0) = -f'(x0)                  (2a)
;
;   where I is the identity matrix.  When nu is large, the overall
;   matrix is diagonally dominant, and the iterations follow steepest
;   descent.  When nu is small, the iterations are quadratically
;   convergent.
;
;   In principle, if f''(x0) and f'(x0) are known then xm-x0 can be
;   determined.  However the Hessian matrix is often difficult or
;   impossible to compute.  The gradient f'(x0) may be easier to
;   compute, if even by finite difference techniques.  So-called
;   quasi-Newton techniques attempt to successively estimate f''(x0)
;   by building up gradient information as the iterations proceed.
;
;   In the least squares problem there are further simplifications
;   which assist in solving eqn (2).  The function to be minimized is
;   a sum of squares:
;
;       f = Sum(hi^2)                                         (3)
;
;   where hi is the ith residual out of m residuals as described
;   above.  This can be substituted back into eqn (2) after computing
;   the derivatives:
;
;       f'  = 2 Sum(hi  hi')     
;       f'' = 2 Sum(hi' hj') + 2 Sum(hi hi'')                (4)
;
;   If one assumes that the parameters are already close enough to a
;   minimum, then one typically finds that the second term in f'' is
;   negligible [or, in any case, is too difficult to compute].  Thus,
;   equation (2) can be solved, at least approximately, using only
;   gradient information.
;
;   In matrix notation, the combination of eqns (2) and (4) becomes:
;
;        hT' . h' . dx = - hT' . h                          (5)
;
;   Where h is the residual vector (length m), hT is its transpose, h'
;   is the Jacobian matrix (dimensions n x m), and dx is (xm-x0).  The
;   user function supplies the residual vector h, and in some cases h'
;   when it is not found by finite differences (see MPFIT_FDJAC2,
;   which finds h and hT').  Even if dx is not the best absolute step
;   to take, it does provide a good estimate of the best *direction*,
;   so often a line minimization will occur along the dx vector
;   direction.
;
;   The method of solution employed by MINPACK is to form the Q . R
;   factorization of h', where Q is an orthogonal matrix such that QT .
;   Q = I, and R is upper right triangular.  Using h' = Q . R and the
;   ortogonality of Q, eqn (5) becomes
;
;        (RT . QT) . (Q . R) . dx = - (RT . QT) . h
;                     RT . R . dx = - RT . QT . h         (6)
;                          R . dx = - QT . h
;
;   where the last statement follows because R is upper triangular.
;   Here, R, QT and h are known so this is a matter of solving for dx.
;   The routine MPFIT_QRFAC provides the QR factorization of h, with
;   pivoting, and MPFIT_QRSOLV provides the solution for dx.
;   
; REFERENCES:
;
;   MINPACK-1, Jorge More', available from netlib (www.netlib.org).
;   "Optimization Software Guide," Jorge More' and Stephen Wright, 
;     SIAM, *Frontiers in Applied Mathematics*, Number 14.
;   More', Jorge J., "The Levenberg-Marquardt Algorithm:
;     Implementation and Theory," in *Numerical Analysis*, ed. Watson,
;     G. A., Lecture Notes in Mathematics 630, Springer-Verlag, 1977.
;
; MODIFICATION HISTORY:
;   Translated from MINPACK-1 in FORTRAN, Apr-Jul 1998, CM
;   Fixed bug in parameter limits (x vs xnew), 04 Aug 1998, CM
;   Added PERROR keyword, 04 Aug 1998, CM
;   Added COVAR keyword, 20 Aug 1998, CM
;   Added NITER output keyword, 05 Oct 1998
;      D.L Windt, Bell Labs, [e-mail];
;   Made each PARINFO component optional, 05 Oct 1998 CM
;   Analytical derivatives allowed via AUTODERIVATIVE keyword, 09 Nov 1998
;   Parameter values can be tied to others, 09 Nov 1998
;   Fixed small bugs (Wayne Landsman), 24 Nov 1998
;   Added better exception error reporting, 24 Nov 1998 CM
;   Cosmetic documentation changes, 02 Jan 1999 CM
;   Changed definition of ITERPROC to be consistent with TNMIN, 19 Jan 1999 CM
;   Fixed bug when AUTDERIVATIVE=0.  Incorrect sign, 02 Feb 1999 CM
;   Added keyboard stop to MPFIT_DEFITER, 28 Feb 1999 CM
;   Cosmetic documentation changes, 14 May 1999 CM
;   IDL optimizations for speed & FASTNORM keyword, 15 May 1999 CM
;   Tried a faster version of mpfit_enorm, 30 May 1999 CM
;   Changed web address to cow.physics.wisc.edu, 14 Jun 1999 CM
;   Found malformation of FDJAC in MPFIT for 1 parm, 03 Aug 1999 CM
;   Factored out user-function call into MPFIT_CALL.  It is possible,
;     but currently disabled, to call procedures.  The calling format
;     is similar to CURVEFIT, 25 Sep 1999, CM
;   Slightly changed mpfit_tie to be less intrusive, 25 Sep 1999, CM
;   Fixed some bugs associated with tied parameters in mpfit_fdjac, 25
;     Sep 1999, CM
;   Reordered documentation; now alphabetical, 02 Oct 1999, CM
;   Added QUERY keyword for more robust error detection in drivers, 29
;     Oct 1999, CM
;   Documented PERROR for unweighted fits, 03 Nov 1999, CM
;   Split out MPFIT_RESETPROF to aid in profiling, 03 Nov 1999, CM
;   Some profiling and speed optimization, 03 Nov 1999, CM
;     Worst offenders, in order: fdjac2, qrfac, qrsolv, enorm.
;     fdjac2 depends on user function, qrfac and enorm seem to be
;     fully optimized.  qrsolv probably could be tweaked a little, but
;     is still <10% of total compute time.
;   Made sure that !err was set to 0 in MPFIT_DEFITER, 10 Jan 2000, CM
;   Fixed small inconsistency in setting of QANYLIM, 28 Jan 2000, CM
;   Added PARINFO field RELSTEP, 28 Jan 2000, CM
;   Converted to MPFIT_ERROR common block for indicating error
;     conditions, 28 Jan 2000, CM
;   Corrected scope of MPFIT_ERROR common block, CM, 07 Mar 2000
;   Minor speed improvement in MPFIT_ENORM, CM 26 Mar 2000
;   Corrected case where ITERPROC changed parameter values and
;     parameter values were TIED, CM 26 Mar 2000
;   Changed MPFIT_CALL to modify NFEV automatically, and to support
;     user procedures more, CM 26 Mar 2000
;   Copying permission terms have been liberalized, 26 Mar 2000, CM
;   Catch zero value of zero a(j,lj) in MPFIT_QRFAC, 20 Jul 2000, CM
;      (thanks to David Schlegel <[e-mail]>)
;   MPFIT_SETMACHAR is called only once at init; only one common block
;     is created (MPFIT_MACHAR); it is now a structure; removed almost
;     all CHECK_MATH calls for compatibility with IDL5 and !EXCEPT;
;     profiling data is now in a structure too; noted some
;     mathematical discrepancies in Linux IDL5.0, 17 Nov 2000, CM
;   Some significant changes.  New PARINFO fields: MPSIDE, MPMINSTEP,
;     MPMAXSTEP.  Improved documentation.  Now PTIED constraints are
;     maintained in the MPCONFIG common block.  A new procedure to
;     parse PARINFO fields.  FDJAC2 now computes a larger variety of
;     one-sided and two-sided finite difference derivatives.  NFEV is
;     stored in the MPCONFIG common now.  17 Dec 2000, CM
;   Added check that PARINFO and XALL have same size, 29 Dec 2000 CM
;   Don't call function in TERMINATE when there is an error, 05 Jan
;     2000
;   Check for float vs. double discrepancies; corrected implementation
;     of MIN/MAXSTEP, which I still am not sure of, but now at least
;     the correct behavior occurs *without* it, CM 08 Jan 2001
;   Added SCALE_FCN keyword, to allow for scaling, as for the CASH
;     statistic; added documentation about the theory of operation,
;     and under the QR factorization; slowly I'm beginning to
;     understand the bowels of this algorithm, CM 10 Jan 2001
;   Remove MPMINSTEP field of PARINFO, for now at least, CM 11 Jan
;     2001
;   Added RESDAMP keyword, CM, 14 Jan 2001
;   Tried to improve the DAMP handling a little, CM, 13 Mar 2001
;   Corrected .PARNAME behavior in _DEFITER, CM, 19 Mar 2001