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
;    
;     return, resid
;    END
;
;  where FGRAD(x, p, i) is a model function which computes the
;  derivative of the model F(x,p) with respect to parameter P(i) at X.
;
;  A quirk in the implementation leaves a stray negative sign in the
;  definition of DP.  The derivative of the *residual* should be
;  "-FGRAD(x,p,i) / err" because of how the residual is defined
;  ("resid = (data - model) / err").  **HOWEVER** because of the
;  implementation quirk, MPFIT expects FGRAD(x,p,i)/err instead,
;  i.e. the opposite sign of the gradient of RESID.
;
;  Derivatives should be returned in the DP array. DP should be an
;  ARRAY(m,n) array, where m is the number of data points and n is the
;  number of parameters.  -DP[i,j] is the derivative of the ith
;  residual with respect to the jth parameter (note the minus sign
;  due to the quirk described above).
;
;  As noted above, MPFIT may not always request derivatives from the
;  user function.  In those cases, the parameter DP is not passed.
;  Therefore functions can use N_PARAMS() to indicate whether they
;  must compute the derivatives or not.
;  
;  The derivatives with respect to fixed parameters are ignored; zero
;  is an appropriate value to insert for those derivatives.  Upon
;  input to the user function, DP is set to a vector with the same
;  length as P, with a value of 1 for a parameter which is free, and a
;  value of zero for a parameter which is fixed (and hence no
;  derivative needs to be calculated).  This input vector may be
;  overwritten as needed.  In the example above, the original DP
;  vector is saved to a variable called REQUESTED, and used as a mask
;  to calculate only those derivatives that are required.
;
;  If the data is higher than one dimensional, then the *last*
;  dimension should be the parameter dimension.  Example: fitting a
;  50x50 image, "dp" should be 50x50xNPAR.
;
; EXPLICIT DERIVATIVES - TESTING and DEBUGGING
;
;  For reasonably complicated user functions, the calculation of
;  explicit derivatives of the correct sign and magnitude can be
;  difficult to get right.  A simple sign error can cause MPFIT to be
;  confused.  MPFIT has a derivative debugging mode which will compute
;  the derivatives *both* numerically and explicitly, and compare the
;  results.
;
;  It is expected that during production usage, derivative debugging
;  should be disabled for all parameters.
;
;  In order to enable derivative debugging mode, set the following
;  PARINFO members for the ith parameter.
;      PARINFO[i].MPSIDE = 3          ; Enable explicit derivatives
;      PARINFO[i].MPDERIV_DEBUG = 1   ; Enable derivative debugging mode
;      PARINFO[i].MPDERIV_RELTOL = ?? ; Relative tolerance for comparison
;      PARINFO[i].MPDERIV_ABSTOL = ?? ; Absolute tolerance for comparison
;  Note that these settings are maintained on a parameter-by-parameter
;  basis using PARINFO, so the user can choose which parameters
;  derivatives will be tested.
;
;  When .MPDERIV_DEBUG is set, then MPFIT first computes the
;  derivative explicitly by requesting them from the user function.
;  Then, it computes the derivatives numerically via finite
;  differencing, and compares the two values.  If the difference
;  exceeds a tolerance threshold, then the values are printed out to 
;  alert the user.  The tolerance level threshold contains both a
;  relative and an absolute component, and is expressed as,
;
;     ABS(DERIV_U - DERIV_N) GE (ABSTOL + RELTOL*ABS(DERIV_U))
;
;  where DERIV_U and DERIV_N are the derivatives computed explicitly
;  and numerically, respectively.  Appropriate values
;  for most users will be: 
;
;      PARINFO[i].MPDERIV_RELTOL = 1d-3 ;; Suggested relative tolerance 
;      PARINFO[i].MPDERIV_ABSTOL = 1d-7 ;; Suggested absolute tolerance
;
;  although these thresholds may have to be adjusted for a particular
;  problem.  When the threshold is exceeded, users can expect to see a
;  tabular report like this one:
;
;    FJAC DEBUG BEGIN
;    #        IPNT       FUNC    DERIV_U    DERIV_N   DIFF_ABS   DIFF_REL
;    FJAC PARM 2
;               80    -0.7308    0.04233    0.04233 -5.543E-07 -1.309E-05
;               99      1.370    0.01417    0.01417 -5.518E-07 -3.895E-05
;              118    0.07187   -0.01400   -0.01400 -5.566E-07  3.977E-05
;              137      1.844   -0.04216   -0.04216 -5.589E-07  1.326E-05
;    FJAC DEBUG END
;
;  The report will be bracketed by FJAC DEBUG BEGIN/END statements.
;  Each parameter will be delimited by the statement FJAC PARM n,
;  where n is the parameter number.  The columns are,
;
;      IPNT - data point number  (0 ... M-1)
;      FUNC - function value at that point
;      DERIV_U - explicit derivative value at that point
;      DERIV_N - numerical derivative estimate at that point
;      DIFF_ABS - absolute difference = (DERIV_U - DERIV_N)
;      DIFF_REL - relative difference = (DIFF_ABS)/(DERIV_U)
;
;  When prints appear in this report, it is most important to check
;  that the derivatives computed in two different ways have the same
;  numerical sign and the same order of magnitude, since these are the
;  most common programming mistakes.
;
;  A line of this form may also appear 
;
;   # FJAC_MASK = 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 
;
;  This line indicates for which parameters explicit derivatives are
;  expected.  A list of all-1s indicates all explicit derivatives for
;  all parameters are requested from the user function.
;    
;  
; CONSTRAINING PARAMETER VALUES WITH THE PARINFO KEYWORD
;
;  The behavior of MPFIT can be modified with respect to each
;  parameter to be fitted.  A parameter value can be fixed; simple
;  boundary constraints can be imposed; limitations on the parameter
;  changes can be imposed; properties of the automatic derivative can
;  be modified; and parameters can be tied to one another.
;
;  These properties are governed by the PARINFO structure, which is
;  passed as a keyword parameter to MPFIT.
;
;  PARINFO should be an array of structures, one for each parameter.
;  Each parameter is associated with one element of the array, in
;  numerical order.  The structure can have the following entries
;  (none are required):
;  
;     .VALUE - the starting parameter value (but see the START_PARAMS
;              parameter for more information).
;  
;     .FIXED - a boolean value, whether the parameter is to be held
;              fixed or not.  Fixed parameters are not varied by
;              MPFIT, but are passed on to MYFUNCT for evaluation.
;  
;     .LIMITED - a two-element boolean array.  If the first/second
;                element is set, then the parameter is bounded on the
;                lower/upper side.  A parameter can be bounded on both
;                sides.  Both LIMITED and LIMITS must be given
;                together.
;  
;     .LIMITS - a two-element float or double array.  Gives the
;               parameter limits on the lower and upper sides,
;               respectively.  Zero, one or two of these values can be
;               set, depending on the values of LIMITED.  Both LIMITED
;               and LIMITS must be given together.
;  
;     .PARNAME - a string, giving the name of the parameter.  The
;                fitting code of MPFIT does not use this tag in any
;                way.  However, the default ITERPROC will print the
;                parameter name if available.
;  
;     .STEP - the step size to be used in calculating the numerical
;             derivatives.  If set to zero, then the step size is
;             computed automatically.  Ignored when AUTODERIVATIVE=0.
;             This value is superceded by the RELSTEP value.
;
;     .RELSTEP - the *relative* step size to be used in calculating
;                the numerical derivatives.  This number is the
;                fractional size of the step, compared to the
;                parameter value.  This value supercedes the STEP
;                setting.  If the parameter is zero, then a default
;                step size is chosen.
;
;     .MPSIDE - selector for type of derivative calculation. This
;               field can take one of five possible values:
;
;                  0 - one-sided derivative computed automatically
;                  1 - one-sided derivative (f(x+h) - f(x)  )/h
;                 -1 - one-sided derivative (f(x)   - f(x-h))/h
;                  2 - two-sided derivative (f(x+h) - f(x-h))/(2*h)
;                  3 - explicit derivative used for this parameter
;
;              In the first four cases, the derivative is approximated
;              numerically by finite difference, with step size
;              H=STEP, where the STEP parameter is defined above.  The
;              last case, MPSIDE=3, indicates to allow the user
;              function to compute the derivative explicitly (see
;              section on "EXPLICIT DERIVATIVES").  AUTODERIVATIVE=0
;              overrides this setting for all parameters, and is
;              equivalent to MPSIDE=3 for all parameters.  For
;              MPSIDE=0, the "automatic" one-sided derivative method
;              will chose a direction for the finite difference which
;              does not violate any constraints.  The other methods
;              (MPSIDE=-1 or MPSIDE=1) do not perform this check.  The
;              two-sided method is in principle more precise, but
;              requires twice as many function evaluations.  Default:
;              0.
;
;     .MPDERIV_DEBUG - set this value to 1 to enable debugging of
;              user-supplied explicit derivatives (see "TESTING and
;              DEBUGGING" section above).  In addition, the
;              user must enable calculation of explicit derivatives by
;              either setting AUTODERIVATIVE=0, or MPSIDE=3 for the
;              desired parameters.  When this option is enabled, a
;              report may be printed to the console, depending on the
;              MPDERIV_ABSTOL and MPDERIV_RELTOL settings.
;              Default: 0 (no debugging)
;
;     
;     .MPDERIV_ABSTOL, .MPDERIV_RELTOL - tolerance settings for
;              print-out of debugging information, for each parameter
;              where debugging is enabled.  See "TESTING and
;              DEBUGGING" section above for the meanings of these two
;              fields.
;
;
;     .MPMAXSTEP - the maximum change to be made in the parameter
;                  value.  During the fitting process, the parameter
;                  will never be changed by more than this value in
;                  one iteration.
;
;                  A value of 0 indicates no maximum.  Default: 0.
;  
;     .TIED - a string expression which "ties" the parameter to other
;             free or fixed parameters as an equality constraint.  Any
;             expression involving constants and the parameter array P
;             are permitted.
;             Example: if parameter 2 is always to be twice parameter
;             1 then use the following: parinfo[2].tied = '2 * P[1]'.
;             Since they are totally constrained, tied parameters are
;             considered to be fixed; no errors are computed for them,
;             and any LIMITS are not obeyed.
;             [ NOTE: the PARNAME can't be used in a TIED expression. ]
;
;     .MPPRINT - if set to 1, then the default ITERPROC will print the
;                parameter value.  If set to 0, the parameter value
;                will not be printed.  This tag can be used to
;                selectively print only a few parameter values out of
;                many.  Default: 1 (all parameters printed)
;
;     .MPFORMAT - IDL format string to print the parameter within
;                 ITERPROC.  Default: '(G20.6)'  (An empty string will
;                 also use the default.)
;
;  Future modifications to the PARINFO structure, if any, will involve
;  adding structure tags beginning with the two letters "MP".
;  Therefore programmers are urged to avoid using tags starting with
;  "MP", but otherwise they are free to include their own fields
;  within the PARINFO structure, which will be ignored by MPFIT.
;  
;  PARINFO Example:
;  parinfo = replicate({value:0.D, fixed:0, limited:[0,0], $
;                       limits:[0.D,0]}, 5)
;  parinfo[0].fixed = 1
;  parinfo[4].limited[0] = 1
;  parinfo[4].limits[0]  = 50.D
;  parinfo[*].value = [5.7D, 2.2, 500., 1.5, 2000.]
;  
;  A total of 5 parameters, with starting values of 5.7,
;  2.2, 500, 1.5, and 2000 are given.  The first parameter
;  is fixed at a value of 5.7, and the last parameter is
;  constrained to be above 50.
;
;
; RECURSION
;
;  Generally, recursion is not allowed.  As of version 1.77, MPFIT has
;  recursion protection which does not allow a model function to
;  itself call MPFIT.  Users who wish to perform multi-level
;  optimization should investigate the 'EXTERNAL' function evaluation
;  methods described below for hard-to-evaluate functions.  That
;  method places more control in the user's hands.  The user can
;  design a "recursive" application by taking care.
;
;  In most cases the recursion protection should be well-behaved.
;  However, if the user is doing debugging, it is possible for the
;  protection system to get "stuck."  In order to reset it, run the
;  procedure:
;     MPFIT_RESET_RECURSION
;  and the protection system should get "unstuck."  It is save to call
;  this procedure at any time.
;
;
; COMPATIBILITY
;
;  This function is designed to work with IDL 5.0 or greater.
;  
;  Because TIED parameters and the "(EXTERNAL)" user-model feature use
;  the EXECUTE() function, they cannot be used with the free version
;  of the IDL Virtual Machine.
;
;
; DETERMINING THE VERSION OF MPFIT
;
;  MPFIT is a changing library.  Users of MPFIT may also depend on a
;  specific version of the library being present.  As of version 1.70
;  of MPFIT, a VERSION keyword has been added which allows the user to
;  query which version is present.  The keyword works like this:
;
;    RESULT = MPFIT(/query, VERSION=version)
;
;  This call uses the /QUERY keyword to query the version number
;  without performing any computations.  Users of MPFIT can call this
;  method to determine which version is in the IDL path before
;  actually using MPFIT to do any numerical work.  Upon return, the
;  VERSION keyword contains the version number of MPFIT, expressed as
;  a string of the form 'X.Y' where X and Y are integers.
;
;  Users can perform their own version checking, or use the built-in
;  error checking of MPFIT.  The MIN_VERSION keyword enforces the
;  requested minimum version number.  For example,
;
;    RESULT = MPFIT(/query, VERSION=version, MIN_VERSION='1.70')
;
;  will check whether the accessed version is 1.70 or greater, without
;  performing any numerical processing.
;
;  The VERSION and MIN_VERSION keywords were added in MPFIT
;  version 1.70 and later.  If the caller attempts to use the VERSION
;  or MIN_VERSION keywords, and an *older* version of the code is
;  present in the caller's path, then IDL will throw an 'unknown
;  keyword' error.  Therefore, in order to be robust, the caller, must
;  use exception handling.  Here is an example demanding at least
;  version 1.70.
;
;    MPFIT_OK = 0  & VERSION = '<unknown>'
;    CATCH, CATCHERR
;    IF CATCHERR EQ 0 THEN MPFIT_OK = MPFIT(/query, VERSION=version, $
;                                         MIN_VERSION='1.70')
;    CATCH, /CANCEL
;
;    IF NOT MPFIT_OK THEN $
;      MESSAGE, 'ERROR: you must have MPFIT version 1.70 or higher in '+$
;             'your path (found version '+version+')'
;
;  Of course, the caller can also do its own version number
;  requirements checking.
;
;
; HARD-TO-COMPUTE FUNCTIONS: "EXTERNAL" EVALUATION
;
;  The normal mode of operation for MPFIT is for the user to pass a
;  function name, and MPFIT will call the user function multiple times
;  as it iterates toward a solution.
;
;  Some user functions are particularly hard to compute using the
;  standard model of MPFIT.  Usually these are functions that depend
;  on a large amount of external data, and so it is not feasible, or
;  at least highly impractical, to have MPFIT call it.  In those cases
;  it may be possible to use the "(EXTERNAL)" evaluation option.
;
;  In this case the user is responsible for making all function *and
;  derivative* evaluations.  The function and Jacobian data are passed
;  in through the EXTERNAL_FVEC and EXTERNAL_FJAC keywords,
;  respectively.  The user indicates the selection of this option by
;  specifying a function name (MYFUNCT) of "(EXTERNAL)".  No
;  user-function calls are made when EXTERNAL evaluation is being
;  used.
;
;  ** SPECIAL NOTE ** For the "(EXTERNAL)" case, the quirk noted above
;     does not apply.  The gradient matrix, EXTERNAL_FJAC, should be
;     comparable to "-FGRAD(x,p)/err", which is the *opposite* sign of
;     the DP matrix described above.  In other words, EXTERNAL_FJAC
;     has the same sign as the derivative of EXTERNAL_FVEC, and the
;     opposite sign of FGRAD.
;
;  At the end of each iteration, control returns to the user, who must
;  reevaluate the function at its new parameter values.  Users should
;  check the return value of the STATUS keyword, where a value of 9
;  indicates the user should supply more data for the next iteration,
;  and re-call MPFIT.  The user may refrain from calling MPFIT
;  further; as usual, STATUS will indicate when the solution has
;  converged and no more iterations are required.
;
;  Because MPFIT must maintain its own data structures between calls,
;  the user must also pass a named variable to the EXTERNAL_STATE
;  keyword.  This variable must be maintained by the user, but not
;  changed, throughout the fitting process.  When no more iterations
;  are desired, the named variable may be discarded.
;
;
; INPUTS:
;   MYFUNCT - a string variable containing the name of the function to
;             be minimized.  The function should return the weighted
;             deviations between the model and the data, as described
;             above.
;
;             For EXTERNAL evaluation of functions, this parameter
;             should be set to a value of "(EXTERNAL)".
;
;   START_PARAMS - An one-dimensional array of starting values for each of the
;                  parameters of the model.  The number of parameters
;                  should be fewer than the number of measurements.
;                  Also, the parameters should have the same data type
;                  as the measurements (double is preferred).
;
;                  This parameter is optional if the PARINFO keyword
;                  is used (but see PARINFO).  The PARINFO keyword
;                  provides a mechanism to fix or constrain individual
;                  parameters.  If both START_PARAMS and PARINFO are
;                  passed, then the starting *value* is taken from
;                  START_PARAMS, but the *constraints* are taken from
;                  PARINFO.
; 
; RETURNS:
;
;   Returns the array of best-fit parameters.
;   Exceptions: 
;      * if /QUERY is set (see QUERY).
;
;
; KEYWORD PARAMETERS:
;
;   AUTODERIVATIVE - If this is set, derivatives of the function will
;                    be computed automatically via a finite
;                    differencing procedure.  If not set, then MYFUNCT
;                    must provide the explicit derivatives.
;                    Default: set (=1) 
;                    NOTE: to supply your own explicit derivatives,
;                      explicitly pass AUTODERIVATIVE=0
;
;   BESTNORM - upon return, the value of the summed squared weighted
;              residuals for the returned parameter values,
;              i.e. TOTAL(DEVIATES^2).
;
;   BEST_FJAC - upon return, BEST_FJAC contains the Jacobian, or
;               partial derivative, matrix for the best-fit model.
;               The values are an array,
;               ARRAY(N_ELEMENTS(DEVIATES),NFREE) where NFREE is the
;               number of free parameters.  This array is only
;               computed if /CALC_FJAC is set, otherwise BEST_FJAC is
;               undefined.
;
;               The returned array is such that BEST_FJAC[I,J] is the
;               partial derivative of DEVIATES[I] with respect to
;               parameter PARMS[PFREE_INDEX[J]].  Note that since
;               deviates are (data-model)*weight, the Jacobian of the
;               *deviates* will have the opposite sign from the
;               Jacobian of the *model*, and may be scaled by a
;               factor.
;
;   BEST_RESID - upon return, an array of best-fit deviates.
;
;   CALC_FJAC - if set, then calculate the Jacobian and return it in
;               BEST_FJAC.  If not set, then the return value of
;               BEST_FJAC is undefined.
;
;   COVAR - the covariance matrix for the set of parameters returned
;           by MPFIT.  The matrix is NxN where N is the number of
;           parameters.  The square root of the diagonal elements
;           gives the formal 1-sigma statistical errors on the
;           parameters IF errors were treated "properly" in MYFUNC.
;           Parameter errors are also returned in PERROR.
;
;           To compute the correlation matrix, PCOR, use this example:
;                  PCOR = COV * 0
;                  FOR i = 0, n-1 DO FOR j = 0, n-1 DO $
;                    PCOR[i,j] = COV[i,j]/sqrt(COV[i,i]*COV[j,j])
;           or equivalently, in vector notation,
;                  PCOR = COV / (PERROR # PERROR)
;
;           If NOCOVAR is set or MPFIT terminated abnormally, then
;           COVAR is set to a scalar with value !VALUES.D_NAN.
;
;   DOF - number of degrees of freedom, computed as
;             DOF = N_ELEMENTS(DEVIATES) - NFREE
;         Note that this doesn't account for pegged parameters (see
;         NPEGGED).  It also does not account for data points which
;         are assigned zero weight by the user function.
;
;   ERRMSG - a string error or warning message is returned.
;
;   EXTERNAL_FVEC - upon input, the function values, evaluated at
;                   START_PARAMS.  This should be an M-vector, where M
;                   is the number of data points.
;
;   EXTERNAL_FJAC - upon input, the Jacobian array of partial
;                   derivative values.  This should be a M x N array,
;                   where M is the number of data points and N is the
;                   number of parameters.  NOTE: that all FIXED or
;                   TIED parameters must *not* be included in this
;                   array.
;
;   EXTERNAL_STATE - a named variable to store MPFIT-related state
;                    information between iterations (used in input and
;                    output to MPFIT).  The user must not manipulate
;                    or discard this data until the final iteration is
;                    performed.
;
;   FASTNORM - set this keyword to select a faster algorithm to
;              compute sum-of-square values internally.  For systems
;              with large numbers of data points, the standard
;              algorithm can become prohibitively slow because it
;              cannot be vectorized well.  By setting this keyword,
;              MPFIT will run faster, but it will be more prone to
;              floating point overflows and underflows.  Thus, setting
;              this keyword may sacrifice some stability in the
;              fitting process.
;              
;   FTOL - a nonnegative input variable. Termination occurs when both
;          the actual and predicted relative reductions in the sum of
;          squares are at most FTOL (and STATUS is accordingly set to
;          1 or 3).  Therefore, FTOL measures the relative error
;          desired in the sum of squares.  Default: 1D-10
;
;   FUNCTARGS - A structure which contains the parameters to be passed
;               to the user-supplied function specified by MYFUNCT via
;               the _EXTRA mechanism.  This is the way you can pass
;               additional data to your user-supplied function without
;               using common blocks.
;
;               Consider the following example:
;                if FUNCTARGS = { XVAL:[1.D,2,3], YVAL:[1.D,4,9],
;                                 ERRVAL:[1.D,1,1] }
;                then the user supplied function should be declared
;                like this:
;                FUNCTION MYFUNCT, P, XVAL=x, YVAL=y, ERRVAL=err
;
;               By default, no extra parameters are passed to the
;               user-supplied function, but your function should
;               accept *at least* one keyword parameter.  [ This is to
;               accomodate a limitation in IDL's _EXTRA
;               parameter-passing mechanism. ]
;
;   GTOL - a nonnegative input variable. Termination occurs when the
;          cosine of the angle between fvec and any column of the
;          jacobian is at most GTOL in absolute value (and STATUS is
;          accordingly set to 4). Therefore, GTOL measures the
;          orthogonality desired between the function vector and the
;          columns of the jacobian.  Default: 1D-10
;
;   ITERARGS - The keyword arguments to be passed to ITERPROC via the
;              _EXTRA mechanism.  This should be a structure, and is
;              similar in operation to FUNCTARGS.
;              Default: no arguments are passed.
;
;   ITERPRINT - The name of an IDL procedure, equivalent to PRINT,
;               that ITERPROC will use to render output.  ITERPRINT
;               should be able to accept at least four positional
;               arguments.  In addition, it should be able to accept
;               the standard FORMAT keyword for output formatting; and
;               the UNIT keyword, to redirect output to a logical file
;               unit (default should be UNIT=1, standard output).
;               These keywords are passed using the ITERARGS keyword
;               above.  The ITERPRINT procedure must accept the _EXTRA
;               keyword.  
;               NOTE: that much formatting can be handled with the 
;                     MPPRINT and MPFORMAT tags.
;               Default: 'MPFIT_DEFPRINT' (default internal formatter)
;
;   ITERPROC - The name of a procedure to be called upon each NPRINT
;              iteration of the MPFIT routine.  ITERPROC is always
;              called in the final iteration.  It should be declared
;              in the following way:
;
;              PRO ITERPROC, MYFUNCT, p, iter, fnorm, FUNCTARGS=fcnargs, $
;                PARINFO=parinfo, QUIET=quiet, DOF=dof, PFORMAT=pformat, $
;                UNIT=unit, ...
;                ; perform custom iteration update
;              END
;         
;              ITERPROC must either accept all three keyword
;              parameters (FUNCTARGS, PARINFO and QUIET), or at least
;              accept them via the _EXTRA keyword.
;          
;              MYFUNCT is the user-supplied function to be minimized,
;              P is the current set of model parameters, ITER is the
;              iteration number, and FUNCTARGS are the arguments to be
;              passed to MYFUNCT.  FNORM should be the chi-squared
;              value.  QUIET is set when no textual output should be
;              printed.  DOF is the number of degrees of freedom,
;              normally the number of points less the number of free
;              parameters.  See below for documentation of PARINFO.
;              PFORMAT is the default parameter value format.  UNIT is
;              passed on to the ITERPRINT procedure, and should
;              indicate the file unit where log output will be sent
;              (default: standard output).
;
;              In implementation, ITERPROC can perform updates to the
;              terminal or graphical user interface, to provide
;              feedback while the fit proceeds.  If the fit is to be
;              stopped for any reason, then ITERPROC should set the
;              common block variable ERROR_CODE to negative value
;              between -15 and -1 (see MPFIT_ERROR common block
;              below).  In principle, ITERPROC should probably not
;              modify the parameter values, because it may interfere
;              with the algorithm's stability.  In practice it is
;              allowed.
;
;              Default: an internal routine is used to print the
;                       parameter values.
;
;   ITERSTOP - Set this keyword if you wish to be able to stop the
;              fitting by hitting the predefined ITERKEYSTOP key on
;              the keyboard.  This only works if you use the default
;              ITERPROC.
;
;   ITERKEYSTOP - A keyboard key which will halt the fit (and if
;                 ITERSTOP is set and the default ITERPROC is used).
;                 ITERSTOPKEY may either be a one-character string
;                 with the desired key, or a scalar integer giving the
;                 ASCII code of the desired key.  
;                 Default: 7b (control-g)
;
;                 NOTE: the default value of ASCI 7 (control-G) cannot
;                 be read in some windowing environments, so you must
;                 change to a printable character like 'q'.
;
;   MAXITER - The maximum number of iterations to perform.  If the
;             number of calculation iterations exceeds MAXITER, then
;             the STATUS value is set to 5 and MPFIT returns.  
;
;             If MAXITER EQ 0, then MPFIT does not iterate to adjust
;             parameter values; however, the user function is evaluated
;             and parameter errors/covariance/Jacobian are estimated
;             before returning.
;             Default: 200 iterations
;
;   MIN_VERSION - The minimum requested version number.  This must be
;                 a scalar string of the form returned by the VERSION
;                 keyword.  If the current version of MPFIT does not
;                 satisfy the minimum requested version number, then,
;                    MPFIT(/query, min_version='...') returns 0
;                    MPFIT(...) returns NAN
;                 Default: no version number check
;                 NOTE: MIN_VERSION was added in MPFIT version 1.70
;
;   NFEV - the number of MYFUNCT function evaluations performed.
;
;   NFREE - the number of free parameters in the fit.  This includes
;           parameters which are not FIXED and not TIED, but it does
;           include parameters which are pegged at LIMITS.
;
;   NITER - the number of iterations completed.
;
;   NOCATCH - if set, then MPFIT will not perform any error trapping.
;             By default (not set), MPFIT will trap errors and report
;             them to the caller.  This keyword will typically be used
;             for debugging.
;
;   NOCOVAR - set this keyword to prevent the calculation of the
;             covariance matrix before returning (see COVAR)
;
;   NPEGGED - the number of free parameters which are pegged at a
;             LIMIT.
;
;   NPRINT - The frequency with which ITERPROC is called.  A value of
;            1 indicates that ITERPROC is called with every iteration,
;            while 2 indicates every other iteration, etc.  Be aware
;            that several Levenberg-Marquardt attempts can be made in
;            a single iteration.  Also, the ITERPROC is *always*
;            called for the final iteration, regardless of the
;            iteration number.
;            Default value: 1
;
;   PARINFO - A one-dimensional array of structures.
;             Provides a mechanism for more sophisticated constraints
;             to be placed on parameter values.  When PARINFO is not
;             passed, then it is assumed that all parameters are free
;             and unconstrained.  Values in PARINFO are never 
;             modified during a call to MPFIT.
;
;             See description above for the structure of PARINFO.
;
;             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
;
;   PFREE_INDEX - upon return, PFREE_INDEX contains an index array
;                 which indicates which parameter were allowed to
;                 vary.  I.e. of all the parameters PARMS, only
;                 PARMS[PFREE_INDEX] were varied.
;
;   QUERY - if set, then MPFIT() will return immediately with one of
;           the following values:
;                 1 - if MIN_VERSION is not set
;                 1 - if MIN_VERSION is set and MPFIT satisfies the minimum
;                 0 - if MIN_VERSION is set and MPFIT does not satisfy it
;           The VERSION output keyword is always set upon return.
;           Default: not set.
;
;   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.
;
;   VERSION - upon return, VERSION will be set to the MPFIT internal
;             version number.  The version number will be a string of
;             the form "X.Y" where X is a major revision number and Y
;             is a minor revision number.
;             NOTE: the VERSION keyword was not present before 
;               MPFIT version number 1.70, therefore, callers must 
;               use exception handling when using this keyword.
;
;   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