mpfitfun.pro

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

;             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
;
;   QUIET - set this keyword when no textual output should be printed
;           by MPFIT
;
;   STATUS - an integer status code is returned.  All values other
;            than zero can represent success.  It can have one of the
;            following values:
;
;	   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.
;
;   WEIGHTS - Array of weights to be used in calculating the
;             chi-squared value.  If WEIGHTS is specified then the ERR
;             parameter is ignored.  The chi-squared value is computed
;             as follows:
;
;                CHISQ = TOTAL( (Y-MYFUNCT(X,P))^2 * ABS(WEIGHTS) )
;
;             Here are common values of WEIGHTS:
;
;                1D/ERR^2 - Normal weighting (ERR is the measurement error)
;                1D/Y     - Poisson weighting (counting statistics)
;                1D       - Unweighted
;
;   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
;
;   YFIT - the best-fit model function, as returned by MYFUNCT.
;
;   
; EXAMPLE:
;
;   ; First, generate some synthetic data
;   npts = 200
;   x  = dindgen(npts) * 0.1 - 10.                  ; Independent variable 
;   yi = gauss1(x, [2.2D, 1.4, 3000.])              ; "Ideal" Y variable
;   y  = yi + randomn(seed, npts) * sqrt(1000. + yi); Measured, w/ noise
;   sy = sqrt(1000.D + y)                           ; Poisson errors
;
;   ; Now fit a Gaussian to see how well we can recover
;   p0 = [1.D, 1., 1000.]                   ; Initial guess (cent, width, area)
;   p = mpfitfun('GAUSS1', x, y, sy, p0)    ; Fit a function
;   print, p
;
;   Generates a synthetic data set with a Gaussian peak, and Poisson
;   statistical uncertainty.  Then the same function is fitted to the
;   data (with different starting parameters) to see how close we can
;   get.
;
;
; 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.  By default the value of ERROR_CODE is
;     zero, indicating a successful function/procedure call.
;
; 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.
;
; MODIFICATION HISTORY:
;   Written, Apr-Jul 1998, CM
;   Added PERROR keyword, 04 Aug 1998, CM
;   Added COVAR keyword, 20 Aug 1998, CM
;   Added ITER output keyword, 05 Oct 1998
;      D.L Windt, Bell Labs, [e-mail];
;   Added ability to return model function in YFIT, 09 Nov 1998
;   Analytical derivatives allowed via AUTODERIVATIVE keyword, 09 Nov 1998
;   Parameter values can be tied to others, 09 Nov 1998
;   Cosmetic documentation updates, 16 Apr 1999, CM
;   More cosmetic documentation updates, 14 May 1999, CM
;   Made sure to update STATUS, 25 Sep 1999, CM
;   Added WEIGHTS keyword, 25 Sep 1999, CM
;   Changed from handles to common blocks, 25 Sep 1999, CM
;     - commons seem much cleaner and more logical in this case.
;   Alphabetized documented keywords, 02 Oct 1999, CM
;   Added QUERY keyword and query checking of MPFIT, 29 Oct 1999, CM
;   Corrected EXAMPLE (offset of 1000), 30 Oct 1999, CM
;   Check to be sure that X and Y are present, 02 Nov 1999, CM
;   Documented PERROR for unweighted fits, 03 Nov 1999, CM
;   Changed to ERROR_CODE for error condition, 28 Jan 2000, CM
;   Corrected errors in EXAMPLE, 26 Mar 2000, CM
;   Copying permission terms have been liberalized, 26 Mar 2000, CM
;   Propagated improvements from MPFIT, 17 Dec 2000, CM
;   Added CASH statistic, 10 Jan 2001
;   Added NFREE and NPEGGED keywords, 11 Sep 2002, CM
;   Documented RELSTEP field of PARINFO (!!), CM, 25 Oct 2002
;   Add DOF keyword to return degrees of freedom, CM, 23 June 2003
;   Convert to IDL 5 array syntax (!), 16 Jul 2006, CM
;   Move STRICTARR compile option inside each function/procedure, 9
;     Oct 2006
;   Add NAN keyword, to ignore non-finite data values, 28 Oct 2006, CM
;   Clarify documentation on user-function, derivatives, and PARINFO,
;     27 May 2007
;   Fix bug in handling of explicit derivatives with errors/weights
;     (the weights were not being applied), CM, 03 Sep 2007
;   Add COMPATIBILITY section, CM, 13 Dec 2007
;
;  $Id: mpfitfun.pro,v 1.13 2007/12/15 14:28:27 craigm Exp $
;-
; Copyright (C) 1997-2002, 2003, 2006, 2007, Craig Markwardt
; This software is provided as is without any warranty whatsoever.
; Permission to use, copy, modify, and distribute modified or
; unmodified copies is granted, provided this copyright and disclaimer
; are included unchanged.
;-

FORWARD_FUNCTION mpfitfun_eval, mpfitfun, mpfit

; This is the call-back function for MPFIT.  It evaluates the
; function, subtracts the data, and returns the residuals.
function mpfitfun_eval, p, dp, _EXTRA=extra

  COMPILE_OPT strictarr
  common mpfitfun_common, fcn, x, y, err, wts, f, fcnargs

  ;; Save the original DP matrix for later use
  if n_params() GT 1 then if n_elements(dp) GT 0 then dp0 = dp

  ;; The function is evaluated here.  There are four choices,
  ;; depending on whether (a) FUNCTARGS was passed to MPFITFUN, which
  ;; is passed to this function as "hf"; or (b) the derivative
  ;; parameter "dp" is passed, meaning that derivatives should be
  ;; calculated analytically by the function itself.
  if n_elements(fcnargs) GT 0 then begin
      if n_params() GT 1 then f = call_function(fcn, x, p, dp, _EXTRA=fcnargs)$
      else                    f = call_function(fcn, x, p, _EXTRA=fcnargs)
  endif else begin
      if n_params() GT 1 then f = call_function(fcn, x, p, dp) $
      else                    f = call_function(fcn, x, p)
  endelse

  np = n_elements(p)
  nf = n_elements(f)

  ;; Compute the deviates, applying either errors or weights
  if n_elements(wts) GT 0 then begin
      result = (y-f)*wts
      if n_elements(dp0) GT 0 AND n_elements(dp) EQ np*nf then begin
          for j = 0, np-1 do dp[j*nf] = dp[j*nf:j*nf+nf-1] * wts
      endif
  endif else if n_elements(err) GT 0 then begin
      result = (y-f)/err
      if n_elements(dp0) GT 0 AND n_elements(dp) EQ np*nf then begin
          for j = 0, np-1 do dp[j*nf] = dp[j*nf:j*nf+nf-1] / err
      endif
  endif else begin
      result = (y-f)
  endelse
      
  ;; Make sure the returned result is one-dimensional.
  result = reform(result, n_elements(result), /overwrite)
  return, result
  
end

;; Implement residual and gradient scaling according to the
;; prescription of Cash (ApJ, 228, 939)
pro mpfitfun_cash, resid, dresid
  COMPILE_OPT strictarr
  common mpfitfun_common, fcn, x, y, err, wts, f, fcnargs

  sz = size(dresid)
  m = sz[1]
  n = sz[2]

  ;; Do rudimentary dimensions checks, so we don't do something stupid
  if n_elements(y) NE m OR n_elements(f) NE m OR n_elements(resid) NE m then begin
      DIM_ERROR:
      message, 'ERROR: dimensions of Y, F, RESID or DRESID are not consistent'
  endif

  ;; Scale gradient by sqrt(y)/f
  gfact = temporary(dresid) * rebin(reform(sqrt(y)/f,m,1),m,n)
  dresid = reform(dresid, m, n, /overwrite)
  
  ;; Scale residuals by 1/sqrt(y)
  resid = temporary(resid)/sqrt(y)

  return
end

function mpfitfun, fcn, x, y, err, p, WEIGHTS=wts, FUNCTARGS=fa, $
                   BESTNORM=bestnorm, nfev=nfev, STATUS=status, $
                   parinfo=parinfo, query=query, CASH=cash, $
                   covar=covar, perror=perror, yfit=yfit, $
                   niter=niter, nfree=nfree, npegged=npegged, dof=dof, $
                   quiet=quiet, ERRMSG=errmsg, NAN=NAN, _EXTRA=extra

  COMPILE_OPT strictarr
  status = 0L
  errmsg = ''

  ;; Detect MPFIT and crash if it was not found
  catch, catcherror
  if catcherror NE 0 then begin
      MPFIT_NOTFOUND:
      catch, /cancel
      message, 'ERROR: the required function MPFIT must be in your IDL path', /info
      return, !values.d_nan
  endif
  if mpfit(/query) NE 1 then goto, MPFIT_NOTFOUND
  catch, /cancel
  if keyword_set(query) then return, 1

  if n_params() EQ 0 then begin
      message, "USAGE: PARMS = MPFITFUN('MYFUNCT', X, Y, ERR, "+ $
        "START_PARAMS, ... )", /info
      return, !values.d_nan
  endif
  if n_elements(x) EQ 0 OR n_elements(y) EQ 0 then begin
      message, 'ERROR: X and Y must be defined', /info
      return, !values.d_nan
  endif

  if n_elements(err) GT 0 OR n_elements(wts) GT 0 AND keyword_set(cash) then begin
      message, 'ERROR: WEIGHTS or ERROR cannot be specified with CASH', /info
      return, !values.d_nan
  endif
  if keyword_set(cash) then begin
      scalfcn = 'mpfitfun_cash'
  endif

  ;; Use common block to pass data back and forth
  common mpfitfun_common, fc, xc, yc, ec, wc, mc, ac
  fc = fcn & xc = x & yc = y & mc = 0L
  ;; These optional parameters must be undefined first
  ac = 0 & dummy = size(temporary(ac))
  ec = 0 & dummy = size(temporary(ec))
  wc = 0 & dummy = size(temporary(wc))

  ;; FUNCTARGS
  if n_elements(fa) GT 0 then ac = fa

  ;; WEIGHTS or ERROR
  if n_elements(wts) GT 0 then begin
      wc = sqrt(abs(wts))
  endif else if n_elements(err) GT 0 then begin
      wh = where(err EQ 0, ct)
      if ct GT 0 then begin
          errmsg = 'ERROR: ERROR value must not be zero.  Use WEIGHTS instead.'
          message, errmsg, /info
          return, !values.d_nan
      endif
      ;; Appropriate weight for gaussian errors
      wc = 1/err
  endif

  ;; Check for weights/errors which do not match the dimension
  ;; of the data points
  if n_elements(wc) GT 0 AND $
    n_elements(wc) NE 1 AND $
    n_elements(wc) NE n_elements(yc) then begin
      errmsg = 'ERROR: ERROR/WEIGHTS must either be a scalar or match the number of Y values'
      message, errmsg, /info
      return, !values.d_nan
  endif

  ;; If the weights/errors are a scalar value, and not finite, then 
  ;; the fit will surely fail
  if n_elements(wc) EQ 1 then begin
      if finite(wc[0]) EQ 0 then begin
          errmsg = 'ERROR: the supplied scalar WEIGHT/ERROR value was not finite'
          message, errmsg, /info
          return, !values.d_nan
      endif
  endif

  ;; Handle the cases of non-finite data points or weights
  if keyword_set(nan) then begin
      ;; Non-finite data points
      wh = where(finite(yc) EQ 0, ct)
      if ct GT 0 then begin
          yc[wh] = 0
          ;; Careful: handle case when weights were a scalar...
          ;;   ... promote to a vector
          if n_elements(wc) EQ 1 then wc = replicate(wc[0], n_elements(yc))
          wc[wh] = 0
      endif

      ;; Non-finite weights
      wh = where(finite(wc) EQ 0, ct)
      if ct GT 0 then wc[wh] = 0
  endif

  result = mpfit('mpfitfun_eval', p, SCALE_FCN=scalfcn, $
                 parinfo=parinfo, STATUS=status, nfev=nfev, BESTNORM=bestnorm,$
                 covar=covar, perror=perror, $
                 niter=niter, nfree=nfree, npegged=npegged, dof=dof, $
                 ERRMSG=errmsg, quiet=quiet, _EXTRA=extra)

  ;; Retrieve the fit value
  yfit = temporary(mc)
  ;; Some cleanup
  xc = 0 & yc = 0 & wc = 0 & ec = 0 & mc = 0 & ac = 0

  ;; Print error message if there is one.
  if NOT keyword_set(quiet) AND errmsg NE '' then $
    message, errmsg, /info

  return, result
end