mpfit.pro

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

;  
;  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.
;
;
; 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.
;
;
; 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.
;
;
; 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 - the value of the summed squared weighted residuals for
;              the returned parameter values, i.e. TOTAL(DEVIATES^2).
;
;   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:
;           IDL> PCOR = COV * 0
;           IDL> 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])
;
;           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).
;
;   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 is exceeded, then the STATUS value is set to 5
;             and MPFIT returns.
;             Default: 200 iterations
;
;   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.