mpfit.pro

;     is zero, indicating a successful function/procedure call.
;
;   COMMON MPFIT_PROFILE
;   COMMON MPFIT_MACHAR
;   COMMON MPFIT_CONFIG
;
;     These are undocumented common blocks are used internally by
;     MPFIT and may change in future implementations.
;
; THEORY OF OPERATION:
;
;   There are many specific strategies for function minimization.  One
;   very popular technique is to use function gradient information to
;   realize the local structure of the function.  Near a local minimum
;   the function value can be taylor expanded about x0 as follows:
;
;      f(x) = f(x0) + f'(x0) . (x-x0) + (1/2) (x-x0) . f''(x0) . (x-x0)
;             -----   ---------------   -------------------------------  (1)
;     Order    0th          1st                      2nd
;
;   Here f'(x) is the gradient vector of f at x, and f''(x) is the
;   Hessian matrix of second derivatives of f at x.  The vector x is
;   the set of function parameters, not the measured data vector.  One
;   can find the minimum of f, f(xm) using Newton's method, and
;   arrives at the following linear equation:
;
;      f''(x0) . (xm-x0) = - f'(x0)                            (2)
;
;   If an inverse can be found for f''(x0) then one can solve for
;   (xm-x0), the step vector from the current position x0 to the new
;   projected minimum.  Here the problem has been linearized (ie, the
;   gradient information is known to first order).  f''(x0) is
;   symmetric n x n matrix, and should be positive definite.
;
;   The Levenberg - Marquardt technique is a variation on this theme.
;   It adds an additional diagonal term to the equation which may aid the
;   convergence properties:
;
;      (f''(x0) + nu I) . (xm-x0) = -f'(x0)                  (2a)
;
;   where I is the identity matrix.  When nu is large, the overall
;   matrix is diagonally dominant, and the iterations follow steepest
;   descent.  When nu is small, the iterations are quadratically
;   convergent.
;
;   In principle, if f''(x0) and f'(x0) are known then xm-x0 can be
;   determined.  However the Hessian matrix is often difficult or
;   impossible to compute.  The gradient f'(x0) may be easier to
;   compute, if even by finite difference techniques.  So-called
;   quasi-Newton techniques attempt to successively estimate f''(x0)
;   by building up gradient information as the iterations proceed.
;
;   In the least squares problem there are further simplifications
;   which assist in solving eqn (2).  The function to be minimized is
;   a sum of squares:
;
;       f = Sum(hi^2)                                         (3)
;
;   where hi is the ith residual out of m residuals as described
;   above.  This can be substituted back into eqn (2) after computing
;   the derivatives:
;
;       f'  = 2 Sum(hi  hi')     
;       f'' = 2 Sum(hi' hj') + 2 Sum(hi hi'')                (4)
;
;   If one assumes that the parameters are already close enough to a
;   minimum, then one typically finds that the second term in f'' is
;   negligible [or, in any case, is too difficult to compute].  Thus,
;   equation (2) can be solved, at least approximately, using only
;   gradient information.
;
;   In matrix notation, the combination of eqns (2) and (4) becomes:
;
;        hT' . h' . dx = - hT' . h                          (5)
;
;   Where h is the residual vector (length m), hT is its transpose, h'
;   is the Jacobian matrix (dimensions n x m), and dx is (xm-x0).  The
;   user function supplies the residual vector h, and in some cases h'
;   when it is not found by finite differences (see MPFIT_FDJAC2,
;   which finds h and hT').  Even if dx is not the best absolute step
;   to take, it does provide a good estimate of the best *direction*,
;   so often a line minimization will occur along the dx vector
;   direction.
;
;   The method of solution employed by MINPACK is to form the Q . R
;   factorization of h', where Q is an orthogonal matrix such that QT .
;   Q = I, and R is upper right triangular.  Using h' = Q . R and the
;   ortogonality of Q, eqn (5) becomes
;
;        (RT . QT) . (Q . R) . dx = - (RT . QT) . h
;                     RT . R . dx = - RT . QT . h         (6)
;                          R . dx = - QT . h
;
;   where the last statement follows because R is upper triangular.
;   Here, R, QT and h are known so this is a matter of solving for dx.
;   The routine MPFIT_QRFAC provides the QR factorization of h, with
;   pivoting, and MPFIT_QRSOL;V provides the solution for dx.
;   
; REFERENCES:
;
;   Markwardt, C. B. 2008, "Non-Linear Least Squares Fitting in IDL
;     with MPFIT," in proc. Astronomical Data Analysis Software and
;     Systems XVIII, Quebec, Canada, ASP Conference Series, Vol. XXX, eds.
;     D. Bohlender, P. Dowler & D. Durand (Astronomical Society of the
;     Pacific: San Francisco), p. 251-254 (ISBN: 978-1-58381-702-5)
;       http://arxiv.org/abs/0902.2850
;       Link to NASA ADS: http://adsabs.harvard.edu/abs/2009ASPC..411..251M
;       Link to ASP: http://aspbooks.org/a/volumes/table_of_contents/411
;
;   Refer to the MPFIT website as:
;       http://purl.com/net/mpfit
;
;   MINPACK-1 software, by Jorge More' et al, available from netlib.
;     http://www.netlib.org/
;
;   "Optimization Software Guide," Jorge More' and Stephen Wright, 
;     SIAM, *Frontiers in Applied Mathematics*, Number 14.
;     (ISBN: 978-0-898713-22-0)
;
;   More', J. 1978, "The Levenberg-Marquardt Algorithm: Implementation
;     and Theory," in Numerical Analysis, vol. 630, ed. G. A. Watson
;     (Springer-Verlag: Berlin), p. 105 (DOI: 10.1007/BFb0067690 )
;
; MODIFICATION HISTORY:
;   Translated from MINPACK-1 in FORTRAN, Apr-Jul 1998, CM
;   Fixed bug in parameter limits (x vs xnew), 04 Aug 1998, CM
;   Added PERROR keyword, 04 Aug 1998, CM
;   Added COVAR keyword, 20 Aug 1998, CM
;   Added NITER output keyword, 05 Oct 1998
;      D.L Windt, Bell Labs, windt@bell-labs.com;
;   Made each PARINFO component optional, 05 Oct 1998 CM
;   Analytical derivatives allowed via AUTODERIVATIVE keyword, 09 Nov 1998
;   Parameter values can be tied to others, 09 Nov 1998
;   Fixed small bugs (Wayne Landsman), 24 Nov 1998
;   Added better exception error reporting, 24 Nov 1998 CM
;   Cosmetic documentation changes, 02 Jan 1999 CM
;   Changed definition of ITERPROC to be consistent with TNMIN, 19 Jan 1999 CM
;   Fixed bug when AUTDERIVATIVE=0.  Incorrect sign, 02 Feb 1999 CM
;   Added keyboard stop to MPFIT_DEFITER, 28 Feb 1999 CM
;   Cosmetic documentation changes, 14 May 1999 CM
;   IDL optimizations for speed & FASTNORM keyword, 15 May 1999 CM
;   Tried a faster version of mpfit_enorm, 30 May 1999 CM
;   Changed web address to cow.physics.wisc.edu, 14 Jun 1999 CM
;   Found malformation of FDJAC in MPFIT for 1 parm, 03 Aug 1999 CM
;   Factored out user-function call into MPFIT_CALL.  It is possible,
;     but currently disabled, to call procedures.  The calling format
;     is similar to CURVEFIT, 25 Sep 1999, CM
;   Slightly changed mpfit_tie to be less intrusive, 25 Sep 1999, CM
;   Fixed some bugs associated with tied parameters in mpfit_fdjac, 25
;     Sep 1999, CM
;   Reordered documentation; now alphabetical, 02 Oct 1999, CM
;   Added QUERY keyword for more robust error detection in drivers, 29
;     Oct 1999, CM
;   Documented PERROR for unweighted fits, 03 Nov 1999, CM
;   Split out MPFIT_RESETPROF to aid in profiling, 03 Nov 1999, CM
;   Some profiling and speed optimization, 03 Nov 1999, CM
;     Worst offenders, in order: fdjac2, qrfac, qrsolv, enorm.
;     fdjac2 depends on user function, qrfac and enorm seem to be
;     fully optimized.  qrsolv probably could be tweaked a little, but
;     is still <10% of total compute time.
;   Made sure that !err was set to 0 in MPFIT_DEFITER, 10 Jan 2000, CM
;   Fixed small inconsistency in setting of QANYLIM, 28 Jan 2000, CM
;   Added PARINFO field RELSTEP, 28 Jan 2000, CM
;   Converted to MPFIT_ERROR common block for indicating error
;     conditions, 28 Jan 2000, CM
;   Corrected scope of MPFIT_ERROR common block, CM, 07 Mar 2000
;   Minor speed improvement in MPFIT_ENORM, CM 26 Mar 2000
;   Corrected case where ITERPROC changed parameter values and
;     parameter values were TIED, CM 26 Mar 2000
;   Changed MPFIT_CALL to modify NFEV automatically, and to support
;     user procedures more, CM 26 Mar 2000
;   Copying permission terms have been liberalized, 26 Mar 2000, CM
;   Catch zero value of zero a(j,lj) in MPFIT_QRFAC, 20 Jul 2000, CM
;      (thanks to David Schlegel <schlegel@astro.princeton.edu>)
;   MPFIT_SETMACHAR is called only once at init; only one common block
;     is created (MPFIT_MACHAR); it is now a structure; removed almost
;     all CHECK_MATH calls for compatibility with IDL5 and !EXCEPT;
;     profiling data is now in a structure too; noted some
;     mathematical discrepancies in Linux IDL5.0, 17 Nov 2000, CM
;   Some significant changes.  New PARINFO fields: MPSIDE, MPMINSTEP,
;     MPMAXSTEP.  Improved documentation.  Now PTIED constraints are
;     maintained in the MPCONFIG common block.  A new procedure to
;     parse PARINFO fields.  FDJAC2 now computes a larger variety of
;     one-sided and two-sided finite difference derivatives.  NFEV is
;     stored in the MPCONFIG common now.  17 Dec 2000, CM
;   Added check that PARINFO and XALL have same size, 29 Dec 2000 CM
;   Don't call function in TERMINATE when there is an error, 05 Jan
;     2000
;   Check for float vs. double discrepancies; corrected implementation
;     of MIN/MAXSTEP, which I still am not sure of, but now at least
;     the correct behavior occurs *without* it, CM 08 Jan 2001
;   Added SCALE_FCN keyword, to allow for scaling, as for the CASH
;     statistic; added documentation about the theory of operation,
;     and under the QR factorization; slowly I'm beginning to
;     understand the bowels of this algorithm, CM 10 Jan 2001
;   Remove MPMINSTEP field of PARINFO, for now at least, CM 11 Jan
;     2001
;   Added RESDAMP keyword, CM, 14 Jan 2001
;   Tried to improve the DAMP handling a little, CM, 13 Mar 2001
;   Corrected .PARNAME behavior in _DEFITER, CM, 19 Mar 2001
;   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
;   Documentation of how TIED interacts with LIMITS, CM, 21 Dec 2008
;   Better documentation of references, CM, 27 Feb 2009
;   If MAXITER=0, then be sure to set STATUS=5, which permits the
;      the covariance matrix to be computed, CM, 14 Apr 2009
;   Avoid numerical underflow while solving for the LM parameter,
;      (thanks to Sergey Koposov) CM, 14 Apr 2009
;   Use individual functions for all possible MPFIT_CALL permutations,
;      (and make sure the syntax is right) CM, 01 Sep 2009
;   Correct behavior of MPMAXSTEP when some parameters are frozen,
;      thanks to Josh Destree, CM, 22 Nov 2009
;   Update the references section, CM, 22 Nov 2009
;   1.70 - Add the VERSION and MIN_VERSION keywords, CM, 22 Nov 2009
;   1.71 - Store pre-calculated revision in common, CM, 23 Nov 2009
;   1.72-1.74 - Documented alternate method to compute correlation matrix,
;          CM, 05 Feb 2010
;   1.75 - Enforce TIED constraints when preparing to terminate the
;          routine, CM, 2010-06-22
;   1.76 - Documented input keywords now are not modified upon output,
;          CM, 2010-07-13
;   1.77 - Upon user request (/CALC_FJAC), compute Jacobian matrix and
;          return in BEST_FJAC; also return best residuals in
;          BEST_RESID; also return an index list of free parameters as
;          PFREE_INDEX; add a fencepost to prevent recursion
;          CM, 2010-10-27
;   1.79 - Documentation corrections.  CM, 2011-08-26
;   1.81 - Fix bug in interaction of AUTODERIVATIVE=0 and .MPSIDE=3;
;          Document FJAC_MASK. CM, 2012-05-08
;
;  $Id: mpfit.pro,v 1.82 2012/09/27 23:59:44 cmarkwar Exp $
;-
; Original MINPACK by More' Garbow and Hillstrom, translated with permission
; Modifications and enhancements are:
; Copyright (C) 1997-2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012 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 with no _EXTRA parameters
function mpfit_call_func_noextra, fcn, x, fjac, _EXTRA=extra
  if n_params() EQ 2 then begin
     return, call_function(fcn, x)
  endif else begin
     return, call_function(fcn, x, fjac)
  endelse
end

; Call user function with _EXTRA parameters
function mpfit_call_func_extra, fcn, x, fjac, _EXTRA=extra
  if n_params() EQ 2 then begin
     return, call_function(fcn, x, _EXTRA=extra)
  endif else begin
     return, call_function(fcn, x, fjac, _EXTRA=extra)
  endelse
end

; Call user procedure with no _EXTRA parameters
function mpfit_call_pro_noextra, fcn, x, fjac, _EXTRA=extra
  if n_params() EQ 2 then begin
     call_procedure, fcn, x, f
  endif else begin
     call_procedure, fcn, x, f, fjac
  endelse
  return, f
end

; Call user procedure with _EXTRA parameters
function mpfit_call_pro_extra, fcn, x, fjac, _EXTRA=extra
  if n_params() EQ 2 then begin
     call_procedure, fcn, x, f, _EXTRA=extra
  endif else begin
     call_procedure, fcn, x, f, fjac, _EXTRA=extra
  endelse
  return, f
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, and 
  ;; with/without FUNCTARGS
  proname = 'MPFIT_CALL'
  proname = proname + ((mpconfig.proc) ? '_PRO' : '_FUNC')
  proname = proname + ((n_elements(extra) GT 0) ? '_EXTRA' : '_NOEXTRA')

  if n_params() EQ 2 then begin
     f = call_function(proname, fcn, x, _EXTRA=extra)
  endif else begin
     f = call_function(proname, fcn, x, fjac, _EXTRA=extra)
  endelse
  mpconfig.nfev = mpconfig.nfev + 1

  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, $
                 deriv_debug=ddebug, deriv_reltol=ddrtol, deriv_abstol=ddatol

  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.
  if n_elements(ddebug) EQ 0 then ddebug = intarr(n_elements(xall))
  if n_elements(ddrtol) EQ 0 then ddrtol = x * 0.
  if n_elements(ddatol) EQ 0 then ddatol = x * 0.
  has_debug_deriv = max(ddebug)

  if keyword_set(has_debug_deriv) then begin
      ;; Header for debugging
      print, 'FJAC DEBUG BEGIN'
      print, "IPNT", "FUNC", "DERIV_U", "DERIV_N", "DIFF_ABS", "DIFF_REL", $
        format='("#  ",A10," ",A10," ",A10," ",A10," ",A10," ",A10)'
  endif

  nall = n_elements(xall)

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

  ;; Compute analytical derivative if requested
  ;; Two ways to enable computation of explicit derivatives:
  ;;   1. AUTODERIVATIVE=0
  ;;   2. AUTODERIVATIVE=1, but P[i].MPSIDE EQ 3

  if keyword_set(autoderiv) EQ 0 OR max(dside[ifree] EQ 3) EQ 1 then begin
      fjac_mask = intarr(nall)
      ;; Specify which parameters need derivatives
      ;;                  ---- Case 2 ------     ----- Case 1 -----
      fjac_mask[ifree] = (dside[ifree] EQ 3) OR (keyword_set(autoderiv) EQ 0)
      if has_debug_deriv then $
        print, fjac_mask, format='("# FJAC_MASK = ",100000(I0," ",:))'

      fjac = fjac_mask  ;; Pass the mask to the calling function as FJAC
      mperr = 0
      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 (WRONG, see below)
      ;; Sign error found (thanks Jesus Fernandez <fernande@irm.chu-caen.fr>)

      ;; ... 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)
      
      ;; Are we done computing derivatives?  The answer is, YES, if we
      ;; computed explicit derivatives for all free parameters, EXCEPT
      ;; when we are going on to compute debugging derivatives.
      if min(fjac_mask[ifree]) EQ 1 AND NOT has_debug_deriv then begin
          return, fjac
      endif
  endif

  ;; Final output array, if it was not already created above
  if n_elements(fjac) EQ 0 then begin
      fjac = make_array(m, n, value=fvec[0]*0.)
      fjac = reform(fjac, m, n, /overwrite)
  endif

  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
      dsidej = dside[ifree[j]]
      ddebugj = ddebug[ifree[j]]

      ;; Skip this parameter if we already computed its derivative
      ;; explicitly, and we are not debugging.
      if (dsidej EQ 3) AND (ddebugj EQ 0) then continue
      if (dsidej EQ 3) AND (ddebugj EQ 1) then $
        print, ifree[j], format='("FJAC PARM ",I0)'

      xp = xall
      xp[ifree[j]] = xp[ifree[j]] + h[j]
      
      mperr = 0
      fp = mpfit_call(fcn, xp, _EXTRA=fcnargs)
      
      iflag = mperr
      if iflag LT 0 then return, !values.d_nan

      if ((dsidej GE -1) AND (dsidej LE 1)) OR (dsidej EQ 3) then begin
          ;; COMPUTE THE ONE-SIDED DERIVATIVE
          ;; Note optimization fjac(0:*,j)
          fjacj = (fp-fvec)/h[j]

      endif else begin
          ;; COMPUTE THE TWO-SIDED DERIVATIVE
          xp[ifree[j]] = xall[ifree[j]] - h[j]

          mperr = 0
          fm = mpfit_call(fcn, xp, _EXTRA=fcnargs)
          
          iflag = mperr
          if iflag LT 0 then return, !values.d_nan
          
          ;; Note optimization fjac(0:*,j)
          fjacj = (fp-fm)/(2*h[j])
      endelse          
      
      ;; Debugging of explicit derivatives
      if (dsidej EQ 3) AND (ddebugj EQ 1) then begin
          ;; Relative and absolute tolerances
          dr = ddrtol[ifree[j]] & da = ddatol[ifree[j]]

          ;; Explicitly calculated
          fjaco = fjac[*,j]
          
          ;; If tolerances are zero, then any value for deriv triggers print...
          if (da EQ 0 AND dr EQ 0) then $
            diffj = (fjaco NE 0 OR fjacj NE 0)
          ;; ... otherwise the difference must be a greater than tolerance
          if (da NE 0 OR dr NE 0) then $
            diffj = (abs(fjaco-fjacj) GT (da+abs(fjaco)*dr))

          for k = 0L, m-1 do if diffj[k] then begin
              print, k, fvec[k], fjaco[k], fjacj[k], fjaco[k]-fjacj[k], $
                (fjaco[k] EQ 0)?(0):((fjaco[k]-fjacj[k])/fjaco[k]), $
                format='("   ",I10," ",G10.4," ",G10.4," ",G10.4," ",G10.4," ",G10.4)'
          endif
      endif

      ;; Store final results in output array
      fjac[0,j] = fjacj
          
  endfor

  if has_debug_deriv then print, 'FJAC DEBUG END'

;  profvals.fdjac2 = profvals.fdjac2 + (systime(1) - prof_start)
  return, fjac
end

function mpfit_enorm, vec

  COMPILE_OPT strictarr
  ;; NOTE: it turns out that, for systems that have a lot of data
  ;; points, this routine is a big computing bottleneck.  The extended
  ;; computations that need to be done cannot be effectively
  ;; vectorized.  The introduction of the FASTNORM configuration
  ;; parameter allows the user to select a faster routine, which is 
  ;; based on TOTAL() alone.
  common mpfit_profile, profvals
;  prof_start = systime(1)

  common mpfit_config, mpconfig
; Very simple-minded sum-of-squares
  if n_elements(mpconfig) GT 0 then if mpconfig.fastnorm then begin
      ans = sqrt(total(vec^2))
      goto, TERMINATE
  endif

  common mpfit_machar, machvals

  agiant = machvals.rgiant / n_elements(vec)
  adwarf = machvals.rdwarf * n_elements(vec)

  ;; This is hopefully a compromise between speed and robustness.
  ;; Need to do this because of the possibility of over- or underflow.
  mx = max(vec, min=mn)
  mx = max(abs([mx,mn]))
  if mx EQ 0 then return, vec[0]*0.

  if mx GT agiant OR mx LT adwarf then ans = mx * sqrt(total((vec/mx)^2))$
  else                                 ans = sqrt( total(vec^2) )

  TERMINATE:
;  profvals.enorm = profvals.enorm + (systime(1) - prof_start)
  return, ans
end

;     **********
;
;     subroutine qrfac
;
;     this subroutine uses householder transformations with column
;     pivoting (optional) to compute a qr factorization of the
;     m by n matrix a. that is, qrfac determines an orthogonal
;     matrix q, a permutation matrix p, and an upper trapezoidal
;     matrix r with diagonal elements of nonincreasing magnitude,
;     such that a*p = q*r. the householder transformation for
;     column k, k = 1,2,...,min(m,n), is of the form
;
;         t
;     i - (1/u(k))*u*u
;
;     where u has zeros in the first k-1 positions. the form of
;     this transformation and the method of pivoting first
;     appeared in the corresponding linpack subroutine.
;
;     the subroutine statement is
;
; subroutine qrfac(m,n,a,lda,pivot,ipvt,lipvt,rdiag,acnorm,wa)
;
;     where
;
; m is a positive integer input variable set to the number
;   of rows of a.
;
; n is a positive integer input variable set to the number
;   of columns of a.
;
; a is an m by n array. on input a contains the matrix for
;   which the qr factorization is to be computed. on output
;   the strict upper trapezoidal part of a contains the strict
;   upper trapezoidal part of r, and the lower trapezoidal
;   part of a contains a factored form of q (the non-trivial
;   elements of the u vectors described above).
;
; lda is a positive integer input variable not less than m
;   which specifies the leading dimension of the array a.
;
; pivot is a logical input variable. if pivot is set true,
;   then column pivoting is enforced. if pivot is set false,
;   then no column pivoting is done.
;
; ipvt is an integer output array of length lipvt. ipvt
;   defines the permutation matrix p such that a*p = q*r.
;   column j of p is column ipvt(j) of the identity matrix.
;   if pivot is false, ipvt is not referenced.
;
; lipvt is a positive integer input variable. if pivot is false,
;   then lipvt may be as small as 1. if pivot is true, then
;   lipvt must be at least n.
;
; rdiag is an output array of length n which contains the
;   diagonal elements of r.
;
; acnorm is an output array of length n which contains the
;   norms of the corresponding columns of the input matrix a.
;   if this information is not needed, then acnorm can coincide
;   with rdiag.
;
; wa is a work array of length n. if pivot is false, then wa
;   can coincide with rdiag.
;
;     subprograms called
;
; minpack-supplied ... dpmpar,enorm
;
; fortran-supplied ... dmax1,dsqrt,min0
;
;     argonne national laboratory. minpack project. march 1980.
;     burton s. garbow, kenneth e. hillstrom, jorge j. more
;
;     **********
;
; PIVOTING / PERMUTING:
;
; Upon return, A(*,*) is in standard parameter order, A(*,IPVT) is in
; permuted order.
;
; RDIAG is in permuted order.
;
; ACNORM is in standard parameter order.
;
; NOTE: in IDL the factors appear slightly differently than described
; above.  The matrix A is still m x n where m >= n.  
;
; The "upper" triangular matrix R is actually stored in the strict
; lower left triangle of A under the standard notation of IDL.
;
; The reflectors that generate Q are in the upper trapezoid of A upon
; output.
;
;  EXAMPLE:  decompose the matrix [[9.,2.,6.],[4.,8.,7.]]
;    aa = [[9.,2.,6.],[4.,8.,7.]]
;    mpfit_qrfac, aa, aapvt, rdiag, aanorm
;     IDL> print, aa
;          1.81818*     0.181818*     0.545455*
;         -8.54545+      1.90160*     0.432573*
;     IDL> print, rdiag
;         -11.0000+     -7.48166+
;
; The components marked with a * are the components of the
; reflectors, and those marked with a + are components of R.
;
; To reconstruct Q and R we proceed as follows.  First R.
;    r = fltarr(m, n)
;    for i = 0, n-1 do r(0:i,i) = aa(0:i,i)  ; fill in lower diag
;    r(lindgen(n)*(m+1)) = rdiag
;
; Next, Q, which are composed from the reflectors.  Each reflector v
; is taken from the upper trapezoid of aa, and converted to a matrix
; via (I - 2 vT . v / (v . vT)).
;
;   hh = ident                                    ;; identity matrix
;   for i = 0, n-1 do begin
;    v = aa(*,i) & if i GT 0 then v(0:i-1) = 0    ;; extract reflector
;    hh = hh ## (ident - 2*(v # v)/total(v * v))  ;; generate matrix
;   endfor
;
; Test the result:
;    IDL> print, hh ## transpose(r)
;          9.00000      4.00000
;          2.00000      8.00000
;          6.00000      7.00000
;
; Note that it is usually never necessary to form the Q matrix
; explicitly, and MPFIT does not.

pro mpfit_qrfac, a, ipvt, rdiag, acnorm, pivot=pivot

  COMPILE_OPT strictarr
  sz = size(a)
  m = sz[1]
  n = sz[2]

  common mpfit_machar, machvals
  common mpfit_profile, profvals
;  prof_start = systime(1)

  MACHEP0 = machvals.machep
  DWARF   = machvals.minnum
  
  ;; Compute the initial column norms and initialize arrays
  acnorm = make_array(n, value=a[0]*0.)
  for j = 0L, n-1 do $
    acnorm[j] = mpfit_enorm(a[*,j])
  rdiag = acnorm
  wa = rdiag
  ipvt = lindgen(n)

  ;; Reduce a to r with householder transformations
  minmn = min([m,n])
  for j = 0L, minmn-1 do begin
      if NOT keyword_set(pivot) then goto, HOUSE1
      
      ;; Bring the column of largest norm into the pivot position
      rmax = max(rdiag[j:*])
      kmax = where(rdiag[j:*] EQ rmax, ct) + j
      if ct LE 0 then goto, HOUSE1
      kmax = kmax[0]
      
      ;; Exchange rows via the pivot only.  Avoid actually exchanging
      ;; the rows, in case there is lots of memory transfer.  The
      ;; exchange occurs later, within the body of MPFIT, after the
      ;; extraneous columns of the matrix have been shed.
      if kmax NE j then begin
          temp     = ipvt[j]   & ipvt[j]    = ipvt[kmax] & ipvt[kmax]  = temp
          rdiag[kmax] = rdiag[j]
          wa[kmax]    = wa[j]
      endif
      
      HOUSE1:

      ;; Compute the householder transformation to reduce the jth
      ;; column of A to a multiple of the jth unit vector
      lj     = ipvt[j]
      ajj    = a[j:*,lj]
      ajnorm = mpfit_enorm(ajj)
      if ajnorm EQ 0 then goto, NEXT_ROW
      if a[j,lj] LT 0 then ajnorm = -ajnorm
      
      ajj     = ajj / ajnorm
      ajj[0]  = ajj[0] + 1
      ;; *** Note optimization a(j:*,j)
      a[j,lj] = ajj
      
      ;; Apply the transformation to the remaining columns
      ;; and update the norms

      ;; NOTE to SELF: tried to optimize this by removing the loop,
      ;; but it actually got slower.  Reverted to "for" loop to keep
      ;; it simple.
      if j+1 LT n then begin
          for k=j+1, n-1 do begin
              lk = ipvt[k]
              ajk = a[j:*,lk]
              ;; *** Note optimization a(j:*,lk) 
              ;; (corrected 20 Jul 2000)
              if a[j,lj] NE 0 then $
                a[j,lk] = ajk - ajj * total(ajk*ajj)/a[j,lj]

              if keyword_set(pivot) AND rdiag[k] NE 0 then begin
                  temp = a[j,lk]/rdiag[k]
                  rdiag[k] = rdiag[k] * sqrt((1.-temp^2) > 0)
                  temp = rdiag[k]/wa[k]
                  if 0.05D*temp*temp LE MACHEP0 then begin
                      rdiag[k] = mpfit_enorm(a[j+1:*,lk])
                      wa[k] = rdiag[k]
                  endif
              endif
          endfor
      endif

      NEXT_ROW:
      rdiag[j] = -ajnorm
  endfor

;  profvals.qrfac = profvals.qrfac + (systime(1) - prof_start)
  return
end

;     **********
;
;     subroutine qrsolv
;
;     given an m by n matrix a, an n by n diagonal matrix d,
;     and an m-vector b, the problem is to determine an x which
;     solves the system
;
;           a*x = b ,     d*x = 0 ,
;
;     in the least squares sense.
;
;     this subroutine completes the solution of the problem
;     if it is provided with the necessary information from the
;     qr factorization, with column pivoting, of a. that is, if
;     a*p = q*r, where p is a permutation matrix, q has orthogonal
;     columns, and r is an upper triangular matrix with diagonal
;     elements of nonincreasing magnitude, then qrsolv expects
;     the full upper triangle of r, the permutation matrix p,
;     and the first n components of (q transpose)*b. the system
;     a*x = b, d*x = 0, is then equivalent to
;
;                  t       t
;           r*z = q *b ,  p *d*p*z = 0 ,
;
;     where x = p*z. if this system does not have full rank,
;     then a least squares solution is obtained. on output qrsolv
;     also provides an upper triangular matrix s such that
;
;            t   t               t
;           p *(a *a + d*d)*p = s *s .
;
;     s is computed within qrsolv and may be of separate interest.
;
;     the subroutine statement is
;
;       subroutine qrsolv(n,r,ldr,ipvt,diag,qtb,x,sdiag,wa)
;
;     where
;
;       n is a positive integer input variable set to the order of r.
;
;       r is an n by n array. on input the full upper triangle
;         must contain the full upper triangle of the matrix r.
;         on output the full upper triangle is unaltered, and the
;         strict lower triangle contains the strict upper triangle
;         (transposed) of the upper triangular matrix s.
;
;       ldr is a positive integer input variable not less than n
;         which specifies the leading dimension of the array r.
;
;       ipvt is an integer input array of length n which defines the
;         permutation matrix p such that a*p = q*r. column j of p
;         is column ipvt(j) of the identity matrix.
;
;       diag is an input array of length n which must contain the
;         diagonal elements of the matrix d.
;
;       qtb is an input array of length n which must contain the first
;         n elements of the vector (q transpose)*b.
;
;       x is an output array of length n which contains the least
;         squares solution of the system a*x = b, d*x = 0.
;
;       sdiag is an output array of length n which contains the
;         diagonal elements of the upper triangular matrix s.
;
;       wa is a work array of length n.
;
;     subprograms called
;
;       fortran-supplied ... dabs,dsqrt
;
;     argonne national laboratory. minpack project. march 1980.
;     burton s. garbow, kenneth e. hillstrom, jorge j. more
;
pro mpfit_qrsolv, r, ipvt, diag, qtb, x, sdiag

  COMPILE_OPT strictarr
  sz = size(r)
  m = sz[1]
  n = sz[2]
  delm = lindgen(n) * (m+1) ;; Diagonal elements of r

  common mpfit_profile, profvals
;  prof_start = systime(1)

  ;; copy r and (q transpose)*b to preserve input and initialize s.
  ;; in particular, save the diagonal elements of r in x.

  for j = 0L, n-1 do $