mpfit.pro

;+
; NAME:
;   MPFIT
;
; AUTHOR:
;   Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
;   craigm@lheamail.gsfc.nasa.gov
;   UPDATED VERSIONs can be found on my WEB PAGE: 
;      http://cow.physics.wisc.edu/~craigm/idl/idl.html
;
; PURPOSE:
;   Perform Levenberg-Marquardt least-squares minimization (MINPACK-1)
;
; MAJOR TOPICS:
;   Curve and Surface Fitting
;
; CALLING SEQUENCE:
;   parms = MPFIT(MYFUNCT, start_parms, FUNCTARGS=fcnargs, NFEV=nfev,
;                 MAXITER=maxiter, ERRMSG=errmsg, NPRINT=nprint, QUIET=quiet, 
;                 FTOL=ftol, XTOL=xtol, GTOL=gtol, NITER=niter, 
;                 STATUS=status, ITERPROC=iterproc, ITERARGS=iterargs,
;                 COVAR=covar, PERROR=perror, BESTNORM=bestnorm,
;                 PARINFO=parinfo)
;
; DESCRIPTION:
;
;  MPFIT uses the Levenberg-Marquardt technique to solve the
;  least-squares problem.  In its typical use, MPFIT will be used to
;  fit a user-supplied function (the "model") to user-supplied data
;  points (the "data") by adjusting a set of parameters.  MPFIT is
;  based upon MINPACK-1 (LMDIF.F) by More' and collaborators.
;
;  For example, a researcher may think that a set of observed data
;  points is best modelled with a Gaussian curve.  A Gaussian curve is
;  parameterized by its mean, standard deviation and normalization.
;  MPFIT will, within certain constraints, find the set of parameters
;  which best fits the data.  The fit is "best" in the least-squares
;  sense; that is, the sum of the weighted squared differences between
;  the model and data is minimized.
;
;  The Levenberg-Marquardt technique is a particular strategy for
;  iteratively searching for the best fit.  This particular
;  implementation is drawn from MINPACK-1 (see NETLIB), and seems to
;  be more robust than routines provided with IDL.  This version
;  allows upper and lower bounding constraints to be placed on each
;  parameter, or the parameter can be held fixed.
;
;  The IDL user-supplied function should return an array of weighted
;  deviations between model and data.  In a typical scientific problem
;  the residuals should be weighted so that each deviate has a
;  gaussian sigma of 1.0.  If X represents values of the independent
;  variable, Y represents a measurement for each value of X, and ERR
;  represents the error in the measurements, then the deviates could
;  be calculated as follows:
;
;    DEVIATES = (Y - F(X)) / ERR
;
;  where F is the function representing the model.  You are
;  recommended to use the convenience functions MPFITFUN and
;  MPFITEXPR, which are driver functions that calculate the deviates
;  for you.  If ERR are the 1-sigma uncertainties in Y, then
;
;    TOTAL( DEVIATES^2 ) 
;
;  will be the total chi-squared value.  MPFIT will minimize the
;  chi-square value.  The values of X, Y and ERR are passed through
;  MPFIT to the user-supplied function via the FUNCTARGS keyword.
;
;  Simple constraints can be placed on parameter values by using the
;  PARINFO keyword to MPFIT.  See below for a description of this
;  keyword.
;
;  MPFIT does not perform more general optimization tasks.  See TNMIN
;  instead.  MPFIT is customized, based on MINPACK-1, to the
;  least-squares minimization problem.
;
; USER FUNCTION
;
;  The user must define a function which returns the appropriate
;  values as specified above.  The function should return the weighted
;  deviations between the model and the data.  For applications which
;  use finite-difference derivatives -- the default -- the user
;  function should be declared in the following way:
;
;    FUNCTION MYFUNCT, p, X=x, Y=y, ERR=err
;     ; Parameter values are passed in "p"
;     model = F(x, p)
;     return, (y-model)/err
;    END
;
;  See below for applications with explicit derivatives.
;
;  The keyword parameters X, Y, and ERR in the example above are
;  suggestive but not required.  Any parameters can be passed to
;  MYFUNCT by using the FUNCTARGS keyword to MPFIT.  Use MPFITFUN and
;  MPFITEXPR if you need ideas on how to do that.  The function *must*
;  accept a parameter list, P.
;  
;  In general there are no restrictions on the number of dimensions in
;  X, Y or ERR.  However the deviates *must* be returned in a
;  one-dimensional array, and must have the same type (float or
;  double) as the input arrays.
;
;  See below for error reporting mechanisms.
;
;
; CHECKING STATUS AND HANNDLING ERRORS
;
;  Upon return, MPFIT will report the status of the fitting operation
;  in the STATUS and ERRMSG keywords.  The STATUS keyword will contain
;  a numerical code which indicates the success or failure status.
;  Generally speaking, any value 1 or greater indicates success, while
;  a value of 0 or less indicates a possible failure.  The ERRMSG
;  keyword will contain a text string which should describe the error
;  condition more fully.
;
;  By default, MPFIT will trap fatal errors and report them to the
;  caller gracefully.  However, during the debugging process, it is
;  often useful to halt execution where the error occurred.  When you
;  set the NOCATCH keyword, MPFIT will not do any special error
;  trapping, and execution will stop whereever the error occurred.
;
;  MPFIT does not explicitly change the !ERROR_STATE variable
;  (although it may be changed implicitly if MPFIT calls MESSAGE).  It
;  is the caller's responsibility to call MESSAGE, /RESET to ensure
;  that the error state is initialized before calling MPFIT.
;
;  User functions may also indicate non-fatal error conditions using
;  the ERROR_CODE common block variable, as described below under the
;  MPFIT_ERROR common block definition (by setting ERROR_CODE to a
;  number between -15 and -1).  When the user function sets an error
;  condition via ERROR_CODE, MPFIT will gracefully exit immediately
;  and report this condition to the caller.  The ERROR_CODE is
;  returned in the STATUS keyword in that case.
;
;
; EXPLICIT DERIVATIVES
; 
;  In the search for the best-fit solution, MPFIT by default
;  calculates derivatives numerically via a finite difference
;  approximation.  The user-supplied function need not calculate the
;  derivatives explicitly.  However, the user function *may* calculate
;  the derivatives if desired, but only if the model function is
;  declared with an additional position parameter, DP, as described
;  below.  If the user function does not accept this additional
;  parameter, MPFIT will report an error.  As a practical matter, it
;  is often sufficient and even faster to allow MPFIT to calculate the
;  derivatives numerically, but this option is available for users who
;  wish more control over the fitting process.
;
;  There are two ways to enable explicit derivatives.  First, the user
;  can set the keyword AUTODERIVATIVE=0, which is a global switch for
;  all parameters.  In this case, MPFIT will request explicit
;  derivatives for every free parameter.  
;
;  Second, the user may request explicit derivatives for specifically
;  selected parameters using the PARINFO.MPSIDE=3 (see "CONSTRAINING
;  PARAMETER VALUES WITH THE PARINFO KEYWORD" below).  In this
;  strategy, the user picks and chooses which parameter derivatives
;  are computed explicitly versus numerically.  When PARINFO[i].MPSIDE
;  EQ 3, then the ith parameter derivative is computed explicitly.
;
;  The keyword setting AUTODERIVATIVE=0 always globally overrides the
;  individual values of PARINFO.MPSIDE.  Setting AUTODERIVATIVE=0 is
;  equivalent to resetting PARINFO.MPSIDE=3 for all parameters.
;
;  Even if the user requests explicit derivatives for some or all
;  parameters, MPFIT will not always request explicit derivatives on
;  every user function call.
;
; EXPLICIT DERIVATIVES - CALLING INTERFACE
;
;  When AUTODERIVATIVE=0, the user function is responsible for
;  calculating the derivatives of the *residuals* with respect to each
;  parameter.  The user function should be declared as follows:
;
;    ;
;    ; MYFUNCT - example user function
;    ;   P - input parameter values (N-element array)
;    ;   DP - upon input, an N-vector indicating which parameters
;    ;          to compute derivatives for; 
;    ;        upon output, the user function must return
;    ;          an ARRAY(M,N) of derivatives in this keyword
;    ;   (keywords) - any other keywords specified by FUNCTARGS
;    ; RETURNS - residual values
;    ;
;    FUNCTION MYFUNCT, p, dp, X=x, Y=y, ERR=err
;     model = F(x, p)         ;; Model function
;     resid = (y - model)/err ;; Residual calculation (for example)
;     
;     if n_params() GT 1 then begin
;       ; Create derivative and compute derivative array
;       requested = dp   ; Save original value of DP
;       dp = make_array(n_elements(x), n_elements(p), value=x[0]*0)
;
;       ; Compute derivative if requested by caller
;       for i = 0, n_elements(p)-1 do if requested(i) NE 0 then $
;         dp(*,i) = FGRAD(x, p, i) / err
;     endif
;