;+
; NAME:
;   DUSTEM_MPFITFUN
;
; ORIGINAL AUTHOR:
;   Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
;   craigm@lheamail.gsfc.nasa.gov
;   UPDATED VERSIONs can be found on my WEB PAGE: 
;      http://cow.physics.wisc.edu/~craigm/idl/idl.html
;
; PURPOSE:
;   *** Note that this is a direct copy of CM's MPFITFUN
;   code. Only the names of subroutines have been modified according
;   to the dustem_*** convention of the DustEMWrap code. Renaming has
;   been done to facilitate debugging and user support in case of
;   issues, and to avoid conflicts with any existing installation of
;   MPFIT IDL library.***
;
;   Perform Levenberg-Marquardt least-squares fit to IDL function
;
; MAJOR TOPICS:
;   Curve and Surface Fitting
;
; CALLING SEQUENCE:
;   parms = DUSTEM_MPFITFUN(MYFUNCT, X, Y, ERR, start_params, ...)
;
; DESCRIPTION:
;
;  MPFITFUN fits a user-supplied model -- in the form of an IDL
;  function -- to a set of user-supplied data.  MPFITFUN calls
;  MPFIT, the MINPACK-1 least-squares minimizer, to do the main
;  work.
;
;  Given the data and their uncertainties, MPFITFUN finds the best set
;  of model parameters which match the data (in a least-squares
;  sense) and returns them in an array.
;  
;  The user must supply the following items:
;   - An array of independent variable values ("X").
;   - An array of "measured" *dependent* variable values ("Y").
;   - An array of "measured" 1-sigma uncertainty values ("ERR").
;   - The name of an IDL function which computes Y given X ("MYFUNCT").
;   - Starting guesses for all of the parameters ("START_PARAMS").
;
;  There are very few restrictions placed on X, Y or MYFUNCT.  Simply
;  put, MYFUNCT must map the "X" values into "Y" values given the
;  model parameters.  The "X" values may represent any independent
;  variable (not just Cartesian X), and indeed may be multidimensional
;  themselves.  For example, in the application of image fitting, X
;  may be a 2xN array of image positions.
;
;  Data values of NaN or Infinity for "Y", "ERR" or "WEIGHTS" will be
;  ignored as missing data if the NAN keyword is set.  Otherwise, they
;  may cause the fitting loop to halt with an error message.  Note
;  that the fit will still halt if the model function, or its
;  derivatives, produces infinite or NaN values.
;
;  MPFITFUN carefully avoids passing large arrays where possible to
;  improve performance.
;
;  See below for an example of usage.
;
; USER FUNCTION
;
;  The user must define a function which returns the model value.  For
;  applications which use finite-difference derivatives -- the default
;  -- the user function should be declared in the following way:
;
;    FUNCTION MYFUNCT, X, P
;     ; The independent variable is X
;     ; Parameter values are passed in "P"
;     YMOD = ... computed model values at X ...
;     return, YMOD
;    END
;
;  The returned array YMOD must have the same dimensions and type as
;  the "measured" Y values.
;
;  User functions may also indicate a fatal error condition
;  using the ERROR_CODE common block variable, as described
;  below under the MPFIT_ERROR common block definition.
;
;  MPFIT by default calculates derivatives numerically via a finite
;  difference approximation.  However, the user function *may*
;  calculate the derivatives if desired, but only if the model
;  function is declared with an additional position parameter, DP, as
;  described below.
;
;  To enable explicit derivatives for all parameters, set
;  AUTODERIVATIVE=0.
;
;  When AUTODERIVATIVE=0, the user function is responsible for
;  calculating the derivatives of the user function with respect to
;  each parameter.  The user function should be declared as follows:
;
;    ;
;    ; MYFUNCT - example user function
;    ;   P - input parameter values (N-element array)
;    ;   DP - upon input, an N-vector indicating which parameters
;    ;          to compute derivatives for; 
;    ;        upon output, the user function must return
;    ;          an ARRAY(M,N) of derivatives in this keyword
;    ;   (keywords) - any other keywords specified by FUNCTARGS
;    ; RETURNS - function values
;    ;
;    FUNCTION MYFUNCT, x, p, dp [, (additional keywords if desired)]
;     model = F(x, p)         ;; Model function
;     
;     if n_params() GT 2 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)
;     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.
;
;  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
;  function value with respect to the jth parameter.
;
;  MPFIT may not always request derivatives from the user function.
;  In those cases, the parameter DP is not passed.  Therefore
;  functions can use N_PARAMS() to indicate whether they must compute
;  the derivatives or not.
;
;  For additional information about explicit derivatives, including
;  additional settings and debugging options, see the discussion under
;  "EXPLICIT DERIVATIVES" and AUTODERIVATIVE in MPFIT.PRO.
;
; 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.
;  
;  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 rely on the EXECUTE() function, they cannot
;  be used with the free version of the IDL Virtual Machine.
;
;
; INPUTS:
;   MYFUNCT - a string variable containing the name of an IDL function.
;             This function computes the "model" Y values given the
;             X values and model parameters, as desribed above.
;
;   X - Array of independent variable values.
;
;   Y - Array of "measured" dependent variable values.  Y should have
;       the same data type as X.  The function MYFUNCT should map
;       X->Y.
;         NOTE: the following special cases apply:
;                * if Y is NaN or Infinite, and the NAN keyword is
;                  set, then the corresponding data point is ignored
;
;   ERR - Array of "measured" 1-sigma uncertainties.  ERR should have
;         the same data type as Y.  ERR is ignored if the WEIGHTS
;         keyword is specified.
;         NOTE: the following special cases apply:
;                * if ERR is zero, then the corresponding data point
;                  is ignored
;                * if ERR is NaN or Infinite, and the NAN keyword is
;                  set, then the corresponding data point is ignored
;                * if ERR is negative, then the absolute value of 
;                  ERR is used.
;
;   START_PARAMS - An 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 (see MPFIT).  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:
;
;   BESTNORM - the value of the summed squared residuals for the
;              returned parameter values.
;
;   BEST_FJAC - upon return, BEST_FJAC contains the Jacobian, or
;               partial derivative, matrix for the best-fit model.
;               The values are an array,
;               ARRAY(N_ELEMENTS(DEVIATES),NFREE) where NFREE is the
;               number of free parameters.  This array is only
;               computed if /CALC_FJAC is set, otherwise BEST_FJAC is
;               undefined.
;
;               The returned array is such that BEST_FJAC[I,J] is the
;               partial derivative of the model with respect to
;               parameter PARMS[PFREE_INDEX[J]].
;
;   BEST_RESID - upon return, an array of best-fit deviates,
;                normalized by the weights or errors.
;
;   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:
;                  PCOR = COV * 0
;                  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])
;           or equivalently, in vector notation,
;                  PCOR = COV / (PERROR # PERROR)
;
;           If NOCOVAR is set or MPFIT terminated abnormally, then
;           COVAR is set to a scalar with value !VALUES.D_NAN.
;
;   CASH - when set, the fit statistic is changed to a derivative of
;          the CASH statistic.  The model function must be strictly
;          positive. WARNING: this option is incomplete and untested.
;
;   DOF - number of degrees of freedom, computed as
;             DOF = N_ELEMENTS(DEVIATES) - NFREE
;         Note that this doesn't account for pegged parameters (see
;         NPEGGED).  It also does not account for data points which
;         are assigned zero weight, for example if :
;           * WEIGHTS[i] EQ 0, or
;           * ERR[i] EQ infinity, or 
;           * any of the values is "undefined" and /NAN is set.
;
;   ERRMSG - a string error or warning message is returned.
;
;   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.
;
;               By default, no extra parameters are passed to the
;               user-supplied function.
;
;   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.
;
;   ITERPROC - The name of a procedure to be called upon each NPRINT
;              iteration of the MPFIT routine.  It should be declared
;              in the following way:
;
;              PRO ITERPROC, MYFUNCT, p, iter, fnorm, FUNCTARGS=fcnargs, $
;                PARINFO=parinfo, QUIET=quiet, ...
;                ; 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.  See below for documentation of
;              PARINFO.
;
;              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 (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.
;
;   MAXITER - The maximum number of iterations to perform.  If the
;             number of calculation iterations exceeds MAXITER, then
;             the STATUS value is set to 5 and MPFIT returns.  
;
;             If MAXITER EQ 0, then MPFIT does not iterate to adjust
;             parameter values; however, the user function is evaluated
;             and parameter errors/covariance/Jacobian are estimated
;             before returning.
;             Default: 200 iterations
;
;   NAN - ignore infinite or NaN values in the Y, ERR or WEIGHTS
;         parameters.  These values will be treated as missing data.
;         However, the fit will still halt with an error condition
;         if the model function becomes infinite.
;
;   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.
;
;   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.
;
;             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
;
;   PFREE_INDEX - upon return, PFREE_INDEX contains an index array
;                 which indicates which parameter were allowed to
;                 vary.  I.e. of all the parameters PARMS, only
;                 PARMS[PFREE_INDEX] were varied.
;
;   QUERY - if set, then MPFIT() will return immediately with one of
;           the following values:
;                 1 - if MIN_VERSION is not set
;                 1 - if MIN_VERSION is set and MPFIT satisfies the minimum
;                 0 - if MIN_VERSION is set and MPFIT does not satisfy it
;           Default: not set.
;
;   QUIET - set this keyword when no textual output should be printed
;           by MPFIT
;
;   STATUS - an integer status code is returned.  All values greater
;            than zero can represent success (however STATUS EQ 5 may
;            indicate failure to converge).  It can have one of the
;            following values:
;
;        -18  a fatal execution error has occurred.  More information
;             may be available in the ERRMSG string.
;
;        -16  a parameter or function value has become infinite or an
;             undefined number.  This is usually a consequence of
;             numerical overflow in the user's model function, which
;             must be avoided.
;
;        -15 to -1 
;             these are error codes that either MYFUNCT or ITERPROC
;             may return to terminate the fitting process (see
;             description of MPFIT_ERROR common below).  If either
;             MYFUNCT or ITERPROC set ERROR_CODE to a negative number,
;             then that number is returned in STATUS.  Values from -15
;             to -1 are reserved for the user functions and will not
;             clash with MPFIT.
;
;	   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 for standard weightings:
;
;                1D/ERR^2 - Normal weighting (ERR is the measurement error)
;                1D/Y     - Poisson weighting (counting statistics)
;                1D       - Unweighted
;
;         NOTE: the following special cases apply:
;                * if WEIGHTS is zero, then the corresponding data point
;                  is ignored
;                * if WEIGHTS is NaN or Infinite, and the NAN keyword is
;                  set, then the corresponding data point is ignored
;                * if WEIGHTS is negative, then the absolute value of 
;                  WEIGHTS is used.
;
;   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, windt@bell-labs.com;
;   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
;   Add documentation about NAN behavior, CM, 30 Mar 2009
;   Add keywords BEST_RESIDS, CALC_FJAC, BEST_FJAC, PFREE_INDEX;
;     update some documentation that had become stale, CM, 2010-10-28
;   Documentation corrections, CM, 2011-08-26
;   Additional documentation about explicit derivatives, CM, 2012-07-23
;
;  $Id: mpfitfun.pro,v 1.19 2012/09/27 23:59:31 cmarkwar Exp $
;-
; Copyright (C) 1997-2002, 2003, 2006, 2007, 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.
;-

FORWARD_FUNCTION dustem_mpfitfun_eval, dustem_mpfitfun, dustem_mpfit

; This is the call-back function for MPFIT.  It evaluates the
; function, subtracts the data, and returns the residuals.
function dustem_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 = 0L, 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 = 0L, 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 dustem_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 dustem_mpfitfun, fcn, x, y, err, p, WEIGHTS=wts, FUNCTARGS=fa, $
                   BESTNORM=bestnorm, nfev=nfev, STATUS=status, $
                   best_resid=best_resid, pfree_index=ifree, $
                   calc_fjac=calc_fjac, best_fjac=best_fjac, $
                   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
     DUSTEM_MPFIT_NOTFOUND:
     catch, /cancel
     message, 'ERROR: the required function DUSTEM_MPFIT must be in your IDL path', /info
     return, !values.d_nan
  endif
  if dustem_mpfit(/query) NE 1 then goto, DUSTEM_MPFIT_NOTFOUND
  
  catch, /cancel
  if keyword_set(query) then return, 1

  if n_params() EQ 0 then begin
      message, "USAGE: PARMS = DUSTEM_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 = 'dustem_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/abs(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 = dustem_mpfit('dustem_mpfitfun_eval', p, SCALE_FCN=scalfcn, $
                 parinfo=parinfo, STATUS=status, nfev=nfev, BESTNORM=bestnorm,$
                 covar=covar, perror=perror, $
                 best_resid=best_resid, pfree_index=ifree, $
                 calc_fjac=calc_fjac, best_fjac=best_fjac, $
                 niter=niter, nfree=nfree, npegged=npegged, dof=dof, $
                 ERRMSG=errmsg, quiet=quiet, _EXTRA=extra)

  ;; Retrieve the fit value
  yfit = temporary(mc)

  ;; Rescale the Jacobian according to parameter uncertainties
  if keyword_set(calc_fjac) AND nfree GT 0 AND status GT 0 then begin
      ec = 1/wc  ;; Per-data-point errors (could be INF or NAN!)
      for i = 0, nfree-1 do best_fjac[*,i] = - best_fjac[*,i] * ec
  endif

  ;; 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