mpfit.pro

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

;   Added checks for parameter and function overflow; a new STATUS
;     value to reflect this; STATUS values of -15 to -1 are reserved
;     for user function errors, CM, 03 Apr 2001
;   DAMP keyword is now a TANH, CM, 03 Apr 2001
;   Added more error checking of float vs. double, CM, 07 Apr 2001
;   Fixed bug in handling of parameter lower limits; moved overflow
;     checking to end of loop, CM, 20 Apr 2001
;   Failure using GOTO, TERMINATE more graceful if FNORM1 not defined,
;     CM, 13 Aug 2001
;   Add MPPRINT tag to PARINFO, CM, 19 Nov 2001
;   Add DOF keyword to DEFITER procedure, and print degrees of
;     freedom, CM, 28 Nov 2001
;   Add check to be sure MYFUNCT is a scalar string, CM, 14 Jan 2002
;   Addition of EXTERNAL_FJAC, EXTERNAL_FVEC keywords; ability to save
;     fitter's state from one call to the next; allow '(EXTERNAL)'
;     function name, which implies that user will supply function and
;     Jacobian at each iteration, CM, 10 Mar 2002
;   Documented EXTERNAL evaluation code, CM, 10 Mar 2002
;   Corrected signficant bug in the way that the STEP parameter, and
;     FIXED parameters interacted (Thanks Andrew Steffl), CM, 02 Apr
;     2002
;   Allow COVAR and PERROR keywords to be computed, even in case of
;     '(EXTERNAL)' function, 26 May 2002
;   Add NFREE and NPEGGED keywords; compute NPEGGED; compute DOF using
;     NFREE instead of n_elements(X), thanks to Kristian Kjaer, CM 11
;     Sep 2002
;   Hopefully PERROR is all positive now, CM 13 Sep 2002
;   Documented RELSTEP field of PARINFO (!!), CM, 25 Oct 2002
;   Error checking to detect missing start pars, CM 12 Apr 2003
;   Add DOF keyword to return degrees of freedom, CM, 30 June 2003
;   Always call ITERPROC in the final iteration; add ITERKEYSTOP
;     keyword, CM, 30 June 2003
;   Correct bug in MPFIT_LMPAR of singularity handling, which might
;     likely be fatal for one-parameter fits, CM, 21 Nov 2003
;     (with thanks to Peter Tuthill for the proper test case)
;   Minor documentation adjustment, 03 Feb 2004, CM
;   Correct small error in QR factorization when pivoting; document
;     the return values of QRFAC when pivoting, 21 May 2004, CM
;   Add MPFORMAT field to PARINFO, and correct behavior of interaction
;     between MPPRINT and PARNAME in MPFIT_DEFITERPROC (thanks to Tim
;     Robishaw), 23 May 2004, CM
;   Add the ITERPRINT keyword to allow redirecting output, 26 Sep
;     2004, CM
;   Correct MAXSTEP behavior in case of a negative parameter, 26 Sep
;     2004, CM
;   Fix bug in the parsing of MINSTEP/MAXSTEP, 10 Apr 2005, CM
;   Fix bug in the handling of upper/lower limits when the limit was
;     negative (the fitting code would never "stick" to the lower
;     limit), 29 Jun 2005, CM
;   Small documentation update for the TIED field, 05 Sep 2005, CM
;   Convert to IDL 5 array syntax (!), 16 Jul 2006, CM
;   If MAXITER equals zero, then do the basic parameter checking and
;     uncertainty analysis, but do not adjust the parameters, 15 Aug
;     2006, CM
;   Added documentation, 18 Sep 2006, CM
;   A few more IDL 5 array syntax changes, 25 Sep 2006, CM
;   Move STRICTARR compile option inside each function/procedure, 9 Oct 2006
;   Bug fix for case of MPMAXSTEP and fixed parameters, thanks
;     to Huib Intema (who found it from the Python translation!), 05 Feb 2007
;   Similar fix for MPFIT_FDJAC2 and the MPSIDE sidedness of
;     derivatives, also thanks to Huib Intema, 07 Feb 2007
;   Clarify documentation on user-function, derivatives, and PARINFO,
;     27 May 2007
;   Change the wording of "Analytic Derivatives" to "Explicit 
;     Derivatives" in the documentation, CM, 03 Sep 2007
;   Further documentation tweaks, CM, 13 Dec 2007
;   Add COMPATIBILITY section and add credits to copyright, CM, 13 Dec
;      2007
;   Document and enforce that START_PARMS and PARINFO are 1-d arrays,
;      CM, 29 Mar 2008
;   Previous change for 1-D arrays wasn't correct for
;      PARINFO.LIMITED/.LIMITS; now fixed, CM, 03 May 2008
;   Documentation adjustments, CM, 20 Aug 2008
;   Change some minor FOR-loop variables to type-long, CM, 03 Sep 2008
;   Change error handling slightly, document NOCATCH keyword,
;      document error handling in general, CM, 01 Oct 2008
;   Special case: when either LIMITS is zero, and a parameter pushes
;      against that limit, the coded that 'pegged' it there would not
;      work since it was a relative condition; now zero is handled
;      properly, CM, 08 Nov 2008
;
;  $Id: mpfit.pro,v 1.59 2008/11/08 16:11:54 craigm Exp $
;-
; Original MINPACK by More' Garbow and Hillstrom, translated with permission
; Modifications and enhancements are:
; Copyright (C) 1997-2003, 2004, 2005, 2006, 2007, 2008, 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.
;-

pro mpfit_dummy
  ;; Enclose in a procedure so these are not defined in the main level
  COMPILE_OPT strictarr
  FORWARD_FUNCTION mpfit_fdjac2, mpfit_enorm, mpfit_lmpar, mpfit_covar, $
    mpfit, mpfit_call

  COMMON mpfit_error, error_code  ;; For error passing to user function
  COMMON mpfit_config, mpconfig   ;; For internal error configrations
end

;; Reset profiling registers for another run.  By default, and when
;; uncommented, the profiling registers simply accumulate.

pro mpfit_resetprof
  COMPILE_OPT strictarr
  common mpfit_profile, mpfit_profile_vals

  mpfit_profile_vals = { status: 1L, fdjac2: 0D, lmpar: 0D, mpfit: 0D, $
                         qrfac: 0D,  qrsolv: 0D, enorm: 0D}
  return
end

;; Following are machine constants that can be loaded once.  I have
;; found that bizarre underflow messages can be produced in each call
;; to MACHAR(), so this structure minimizes the number of calls to
;; one.

pro mpfit_setmachar, double=isdouble
  COMPILE_OPT strictarr
  common mpfit_profile, profvals
  if n_elements(profvals) EQ 0 then mpfit_resetprof

  common mpfit_machar, mpfit_machar_vals

  ;; In earlier versions of IDL, MACHAR itself could produce a load of
  ;; error messages.  We try to mask some of that out here.
  if (!version.release) LT 5 then dummy = check_math(1, 1)

  mch = 0.
  mch = machar(double=keyword_set(isdouble))
  dmachep = mch.eps
  dmaxnum = mch.xmax
  dminnum = mch.xmin
  dmaxlog = alog(mch.xmax)
  dminlog = alog(mch.xmin)
  if keyword_set(isdouble) then $
    dmaxgam = 171.624376956302725D $
  else $
    dmaxgam = 171.624376956302725
  drdwarf = sqrt(dminnum*1.5) * 10
  drgiant = sqrt(dmaxnum) * 0.1

  mpfit_machar_vals = {machep: dmachep, maxnum: dmaxnum, minnum: dminnum, $
                       maxlog: dmaxlog, minlog: dminlog, maxgam: dmaxgam, $
                       rdwarf: drdwarf, rgiant: drgiant}

  if (!version.release) LT 5 then dummy = check_math(0, 0)

  return
end



;; Call user function or procedure, with _EXTRA or not, with
;; derivatives or not.
function mpfit_call, fcn, x, fjac, _EXTRA=extra

  COMPILE_OPT strictarr
  common mpfit_config, mpconfig

  if keyword_set(mpconfig.qanytied) then mpfit_tie, x, mpconfig.ptied
    
  ;; Decide whether we are calling a procedure or function
  if mpconfig.proc then proc = 1 else proc = 0
  mpconfig.nfev = mpconfig.nfev + 1

  if proc then begin
      if n_params() EQ 3 then begin
          if n_elements(extra) GT 0 then $
            call_procedure, fcn, x, f, fjac, _EXTRA=extra $
          else $
            call_procedure, fcn, x, f, fjac
      endif else begin
          if n_elements(extra) GT 0 then $
            call_procedure, fcn, x, f, _EXTRA=extra $
          else $
            call_procedure, fcn, x, f
      endelse
  endif else begin
      if n_params() EQ 3 then begin
          if n_elements(extra) GT 0 then $
            f = call_function(fcn, x, fjac, _EXTRA=extra) $
          else $
            f = call_function(fcn, x, fjac)
      endif else begin
          if n_elements(extra) GT 0 then $
            f = call_function(fcn, x, _EXTRA=extra) $
          else $
            f = call_function(fcn, x)
      endelse
  endelse  

  if n_params() EQ 2 AND mpconfig.damp GT 0 then begin
      damp = mpconfig.damp[0]
      
      ;; Apply the damping if requested.  This replaces the residuals
      ;; with their hyperbolic tangent.  Thus residuals larger than
      ;; DAMP are essentially clipped.
      f = tanh(f/damp)
  endif

  return, f
end

function mpfit_fdjac2, fcn, x, fvec, step, ulimited, ulimit, dside, $
                 iflag=iflag, epsfcn=epsfcn, autoderiv=autoderiv, $
                 FUNCTARGS=fcnargs, xall=xall, ifree=ifree, dstep=dstep

  COMPILE_OPT strictarr
  common mpfit_machar, machvals
  common mpfit_profile, profvals
  common mpfit_error, mperr

;  prof_start = systime(1)
  MACHEP0 = machvals.machep
  DWARF   = machvals.minnum

  if n_elements(epsfcn) EQ 0 then epsfcn = MACHEP0
  if n_elements(xall)   EQ 0 then xall = x
  if n_elements(ifree)  EQ 0 then ifree = lindgen(n_elements(xall))
  if n_elements(step)   EQ 0 then step = x * 0.
  nall = n_elements(xall)

  eps = sqrt(max([epsfcn, MACHEP0]));
  m = n_elements(fvec)
  n = n_elements(x)

  ;; Compute analytical derivative if requested
  if NOT keyword_set(autoderiv) then begin
      mperr = 0
      fjac = intarr(nall)
      fjac[ifree] = 1      ;; Specify which parameters need derivatives
      fp = mpfit_call(fcn, xall, fjac, _EXTRA=fcnargs)
      iflag = mperr

      if n_elements(fjac) NE m*nall then begin
          message, /cont, 'ERROR: Derivative matrix was not computed properly.'
          iflag = 1
;          profvals.fdjac2 = profvals.fdjac2 + (systime(1) - prof_start)
          return, 0
      endif

      ;; This definition is consistent with CURVEFIT
      ;; Sign error found (thanks Jesus Fernandez <[e-mail]>)

      ;; ... and now I regret doing this sign flip since it's not
      ;; strictly correct.  The definition should be RESID =
      ;; (Y-F)/SIGMA, so d(RESID)/dP should be -dF/dP.  My response to
      ;; Fernandez was unfounded because he was trying to supply
      ;; dF/dP.  Sigh. (CM 31 Aug 2007)

      fjac = reform(-temporary(fjac), m, nall, /overwrite)

      ;; Select only the free parameters
      if n_elements(ifree) LT nall then $
        fjac = reform(fjac[*,ifree], m, n, /overwrite)
;      profvals.fdjac2 = profvals.fdjac2 + (systime(1) - prof_start)
      return, fjac
  endif

  fjac = make_array(m, n, value=fvec[0]*0.)
  fjac = reform(fjac, m, n, /overwrite)

  h = eps * abs(x)

  ;; if STEP is given, use that
  ;; STEP includes the fixed parameters
  if n_elements(step) GT 0 then begin
      stepi = step[ifree]
      wh = where(stepi GT 0, ct)
      if ct GT 0 then h[wh] = stepi[wh]
  endif

  ;; if relative step is given, use that
  ;; DSTEP includes the fixed parameters
  if n_elements(dstep) GT 0 then begin
      dstepi = dstep[ifree]
      wh = where(dstepi GT 0, ct)
      if ct GT 0 then h[wh] = abs(dstepi[wh]*x[wh])
  endif

  ;; In case any of the step values are zero
  wh = where(h EQ 0, ct)
  if ct GT 0 then h[wh] = eps

  ;; Reverse the sign of the step if we are up against the parameter
  ;; limit, or if the user requested it.
  ;; DSIDE includes the fixed parameters (ULIMITED/ULIMIT have only
  ;; varying ones)
  mask = dside[ifree] EQ -1
  if n_elements(ulimited) GT 0 AND n_elements(ulimit) GT 0 then $
    mask = mask OR (ulimited AND (x GT ulimit-h))
  wh = where(mask, ct)
  if ct GT 0 then h[wh] = -h[wh]

  ;; Loop through parameters, computing the derivative for each
  for j=0L, n-1 do begin
      xp = xall
      xp[ifree[j]] = xp[ifree[j]] + h[j]
      
      mperr = 0
      fp = mpfit_call(fcn, xp, _EXTRA=fcnargs)
      
      iflag = mperr