mpfit.pro

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

;+
; NAME:
;   MPFIT
;
; AUTHOR:
;   Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
;   [e-mail]
;   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 where ever 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, if you desire to compute them
;  explicitly, then the AUTODERIVATIVE=0 keyword must be passed to
;  MPFIT, and then the model function must be declared with an
;  additional positional parameter as described below.  As a practical
;  matter, it is often sufficient and even faster to allow MPFIT to
;  calculate the derivatives numerically, and so AUTODERIVATIVE=0 is
;  not necessary.
;
;  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:
;
;    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).
;
;  When finite differencing is used for computing derivatives (ie,
;  when AUTODERIVATIVE=1), 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.
;  
; 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 - the sidedness of the finite difference when computing
;               numerical derivatives.  This field can take four
;               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)
;
;              Where H is the STEP parameter described above.  The
;              "automatic" one-sided derivative method will chose a
;              direction for the finite difference which does not
;              violate any constraints.  The other methods do not
;              perform this check.  The two-sided method is in
;              principle more precise, but requires twice as many
;              function evaluations.  Default: 0.
;
;     .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.
;             [ 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.