;+ ; NAME: ; DUSTEM_MPFIT ; ; 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 MPFIT ; 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 minimization (MINPACK-1) ; ; MAJOR TOPICS: ; Curve and Surface Fitting ; ; CALLING SEQUENCE: ; parms = MPFIT(MYFUNCT, start_parms, FUNCTARGS=fcnargs, NFEV=nfev, ; MAXITER=maxiter, ERRMSG=errmsg, NPRINT=nprint, QUIET=quiet, ; FTOL=ftol, XTOL=xtol, GTOL=gtol, NITER=niter, ; STATUS=status, ITERPROC=iterproc, ITERARGS=iterargs, ; COVAR=covar, PERROR=perror, BESTNORM=bestnorm, ; PARINFO=parinfo) ; ; DESCRIPTION: ; ; MPFIT uses the Levenberg-Marquardt technique to solve the ; least-squares problem. In its typical use, MPFIT will be used to ; fit a user-supplied function (the "model") to user-supplied data ; points (the "data") by adjusting a set of parameters. MPFIT is ; based upon MINPACK-1 (LMDIF.F) by More' and collaborators. ; ; For example, a researcher may think that a set of observed data ; points is best modelled with a Gaussian curve. A Gaussian curve is ; parameterized by its mean, standard deviation and normalization. ; MPFIT will, within certain constraints, find the set of parameters ; which best fits the data. The fit is "best" in the least-squares ; sense; that is, the sum of the weighted squared differences between ; the model and data is minimized. ; ; The Levenberg-Marquardt technique is a particular strategy for ; iteratively searching for the best fit. This particular ; implementation is drawn from MINPACK-1 (see NETLIB), and seems to ; be more robust than routines provided with IDL. This version ; allows upper and lower bounding constraints to be placed on each ; parameter, or the parameter can be held fixed. ; ; The IDL user-supplied function should return an array of weighted ; deviations between model and data. In a typical scientific problem ; the residuals should be weighted so that each deviate has a ; gaussian sigma of 1.0. If X represents values of the independent ; variable, Y represents a measurement for each value of X, and ERR ; represents the error in the measurements, then the deviates could ; be calculated as follows: ; ; DEVIATES = (Y - F(X)) / ERR ; ; where F is the function representing the model. You are ; recommended to use the convenience functions MPFITFUN and ; MPFITEXPR, which are driver functions that calculate the deviates ; for you. If ERR are the 1-sigma uncertainties in Y, then ; ; TOTAL( DEVIATES^2 ) ; ; will be the total chi-squared value. MPFIT will minimize the ; chi-square value. The values of X, Y and ERR are passed through ; MPFIT to the user-supplied function via the FUNCTARGS keyword. ; ; Simple constraints can be placed on parameter values by using the ; PARINFO keyword to MPFIT. See below for a description of this ; keyword. ; ; MPFIT does not perform more general optimization tasks. See TNMIN ; instead. MPFIT is customized, based on MINPACK-1, to the ; least-squares minimization problem. ; ; USER FUNCTION ; ; The user must define a function which returns the appropriate ; values as specified above. The function should return the weighted ; deviations between the model and the data. For applications which ; use finite-difference derivatives -- the default -- the user ; function should be declared in the following way: ; ; FUNCTION MYFUNCT, p, X=x, Y=y, ERR=err ; ; Parameter values are passed in "p" ; model = F(x, p) ; return, (y-model)/err ; END ; ; See below for applications with explicit derivatives. ; ; The keyword parameters X, Y, and ERR in the example above are ; suggestive but not required. Any parameters can be passed to ; MYFUNCT by using the FUNCTARGS keyword to MPFIT. Use MPFITFUN and ; MPFITEXPR if you need ideas on how to do that. The function *must* ; accept a parameter list, P. ; ; In general there are no restrictions on the number of dimensions in ; X, Y or ERR. However the deviates *must* be returned in a ; one-dimensional array, and must have the same type (float or ; double) as the input arrays. ; ; See below for error reporting mechanisms. ; ; ; CHECKING STATUS AND HANNDLING ERRORS ; ; Upon return, MPFIT will report the status of the fitting operation ; in the STATUS and ERRMSG keywords. The STATUS keyword will contain ; a numerical code which indicates the success or failure status. ; Generally speaking, any value 1 or greater indicates success, while ; a value of 0 or less indicates a possible failure. The ERRMSG ; keyword will contain a text string which should describe the error ; condition more fully. ; ; By default, MPFIT will trap fatal errors and report them to the ; caller gracefully. However, during the debugging process, it is ; often useful to halt execution where the error occurred. When you ; set the NOCATCH keyword, MPFIT will not do any special error ; trapping, and execution will stop whereever the error occurred. ; ; MPFIT does not explicitly change the !ERROR_STATE variable ; (although it may be changed implicitly if MPFIT calls MESSAGE). It ; is the caller's responsibility to call MESSAGE, /RESET to ensure ; that the error state is initialized before calling MPFIT. ; ; User functions may also indicate non-fatal error conditions using ; the ERROR_CODE common block variable, as described below under the ; MPFIT_ERROR common block definition (by setting ERROR_CODE to a ; number between -15 and -1). When the user function sets an error ; condition via ERROR_CODE, MPFIT will gracefully exit immediately ; and report this condition to the caller. The ERROR_CODE is ; returned in the STATUS keyword in that case. ; ; ; EXPLICIT DERIVATIVES ; ; In the search for the best-fit solution, MPFIT by default ; calculates derivatives numerically via a finite difference ; approximation. The user-supplied function need not calculate the ; derivatives explicitly. However, 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. If the user function does not accept this additional ; parameter, MPFIT will report an error. As a practical matter, it ; is often sufficient and even faster to allow MPFIT to calculate the ; derivatives numerically, but this option is available for users who ; wish more control over the fitting process. ; ; There are two ways to enable explicit derivatives. First, the user ; can set the keyword AUTODERIVATIVE=0, which is a global switch for ; all parameters. In this case, MPFIT will request explicit ; derivatives for every free parameter. ; ; Second, the user may request explicit derivatives for specifically ; selected parameters using the PARINFO.MPSIDE=3 (see "CONSTRAINING ; PARAMETER VALUES WITH THE PARINFO KEYWORD" below). In this ; strategy, the user picks and chooses which parameter derivatives ; are computed explicitly versus numerically. When PARINFO[i].MPSIDE ; EQ 3, then the ith parameter derivative is computed explicitly. ; ; The keyword setting AUTODERIVATIVE=0 always globally overrides the ; individual values of PARINFO.MPSIDE. Setting AUTODERIVATIVE=0 is ; equivalent to resetting PARINFO.MPSIDE=3 for all parameters. ; ; Even if the user requests explicit derivatives for some or all ; parameters, MPFIT will not always request explicit derivatives on ; every user function call. ; ; EXPLICIT DERIVATIVES - CALLING INTERFACE ; ; When AUTODERIVATIVE=0, the user function is responsible for ; calculating the derivatives of the *residuals* with respect to each ; parameter. The user function should be declared as follows: ; ; ; ; ; 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 - residual values ; ; ; FUNCTION MYFUNCT, p, dp, X=x, Y=y, ERR=err ; model = F(x, p) ;; Model function ; resid = (y - model)/err ;; Residual calculation (for example) ; ; if n_params() GT 1 then begin ; ; Create derivative and compute derivative array ; requested = dp ; Save original value of DP ; dp = make_array(n_elements(x), n_elements(p), value=x[0]*0) ; ; ; Compute derivative if requested by caller ; for i = 0, n_elements(p)-1 do if requested(i) NE 0 then $ ; dp(*,i) = FGRAD(x, p, i) / err ; endif ; ; return, resid ; END ; ; where FGRAD(x, p, i) is a model function which computes the ; derivative of the model F(x,p) with respect to parameter P(i) at X. ; ; A quirk in the implementation leaves a stray negative sign in the ; definition of DP. The derivative of the *residual* should be ; "-FGRAD(x,p,i) / err" because of how the residual is defined ; ("resid = (data - model) / err"). **HOWEVER** because of the ; implementation quirk, MPFIT expects FGRAD(x,p,i)/err instead, ; i.e. the opposite sign of the gradient of RESID. ; ; Derivatives should be returned in the DP array. DP should be an ; ARRAY(m,n) array, where m is the number of data points and n is the ; number of parameters. -DP[i,j] is the derivative of the ith ; residual with respect to the jth parameter (note the minus sign ; due to the quirk described above). ; ; As noted above, 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. ; ; The derivatives with respect to fixed parameters are ignored; zero ; is an appropriate value to insert for those derivatives. Upon ; input to the user function, DP is set to a vector with the same ; length as P, with a value of 1 for a parameter which is free, and a ; value of zero for a parameter which is fixed (and hence no ; derivative needs to be calculated). This input vector may be ; overwritten as needed. In the example above, the original DP ; vector is saved to a variable called REQUESTED, and used as a mask ; to calculate only those derivatives that are required. ; ; If the data is higher than one dimensional, then the *last* ; dimension should be the parameter dimension. Example: fitting a ; 50x50 image, "dp" should be 50x50xNPAR. ; ; EXPLICIT DERIVATIVES - TESTING and DEBUGGING ; ; For reasonably complicated user functions, the calculation of ; explicit derivatives of the correct sign and magnitude can be ; difficult to get right. A simple sign error can cause MPFIT to be ; confused. MPFIT has a derivative debugging mode which will compute ; the derivatives *both* numerically and explicitly, and compare the ; results. ; ; It is expected that during production usage, derivative debugging ; should be disabled for all parameters. ; ; In order to enable derivative debugging mode, set the following ; PARINFO members for the ith parameter. ; PARINFO[i].MPSIDE = 3 ; Enable explicit derivatives ; PARINFO[i].MPDERIV_DEBUG = 1 ; Enable derivative debugging mode ; PARINFO[i].MPDERIV_RELTOL = ?? ; Relative tolerance for comparison ; PARINFO[i].MPDERIV_ABSTOL = ?? ; Absolute tolerance for comparison ; Note that these settings are maintained on a parameter-by-parameter ; basis using PARINFO, so the user can choose which parameters ; derivatives will be tested. ; ; When .MPDERIV_DEBUG is set, then MPFIT first computes the ; derivative explicitly by requesting them from the user function. ; Then, it computes the derivatives numerically via finite ; differencing, and compares the two values. If the difference ; exceeds a tolerance threshold, then the values are printed out to ; alert the user. The tolerance level threshold contains both a ; relative and an absolute component, and is expressed as, ; ; ABS(DERIV_U - DERIV_N) GE (ABSTOL + RELTOL*ABS(DERIV_U)) ; ; where DERIV_U and DERIV_N are the derivatives computed explicitly ; and numerically, respectively. Appropriate values ; for most users will be: ; ; PARINFO[i].MPDERIV_RELTOL = 1d-3 ;; Suggested relative tolerance ; PARINFO[i].MPDERIV_ABSTOL = 1d-7 ;; Suggested absolute tolerance ; ; although these thresholds may have to be adjusted for a particular ; problem. When the threshold is exceeded, users can expect to see a ; tabular report like this one: ; ; FJAC DEBUG BEGIN ; # IPNT FUNC DERIV_U DERIV_N DIFF_ABS DIFF_REL ; FJAC PARM 2 ; 80 -0.7308 0.04233 0.04233 -5.543E-07 -1.309E-05 ; 99 1.370 0.01417 0.01417 -5.518E-07 -3.895E-05 ; 118 0.07187 -0.01400 -0.01400 -5.566E-07 3.977E-05 ; 137 1.844 -0.04216 -0.04216 -5.589E-07 1.326E-05 ; FJAC DEBUG END ; ; The report will be bracketed by FJAC DEBUG BEGIN/END statements. ; Each parameter will be delimited by the statement FJAC PARM n, ; where n is the parameter number. The columns are, ; ; IPNT - data point number (0 ... M-1) ; FUNC - function value at that point ; DERIV_U - explicit derivative value at that point ; DERIV_N - numerical derivative estimate at that point ; DIFF_ABS - absolute difference = (DERIV_U - DERIV_N) ; DIFF_REL - relative difference = (DIFF_ABS)/(DERIV_U) ; ; When prints appear in this report, it is most important to check ; that the derivatives computed in two different ways have the same ; numerical sign and the same order of magnitude, since these are the ; most common programming mistakes. ; ; A line of this form may also appear ; ; # FJAC_MASK = 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ; ; This line indicates for which parameters explicit derivatives are ; expected. A list of all-1s indicates all explicit derivatives for ; all parameters are requested from the user function. ; ; ; 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 - selector for type of derivative calculation. This ; field can take one of five possible 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) ; 3 - explicit derivative used for this parameter ; ; In the first four cases, the derivative is approximated ; numerically by finite difference, with step size ; H=STEP, where the STEP parameter is defined above. The ; last case, MPSIDE=3, indicates to allow the user ; function to compute the derivative explicitly (see ; section on "EXPLICIT DERIVATIVES"). AUTODERIVATIVE=0 ; overrides this setting for all parameters, and is ; equivalent to MPSIDE=3 for all parameters. For ; MPSIDE=0, the "automatic" one-sided derivative method ; will chose a direction for the finite difference which ; does not violate any constraints. The other methods ; (MPSIDE=-1 or MPSIDE=1) do not perform this check. The ; two-sided method is in principle more precise, but ; requires twice as many function evaluations. Default: ; 0. ; ; .MPDERIV_DEBUG - set this value to 1 to enable debugging of ; user-supplied explicit derivatives (see "TESTING and ; DEBUGGING" section above). In addition, the ; user must enable calculation of explicit derivatives by ; either setting AUTODERIVATIVE=0, or MPSIDE=3 for the ; desired parameters. When this option is enabled, a ; report may be printed to the console, depending on the ; MPDERIV_ABSTOL and MPDERIV_RELTOL settings. ; Default: 0 (no debugging) ; ; ; .MPDERIV_ABSTOL, .MPDERIV_RELTOL - tolerance settings for ; print-out of debugging information, for each parameter ; where debugging is enabled. See "TESTING and ; DEBUGGING" section above for the meanings of these two ; fields. ; ; ; .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, ; and any LIMITS are not obeyed. ; [ 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. ; ; ; RECURSION ; ; Generally, recursion is not allowed. As of version 1.77, MPFIT has ; recursion protection which does not allow a model function to ; itself call MPFIT. Users who wish to perform multi-level ; optimization should investigate the 'EXTERNAL' function evaluation ; methods described below for hard-to-evaluate functions. That ; method places more control in the user's hands. The user can ; design a "recursive" application by taking care. ; ; In most cases the recursion protection should be well-behaved. ; However, if the user is doing debugging, it is possible for the ; protection system to get "stuck." In order to reset it, run the ; procedure: ; MPFIT_RESET_RECURSION ; and the protection system should get "unstuck." It is save to call ; this procedure at any time. ; ; ; COMPATIBILITY ; ; This function is designed to work with IDL 5.0 or greater. ; ; Because TIED parameters and the "(EXTERNAL)" user-model feature use ; the EXECUTE() function, they cannot be used with the free version ; of the IDL Virtual Machine. ; ; ; DETERMINING THE VERSION OF MPFIT ; ; MPFIT is a changing library. Users of MPFIT may also depend on a ; specific version of the library being present. As of version 1.70 ; of MPFIT, a VERSION keyword has been added which allows the user to ; query which version is present. The keyword works like this: ; ; RESULT = MPFIT(/query, VERSION=version) ; ; This call uses the /QUERY keyword to query the version number ; without performing any computations. Users of MPFIT can call this ; method to determine which version is in the IDL path before ; actually using MPFIT to do any numerical work. Upon return, the ; VERSION keyword contains the version number of MPFIT, expressed as ; a string of the form 'X.Y' where X and Y are integers. ; ; Users can perform their own version checking, or use the built-in ; error checking of MPFIT. The MIN_VERSION keyword enforces the ; requested minimum version number. For example, ; ; RESULT = MPFIT(/query, VERSION=version, MIN_VERSION='1.70') ; ; will check whether the accessed version is 1.70 or greater, without ; performing any numerical processing. ; ; The VERSION and MIN_VERSION keywords were added in MPFIT ; version 1.70 and later. If the caller attempts to use the VERSION ; or MIN_VERSION keywords, and an *older* version of the code is ; present in the caller's path, then IDL will throw an 'unknown ; keyword' error. Therefore, in order to be robust, the caller, must ; use exception handling. Here is an example demanding at least ; version 1.70. ; ; MPFIT_OK = 0 & VERSION = '' ; CATCH, CATCHERR ; IF CATCHERR EQ 0 THEN MPFIT_OK = MPFIT(/query, VERSION=version, $ ; MIN_VERSION='1.70') ; CATCH, /CANCEL ; ; IF NOT MPFIT_OK THEN $ ; MESSAGE, 'ERROR: you must have MPFIT version 1.70 or higher in '+$ ; 'your path (found version '+version+')' ; ; Of course, the caller can also do its own version number ; requirements checking. ; ; ; HARD-TO-COMPUTE FUNCTIONS: "EXTERNAL" EVALUATION ; ; The normal mode of operation for MPFIT is for the user to pass a ; function name, and MPFIT will call the user function multiple times ; as it iterates toward a solution. ; ; Some user functions are particularly hard to compute using the ; standard model of MPFIT. Usually these are functions that depend ; on a large amount of external data, and so it is not feasible, or ; at least highly impractical, to have MPFIT call it. In those cases ; it may be possible to use the "(EXTERNAL)" evaluation option. ; ; In this case the user is responsible for making all function *and ; derivative* evaluations. The function and Jacobian data are passed ; in through the EXTERNAL_FVEC and EXTERNAL_FJAC keywords, ; respectively. The user indicates the selection of this option by ; specifying a function name (MYFUNCT) of "(EXTERNAL)". No ; user-function calls are made when EXTERNAL evaluation is being ; used. ; ; ** SPECIAL NOTE ** For the "(EXTERNAL)" case, the quirk noted above ; does not apply. The gradient matrix, EXTERNAL_FJAC, should be ; comparable to "-FGRAD(x,p)/err", which is the *opposite* sign of ; the DP matrix described above. In other words, EXTERNAL_FJAC ; has the same sign as the derivative of EXTERNAL_FVEC, and the ; opposite sign of FGRAD. ; ; At the end of each iteration, control returns to the user, who must ; reevaluate the function at its new parameter values. Users should ; check the return value of the STATUS keyword, where a value of 9 ; indicates the user should supply more data for the next iteration, ; and re-call MPFIT. The user may refrain from calling MPFIT ; further; as usual, STATUS will indicate when the solution has ; converged and no more iterations are required. ; ; Because MPFIT must maintain its own data structures between calls, ; the user must also pass a named variable to the EXTERNAL_STATE ; keyword. This variable must be maintained by the user, but not ; changed, throughout the fitting process. When no more iterations ; are desired, the named variable may be discarded. ; ; ; INPUTS: ; MYFUNCT - a string variable containing the name of the function to ; be minimized. The function should return the weighted ; deviations between the model and the data, as described ; above. ; ; For EXTERNAL evaluation of functions, this parameter ; should be set to a value of "(EXTERNAL)". ; ; START_PARAMS - An one-dimensional 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 (but see PARINFO). 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. ; Exceptions: ; * if /QUERY is set (see QUERY). ; ; ; KEYWORD PARAMETERS: ; ; AUTODERIVATIVE - If this is set, derivatives of the function will ; be computed automatically via a finite ; differencing procedure. If not set, then MYFUNCT ; must provide the explicit derivatives. ; Default: set (=1) ; NOTE: to supply your own explicit derivatives, ; explicitly pass AUTODERIVATIVE=0 ; ; BESTNORM - upon return, the value of the summed squared weighted ; residuals for the returned parameter values, ; i.e. TOTAL(DEVIATES^2). ; ; 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 DEVIATES[I] with respect to ; parameter PARMS[PFREE_INDEX[J]]. Note that since ; deviates are (data-model)*weight, the Jacobian of the ; *deviates* will have the opposite sign from the ; Jacobian of the *model*, and may be scaled by a ; factor. ; ; BEST_RESID - upon return, an array of best-fit deviates. ; ; CALC_FJAC - if set, then calculate the Jacobian and return it in ; BEST_FJAC. If not set, then the return value of ; BEST_FJAC is undefined. ; ; 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. ; ; 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 by the user function. ; ; ERRMSG - a string error or warning message is returned. ; ; EXTERNAL_FVEC - upon input, the function values, evaluated at ; START_PARAMS. This should be an M-vector, where M ; is the number of data points. ; ; EXTERNAL_FJAC - upon input, the Jacobian array of partial ; derivative values. This should be a M x N array, ; where M is the number of data points and N is the ; number of parameters. NOTE: that all FIXED or ; TIED parameters must *not* be included in this ; array. ; ; EXTERNAL_STATE - a named variable to store MPFIT-related state ; information between iterations (used in input and ; output to MPFIT). The user must not manipulate ; or discard this data until the final iteration is ; performed. ; ; FASTNORM - set this keyword to select a faster algorithm to ; compute sum-of-square values internally. For systems ; with large numbers of data points, the standard ; algorithm can become prohibitively slow because it ; cannot be vectorized well. By setting this keyword, ; MPFIT will run faster, but it will be more prone to ; floating point overflows and underflows. Thus, setting ; this keyword may sacrifice some stability in the ; fitting process. ; ; 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. ; ; Consider the following example: ; if FUNCTARGS = { XVAL:[1.D,2,3], YVAL:[1.D,4,9], ; ERRVAL:[1.D,1,1] } ; then the user supplied function should be declared ; like this: ; FUNCTION MYFUNCT, P, XVAL=x, YVAL=y, ERRVAL=err ; ; By default, no extra parameters are passed to the ; user-supplied function, but your function should ; accept *at least* one keyword parameter. [ This is to ; accomodate a limitation in IDL's _EXTRA ; parameter-passing mechanism. ] ; ; 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. ; ; ITERPRINT - The name of an IDL procedure, equivalent to PRINT, ; that ITERPROC will use to render output. ITERPRINT ; should be able to accept at least four positional ; arguments. In addition, it should be able to accept ; the standard FORMAT keyword for output formatting; and ; the UNIT keyword, to redirect output to a logical file ; unit (default should be UNIT=1, standard output). ; These keywords are passed using the ITERARGS keyword ; above. The ITERPRINT procedure must accept the _EXTRA ; keyword. ; NOTE: that much formatting can be handled with the ; MPPRINT and MPFORMAT tags. ; Default: 'MPFIT_DEFPRINT' (default internal formatter) ; ; ITERPROC - The name of a procedure to be called upon each NPRINT ; iteration of the MPFIT routine. ITERPROC is always ; called in the final iteration. It should be declared ; in the following way: ; ; PRO ITERPROC, MYFUNCT, p, iter, fnorm, FUNCTARGS=fcnargs, $ ; PARINFO=parinfo, QUIET=quiet, DOF=dof, PFORMAT=pformat, $ ; UNIT=unit, ... ; ; 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. DOF is the number of degrees of freedom, ; normally the number of points less the number of free ; parameters. See below for documentation of PARINFO. ; PFORMAT is the default parameter value format. UNIT is ; passed on to the ITERPRINT procedure, and should ; indicate the file unit where log output will be sent ; (default: standard output). ; ; 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 ; between -15 and -1 (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. ; ; ITERSTOP - Set this keyword if you wish to be able to stop the ; fitting by hitting the predefined ITERKEYSTOP key on ; the keyboard. This only works if you use the default ; ITERPROC. ; ; ITERKEYSTOP - A keyboard key which will halt the fit (and if ; ITERSTOP is set and the default ITERPROC is used). ; ITERSTOPKEY may either be a one-character string ; with the desired key, or a scalar integer giving the ; ASCII code of the desired key. ; Default: 7b (control-g) ; ; NOTE: the default value of ASCI 7 (control-G) cannot ; be read in some windowing environments, so you must ; change to a printable character like 'q'. ; ; 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 ; ; MIN_VERSION - The minimum requested version number. This must be ; a scalar string of the form returned by the VERSION ; keyword. If the current version of MPFIT does not ; satisfy the minimum requested version number, then, ; MPFIT(/query, min_version='...') returns 0 ; MPFIT(...) returns NAN ; Default: no version number check ; NOTE: MIN_VERSION was added in MPFIT version 1.70 ; ; 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. ; ; NOCATCH - if set, then MPFIT will not perform any error trapping. ; By default (not set), MPFIT will trap errors and report ; them to the caller. This keyword will typically be used ; for debugging. ; ; 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 ; The VERSION output keyword is always set upon return. ; Default: not set. ; ; QUIET - set this keyword when no textual output should be printed ; by MPFIT ; ; RESDAMP - a scalar number, indicating the cut-off value of ; residuals where "damping" will occur. Residuals with ; magnitudes greater than this number will be replaced by ; their logarithm. This partially mitigates the so-called ; large residual problem inherent in least-squares solvers ; (as for the test problem CURVI, http://www.maxthis.com/- ; curviex.htm). A value of 0 indicates no damping. ; Default: 0 ; ; Note: RESDAMP doesn't work with AUTODERIV=0 ; ; 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. ; ; 9 A successful single iteration has been completed, and ; the user must supply another "EXTERNAL" evaluation of ; the function and its derivatives. This status indicator ; is neither an error nor a convergence indicator. ; ; VERSION - upon return, VERSION will be set to the MPFIT internal ; version number. The version number will be a string of ; the form "X.Y" where X is a major revision number and Y ; is a minor revision number. ; NOTE: the VERSION keyword was not present before ; MPFIT version number 1.70, therefore, callers must ; use exception handling when using this keyword. ; ; 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 ; ; ; EXAMPLE: ; ; p0 = [5.7D, 2.2, 500., 1.5, 2000.] ; fa = {X:x, Y:y, ERR:err} ; p = mpfit('MYFUNCT', p0, functargs=fa) ; ; Minimizes sum of squares of MYFUNCT. MYFUNCT is called with the X, ; Y, and ERR keyword parameters that are given by FUNCTARGS. The ; resulting parameter values are returned in p. ; ; ; 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. This value is also returned in the ; STATUS keyword: values of -1 through -15 are reserved error ; codes for the user routines. By default the value of ERROR_CODE ; 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 ) ; 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 dustem_mpfit_dummy ;; Enclose in a procedure so these are not defined in the main level COMPILE_OPT strictarr FORWARD_FUNCTION dustem_mpfit_fdjac2, dustem_mpfit_enorm, dustem_mpfit_lmpar, dustem_mpfit_covar, $ dustem_mpfit, dustem_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 dustem_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 dustem_mpfit_setmachar, double=isdouble COMPILE_OPT strictarr common mpfit_profile, profvals if n_elements(profvals) EQ 0 then dustem_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 dustem_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 dustem_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 dustem_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 dustem_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 dustem_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 = 'DUSTEM_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 dustem_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 = dustem_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 ) ;; ... 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 = dustem_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 = dustem_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 dustem_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 dustem_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] = dustem_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 = dustem_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] = dustem_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 dustem_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 $ r[j:n-1,j] = r[j,j:n-1] x = r[delm] wa = qtb ;; Below may look strange, but it's so we can keep the right precision zero = qtb[0]*0. half = zero + 0.5 quart = zero + 0.25 ;; Eliminate the diagonal matrix d using a givens rotation for j = 0L, n-1 do begin l = ipvt[j] if diag[l] EQ 0 then goto, STORE_RESTORE sdiag[j:*] = 0 sdiag[j] = diag[l] ;; The transformations to eliminate the row of d modify only a ;; single element of (q transpose)*b beyond the first n, which ;; is initially zero. qtbpj = zero for k = j, n-1 do begin if sdiag[k] EQ 0 then goto, ELIM_NEXT_LOOP if abs(r[k,k]) LT abs(sdiag[k]) then begin cotan = r[k,k]/sdiag[k] sine = half/sqrt(quart + quart*cotan*cotan) cosine = sine*cotan endif else begin tang = sdiag[k]/r[k,k] cosine = half/sqrt(quart + quart*tang*tang) sine = cosine*tang endelse ;; Compute the modified diagonal element of r and the ;; modified element of ((q transpose)*b,0). r[k,k] = cosine*r[k,k] + sine*sdiag[k] temp = cosine*wa[k] + sine*qtbpj qtbpj = -sine*wa[k] + cosine*qtbpj wa[k] = temp ;; Accumulate the transformation in the row of s if n GT k+1 then begin temp = cosine*r[k+1:n-1,k] + sine*sdiag[k+1:n-1] sdiag[k+1:n-1] = -sine*r[k+1:n-1,k] + cosine*sdiag[k+1:n-1] r[k+1:n-1,k] = temp endif ELIM_NEXT_LOOP: endfor STORE_RESTORE: sdiag[j] = r[j,j] r[j,j] = x[j] endfor ;; Solve the triangular system for z. If the system is singular ;; then obtain a least squares solution nsing = n wh = where(sdiag EQ 0, ct) if ct GT 0 then begin nsing = wh[0] wa[nsing:*] = 0 endif if nsing GE 1 then begin wa[nsing-1] = wa[nsing-1]/sdiag[nsing-1] ;; Degenerate case ;; *** Reverse loop *** for j=nsing-2,0,-1 do begin sum = total(r[j+1:nsing-1,j]*wa[j+1:nsing-1]) wa[j] = (wa[j]-sum)/sdiag[j] endfor endif ;; Permute the components of z back to components of x x[ipvt] = wa ; profvals.qrsolv = profvals.qrsolv + (systime(1) - prof_start) return end ; ; subroutine lmpar ; ; given an m by n matrix a, an n by n nonsingular diagonal ; matrix d, an m-vector b, and a positive number delta, ; the problem is to determine a value for the parameter ; par such that if x solves the system ; ; a*x = b , sqrt(par)*d*x = 0 , ; ; in the least squares sense, and dxnorm is the euclidean ; norm of d*x, then either par is zero and ; ; (dxnorm-delta) .le. 0.1*delta , ; ; or par is positive and ; ; abs(dxnorm-delta) .le. 0.1*delta . ; ; 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 lmpar expects ; the full upper triangle of r, the permutation matrix p, ; and the first n components of (q transpose)*b. on output ; lmpar also provides an upper triangular matrix s such that ; ; t t t ; p *(a *a + par*d*d)*p = s *s . ; ; s is employed within lmpar and may be of separate interest. ; ; only a few iterations are generally needed for convergence ; of the algorithm. if, however, the limit of 10 iterations ; is reached, then the output par will contain the best ; value obtained so far. ; ; the subroutine statement is ; ; subroutine lmpar(n,r,ldr,ipvt,diag,qtb,delta,par,x,sdiag, ; wa1,wa2) ; ; 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. ; ; delta is a positive input variable which specifies an upper ; bound on the euclidean norm of d*x. ; ; par is a nonnegative variable. on input par contains an ; initial estimate of the levenberg-marquardt parameter. ; on output par contains the final estimate. ; ; x is an output array of length n which contains the least ; squares solution of the system a*x = b, sqrt(par)*d*x = 0, ; for the output par. ; ; sdiag is an output array of length n which contains the ; diagonal elements of the upper triangular matrix s. ; ; wa1 and wa2 are work arrays of length n. ; ; subprograms called ; ; minpack-supplied ... dpmpar,enorm,qrsolv ; ; fortran-supplied ... dabs,dmax1,dmin1,dsqrt ; ; argonne national laboratory. minpack project. march 1980. ; burton s. garbow, kenneth e. hillstrom, jorge j. more ; function dustem_mpfit_lmpar, r, ipvt, diag, qtb, delta, x, sdiag, par=par COMPILE_OPT strictarr common mpfit_machar, machvals common mpfit_profile, profvals ; prof_start = systime(1) MACHEP0 = machvals.machep DWARF = machvals.minnum sz = size(r) m = sz[1] n = sz[2] delm = lindgen(n) * (m+1) ;; Diagonal elements of r ;; Compute and store in x the gauss-newton direction. If the ;; jacobian is rank-deficient, obtain a least-squares solution nsing = n wa1 = qtb rthresh = max(abs(r[delm]))*MACHEP0 wh = where(abs(r[delm]) LT rthresh, ct) if ct GT 0 then begin nsing = wh[0] wa1[wh[0]:*] = 0 endif if nsing GE 1 then begin ;; *** Reverse loop *** for j=nsing-1,0,-1 do begin wa1[j] = wa1[j]/r[j,j] if (j-1 GE 0) then $ wa1[0:(j-1)] = wa1[0:(j-1)] - r[0:(j-1),j]*wa1[j] endfor endif ;; Note: ipvt here is a permutation array x[ipvt] = wa1 ;; Initialize the iteration counter. Evaluate the function at the ;; origin, and test for acceptance of the gauss-newton direction iter = 0L wa2 = diag * x dxnorm = dustem_mpfit_enorm(wa2) fp = dxnorm - delta if fp LE 0.1*delta then goto, TERMINATE ;; If the jacobian is not rank deficient, the newton step provides a ;; lower bound, parl, for the zero of the function. Otherwise set ;; this bound to zero. zero = wa2[0]*0. parl = zero if nsing GE n then begin wa1 = diag[ipvt]*wa2[ipvt]/dxnorm wa1[0] = wa1[0] / r[0,0] ;; Degenerate case for j=1L, n-1 do begin ;; Note "1" here, not zero sum = total(r[0:(j-1),j]*wa1[0:(j-1)]) wa1[j] = (wa1[j] - sum)/r[j,j] endfor temp = dustem_mpfit_enorm(wa1) parl = ((fp/delta)/temp)/temp endif ;; Calculate an upper bound, paru, for the zero of the function for j=0L, n-1 do begin sum = total(r[0:j,j]*qtb[0:j]) wa1[j] = sum/diag[ipvt[j]] endfor gnorm = dustem_mpfit_enorm(wa1) paru = gnorm/delta if paru EQ 0 then paru = DWARF/min([delta,0.1]) ;; If the input par lies outside of the interval (parl,paru), set ;; par to the closer endpoint par = max([par,parl]) par = min([par,paru]) if par EQ 0 then par = gnorm/dxnorm ;; Beginning of an interation ITERATION: iter = iter + 1 ;; Evaluate the function at the current value of par if par EQ 0 then par = max([DWARF, paru*0.001]) temp = sqrt(par) wa1 = temp * diag dustem_mpfit_qrsolv, r, ipvt, wa1, qtb, x, sdiag wa2 = diag*x dxnorm = dustem_mpfit_enorm(wa2) temp = fp fp = dxnorm - delta if (abs(fp) LE 0.1D*delta) $ OR ((parl EQ 0) AND (fp LE temp) AND (temp LT 0)) $ OR (iter EQ 10) then goto, TERMINATE ;; Compute the newton correction wa1 = diag[ipvt]*wa2[ipvt]/dxnorm for j=0L,n-2 do begin wa1[j] = wa1[j]/sdiag[j] wa1[j+1:n-1] = wa1[j+1:n-1] - r[j+1:n-1,j]*wa1[j] endfor wa1[n-1] = wa1[n-1]/sdiag[n-1] ;; Degenerate case temp = dustem_mpfit_enorm(wa1) parc = ((fp/delta)/temp)/temp ;; Depending on the sign of the function, update parl or paru if fp GT 0 then parl = max([parl,par]) if fp LT 0 then paru = min([paru,par]) ;; Compute an improved estimate for par par = max([parl, par+parc]) ;; End of an iteration goto, ITERATION TERMINATE: ;; Termination ; profvals.lmpar = profvals.lmpar + (systime(1) - prof_start) if iter EQ 0 then return, par[0]*0. return, par end ;; Procedure to tie one parameter to another. pro dustem_mpfit_tie, p, _ptied COMPILE_OPT strictarr if n_elements(_ptied) EQ 0 then return if n_elements(_ptied) EQ 1 then if _ptied[0] EQ '' then return for _i = 0L, n_elements(_ptied)-1 do begin if _ptied[_i] EQ '' then goto, NEXT_TIE _cmd = 'p['+strtrim(_i,2)+'] = '+_ptied[_i] _err = execute(_cmd) if _err EQ 0 then begin message, 'ERROR: Tied expression "'+_cmd+'" failed.' return endif NEXT_TIE: endfor end ;; Default print procedure pro dustem_mpfit_defprint, p1, p2, p3, p4, p5, p6, p7, p8, p9, p10, $ p11, p12, p13, p14, p15, p16, p17, p18, $ format=format, unit=unit0, _EXTRA=extra COMPILE_OPT strictarr if n_elements(unit0) EQ 0 then unit = -1 else unit = round(unit0[0]) if n_params() EQ 0 then printf, unit, '' $ else if n_params() EQ 1 then printf, unit, p1, format=format $ else if n_params() EQ 2 then printf, unit, p1, p2, format=format $ else if n_params() EQ 3 then printf, unit, p1, p2, p3, format=format $ else if n_params() EQ 4 then printf, unit, p1, p2, p4, format=format return end ;; Default procedure to be called every iteration. It simply prints ;; the parameter values. pro dustem_mpfit_defiter, fcn, x, iter, fnorm, FUNCTARGS=fcnargs, $ quiet=quiet, iterstop=iterstop, iterkeybyte=iterkeybyte, $ parinfo=parinfo, iterprint=iterprint0, $ format=fmt, pformat=pformat, dof=dof0, _EXTRA=iterargs COMPILE_OPT strictarr common mpfit_error, mperr mperr = 0 if keyword_set(quiet) then goto, DO_ITERSTOP if n_params() EQ 3 then begin fvec = dustem_mpfit_call(fcn, x, _EXTRA=fcnargs) fnorm = dustem_mpfit_enorm(fvec)^2 endif ;; Determine which parameters to print nprint = n_elements(x) iprint = lindgen(nprint) if n_elements(iterprint0) EQ 0 then iterprint = 'DUSTEM_MPFIT_DEFPRINT' $ else iterprint = strtrim(iterprint0[0],2) if n_elements(dof0) EQ 0 then dof = 1L else dof = floor(dof0[0]) call_procedure, iterprint, iter, fnorm, dof, $ format='("Iter ",I6," CHI-SQUARE = ",G15.8," DOF = ",I0)', $ _EXTRA=iterargs if n_elements(fmt) GT 0 then begin call_procedure, iterprint, x, format=fmt, _EXTRA=iterargs endif else begin if n_elements(pformat) EQ 0 then pformat = '(G40.6)' parname = 'P('+strtrim(iprint,2)+')' pformats = strarr(nprint) + pformat if n_elements(parinfo) GT 0 then begin parinfo_tags = tag_names(parinfo) wh = where(parinfo_tags EQ 'PARNAME', ct) if ct EQ 1 then begin wh = where(parinfo.parname NE '', ct) if ct GT 0 then $ parname[wh] = strmid(parinfo[wh].parname,0,25) endif wh = where(parinfo_tags EQ 'MPPRINT', ct) if ct EQ 1 then begin iprint = where(parinfo.mpprint EQ 1, nprint) if nprint EQ 0 then goto, DO_ITERSTOP endif wh = where(parinfo_tags EQ 'MPFORMAT', ct) if ct EQ 1 then begin wh = where(parinfo.mpformat NE '', ct) if ct GT 0 then pformats[wh] = parinfo[wh].mpformat endif endif for i = 0L, nprint-1 do begin call_procedure, iterprint, parname[iprint[i]], x[iprint[i]], $ format='(" ",A0," = ",'+pformats[iprint[i]]+')', $ _EXTRA=iterargs endfor endelse DO_ITERSTOP: if n_elements(iterkeybyte) EQ 0 then iterkeybyte = 7b if keyword_set(iterstop) then begin k = get_kbrd(0) if k EQ string(iterkeybyte[0]) then begin message, 'WARNING: minimization not complete', /info print, 'Do you want to terminate this procedure? (y/n)', $ format='(A,$)' k = '' read, k if strupcase(strmid(k,0,1)) EQ 'Y' then begin message, 'WARNING: Procedure is terminating.', /info mperr = -1 endif endif endif return end ;; Procedure to parse the parameter values in PARINFO pro dustem_mpfit_parinfo, parinfo, tnames, tag, values, default=def, status=status, $ n_param=n COMPILE_OPT strictarr status = 0 if n_elements(n) EQ 0 then n = n_elements(parinfo) if n EQ 0 then begin if n_elements(def) EQ 0 then return values = def status = 1 return endif if n_elements(parinfo) EQ 0 then goto, DO_DEFAULT if n_elements(tnames) EQ 0 then tnames = tag_names(parinfo) wh = where(tnames EQ tag, ct) if ct EQ 0 then begin DO_DEFAULT: if n_elements(def) EQ 0 then return values = make_array(n, value=def[0]) values[0] = def endif else begin values = parinfo.(wh[0]) np = n_elements(parinfo) nv = n_elements(values) values = reform(values[*], nv/np, np) endelse status = 1 return end ; ********** ; ; subroutine covar ; ; given an m by n matrix a, the problem is to determine ; the covariance matrix corresponding to a, defined as ; ; t ; inverse(a *a) . ; ; 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 covar expects ; the full upper triangle of r and the permutation matrix p. ; the covariance matrix is then computed as ; ; t t ; p*inverse(r *r)*p . ; ; if a is nearly rank deficient, it may be desirable to compute ; the covariance matrix corresponding to the linearly independent ; columns of a. to define the numerical rank of a, covar uses ; the tolerance tol. if l is the largest integer such that ; ; abs(r(l,l)) .gt. tol*abs(r(1,1)) , ; ; then covar computes the covariance matrix corresponding to ; the first l columns of r. for k greater than l, column ; and row ipvt(k) of the covariance matrix are set to zero. ; ; the subroutine statement is ; ; subroutine covar(n,r,ldr,ipvt,tol,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 ; r contains the square symmetric covariance matrix. ; ; 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. ; ; tol is a nonnegative input variable used to define the ; numerical rank of a in the manner described above. ; ; wa is a work array of length n. ; ; subprograms called ; ; fortran-supplied ... dabs ; ; argonne national laboratory. minpack project. august 1980. ; burton s. garbow, kenneth e. hillstrom, jorge j. more ; ; ********** function dustem_mpfit_covar, rr, ipvt, tol=tol COMPILE_OPT strictarr sz = size(rr) if sz[0] NE 2 then begin message, 'ERROR: r must be a two-dimensional matrix' return, -1L endif n = sz[1] if n NE sz[2] then begin message, 'ERROR: r must be a square matrix' return, -1L endif zero = rr[0] * 0. one = zero + 1. if n_elements(ipvt) EQ 0 then ipvt = lindgen(n) r = rr r = reform(rr, n, n, /overwrite) ;; Form the inverse of r in the full upper triangle of r l = -1L if n_elements(tol) EQ 0 then tol = one*1.E-14 tolr = tol * abs(r[0,0]) for k = 0L, n-1 do begin if abs(r[k,k]) LE tolr then goto, INV_END_LOOP r[k,k] = one/r[k,k] for j = 0L, k-1 do begin temp = r[k,k] * r[j,k] r[j,k] = zero r[0,k] = r[0:j,k] - temp*r[0:j,j] endfor l = k endfor INV_END_LOOP: ;; Form the full upper triangle of the inverse of (r transpose)*r ;; in the full upper triangle of r if l GE 0 then $ for k = 0L, l do begin for j = 0L, k-1 do begin temp = r[j,k] r[0,j] = r[0:j,j] + temp*r[0:j,k] endfor temp = r[k,k] r[0,k] = temp * r[0:k,k] endfor ;; Form the full lower triangle of the covariance matrix ;; in the strict lower triangle of r and in wa wa = replicate(r[0,0], n) for j = 0L, n-1 do begin jj = ipvt[j] sing = j GT l for i = 0L, j do begin if sing then r[i,j] = zero ii = ipvt[i] if ii GT jj then r[ii,jj] = r[i,j] if ii LT jj then r[jj,ii] = r[i,j] endfor wa[jj] = r[j,j] endfor ;; Symmetrize the covariance matrix in r for j = 0L, n-1 do begin r[0:j,j] = r[j,0:j] r[j,j] = wa[j] endfor return, r end ;; Parse the RCSID revision number function dustem_mpfit_revision ;; NOTE: this string is changed every time an RCS check-in occurs revision = '$Revision: 1.82 $' ;; Parse just the version number portion revision = stregex(revision,'\$'+'Revision: *([0-9.]+) *'+'\$',/extract,/sub) revision = revision[1] return, revision end ;; Parse version numbers of the form 'X.Y' function dustem_mpfit_parse_version, version sz = size(version) if sz[sz[0]+1] NE 7 then return, 0 s = stregex(version[0], '^([0-9]+)\.([0-9]+)$', /extract,/sub) if s[0] NE version[0] then return, 0 return, long(s[1:2]) end ;; Enforce a minimum version number function dustem_mpfit_min_version, version, min_version mv = dustem_mpfit_parse_version(min_version) if mv[0] EQ 0 then return, 0 v =dustem_mpfit_parse_version(version) ;; Compare version components if v[0] LT mv[0] then return, 0 if v[1] LT mv[1] then return, 0 return, 1 end ; Manually reset recursion fencepost if the user gets in trouble pro dustem_mpfit_reset_recursion common mpfit_fencepost, mpfit_fencepost_active mpfit_fencepost_active = 0 end ; ********** ; ; subroutine lmdif ; ; the purpose of lmdif is to minimize the sum of the squares of ; m nonlinear functions in n variables by a modification of ; the levenberg-marquardt algorithm. the user must provide a ; subroutine which calculates the functions. the jacobian is ; then calculated by a forward-difference approximation. ; ; the subroutine statement is ; ; subroutine lmdif(fcn,m,n,x,fvec,ftol,xtol,gtol,maxfev,epsfcn, ; diag,mode,factor,nprint,info,nfev,fjac, ; ldfjac,ipvt,qtf,wa1,wa2,wa3,wa4) ; ; where ; ; fcn is the name of the user-supplied subroutine which ; calculates the functions. fcn must be declared ; in an external statement in the user calling ; program, and should be written as follows. ; ; subroutine fcn(m,n,x,fvec,iflag) ; integer m,n,iflag ; double precision x(n),fvec(m) ; ---------- ; calculate the functions at x and ; return this vector in fvec. ; ---------- ; return ; end ; ; the value of iflag should not be changed by fcn unless ; the user wants to terminate execution of lmdif. ; in this case set iflag to a negative integer. ; ; m is a positive integer input variable set to the number ; of functions. ; ; n is a positive integer input variable set to the number ; of variables. n must not exceed m. ; ; x is an array of length n. on input x must contain ; an initial estimate of the solution vector. on output x ; contains the final estimate of the solution vector. ; ; fvec is an output array of length m which contains ; the functions evaluated at the output x. ; ; ftol is a nonnegative input variable. termination ; occurs when both the actual and predicted relative ; reductions in the sum of squares are at most ftol. ; therefore, ftol measures the relative error desired ; in the sum of squares. ; ; xtol is a nonnegative input variable. termination ; occurs when the relative error between two consecutive ; iterates is at most xtol. therefore, xtol measures the ; relative error desired in the approximate solution. ; ; gtol is 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. therefore, gtol measures the orthogonality ; desired between the function vector and the columns ; of the jacobian. ; ; maxfev is a positive integer input variable. termination ; occurs when the number of calls to fcn is at least ; maxfev by the end of an iteration. ; ; epsfcn is an input variable used in determining a suitable ; step length for the forward-difference approximation. this ; approximation assumes that the relative errors in the ; functions are of the order of epsfcn. if epsfcn is less ; than the machine precision, it is assumed that the relative ; errors in the functions are of the order of the machine ; precision. ; ; diag is an array of length n. if mode = 1 (see ; below), diag is internally set. if mode = 2, diag ; must contain positive entries that serve as ; multiplicative scale factors for the variables. ; ; mode is an integer input variable. if mode = 1, the ; variables will be scaled internally. if mode = 2, ; the scaling is specified by the input diag. other ; values of mode are equivalent to mode = 1. ; ; factor is a positive input variable used in determining the ; initial step bound. this bound is set to the product of ; factor and the euclidean norm of diag*x if nonzero, or else ; to factor itself. in most cases factor should lie in the ; interval (.1,100.). 100. is a generally recommended value. ; ; nprint is an integer input variable that enables controlled ; printing of iterates if it is positive. in this case, ; fcn is called with iflag = 0 at the beginning of the first ; iteration and every nprint iterations thereafter and ; immediately prior to return, with x and fvec available ; for printing. if nprint is not positive, no special calls ; of fcn with iflag = 0 are made. ; ; info is an integer output variable. if the user has ; terminated execution, info is set to the (negative) ; value of iflag. see description of fcn. otherwise, ; info is set as follows. ; ; info = 0 improper input parameters. ; ; info = 1 both actual and predicted relative reductions ; in the sum of squares are at most ftol. ; ; info = 2 relative error between two consecutive iterates ; is at most xtol. ; ; info = 3 conditions for info = 1 and info = 2 both hold. ; ; info = 4 the cosine of the angle between fvec and any ; column of the jacobian is at most gtol in ; absolute value. ; ; info = 5 number of calls to fcn has reached or ; exceeded maxfev. ; ; info = 6 ftol is too small. no further reduction in ; the sum of squares is possible. ; ; info = 7 xtol is too small. no further improvement in ; the approximate solution x is possible. ; ; info = 8 gtol is too small. fvec is orthogonal to the ; columns of the jacobian to machine precision. ; ; nfev is an integer output variable set to the number of ; calls to fcn. ; ; fjac is an output m by n array. the upper n by n submatrix ; of fjac contains an upper triangular matrix r with ; diagonal elements of nonincreasing magnitude such that ; ; t t t ; p *(jac *jac)*p = r *r, ; ; where p is a permutation matrix and jac is the final ; calculated jacobian. column j of p is column ipvt(j) ; (see below) of the identity matrix. the lower trapezoidal ; part of fjac contains information generated during ; the computation of r. ; ; ldfjac is a positive integer input variable not less than m ; which specifies the leading dimension of the array fjac. ; ; ipvt is an integer output array of length n. ipvt ; defines a permutation matrix p such that jac*p = q*r, ; where jac is the final calculated jacobian, q is ; orthogonal (not stored), and r is upper triangular ; with diagonal elements of nonincreasing magnitude. ; column j of p is column ipvt(j) of the identity matrix. ; ; qtf is an output array of length n which contains ; the first n elements of the vector (q transpose)*fvec. ; ; wa1, wa2, and wa3 are work arrays of length n. ; ; wa4 is a work array of length m. ; ; subprograms called ; ; user-supplied ...... fcn ; ; minpack-supplied ... dpmpar,enorm,fdjac2,lmpar,qrfac ; ; fortran-supplied ... dabs,dmax1,dmin1,dsqrt,mod ; ; argonne national laboratory. minpack project. march 1980. ; burton s. garbow, kenneth e. hillstrom, jorge j. more ; ; ********** function dustem_mpfit, fcn, xall, FUNCTARGS=fcnargs, SCALE_FCN=scalfcn, $ ftol=ftol0, xtol=xtol0, gtol=gtol0, epsfcn=epsfcn, $ resdamp=damp0, $ nfev=nfev, maxiter=maxiter, errmsg=errmsg, $ factor=factor0, nprint=nprint0, STATUS=info, $ iterproc=iterproc0, iterargs=iterargs, iterstop=ss,$ iterkeystop=iterkeystop, $ niter=iter, nfree=nfree, npegged=npegged, dof=dof, $ diag=diag, rescale=rescale, autoderivative=autoderiv0, $ pfree_index=ifree, $ perror=perror, covar=covar, nocovar=nocovar, $ bestnorm=fnorm, best_resid=fvec, $ best_fjac=output_fjac, calc_fjac=calc_fjac, $ parinfo=parinfo, quiet=quiet, nocatch=nocatch, $ fastnorm=fastnorm0, proc=proc, query=query, $ external_state=state, external_init=extinit, $ external_fvec=efvec, external_fjac=efjac, $ version=version, min_version=min_version0 COMPILE_OPT strictarr info = 0L errmsg = '' ;; Compute the revision number, to be returned in the VERSION and ;; QUERY keywords. common mpfit_revision_common, mpfit_revision_str if n_elements(mpfit_revision_str) EQ 0 then $ mpfit_revision_str = dustem_mpfit_revision() version = mpfit_revision_str if keyword_set(query) then begin if n_elements(min_version0) GT 0 then $ if dustem_mpfit_min_version(version, min_version0[0]) EQ 0 then $ return, 0 return, 1 endif if n_elements(min_version0) GT 0 then $ if mpfit_min_version(version, min_version0[0]) EQ 0 then begin message, 'ERROR: minimum required version '+min_version0[0]+' not satisfied', /info return, !values.d_nan endif if n_params() EQ 0 then begin message, "USAGE: PARMS = DUSTEM_MPFIT('MYFUNCT', START_PARAMS, ... )", /info return, !values.d_nan endif ;; Use of double here not a problem since f/x/gtol are all only used ;; in comparisons if n_elements(ftol0) EQ 0 then ftol = 1.D-10 else ftol = ftol0[0] if n_elements(xtol0) EQ 0 then xtol = 1.D-10 else xtol = xtol0[0] if n_elements(gtol0) EQ 0 then gtol = 1.D-10 else gtol = gtol0[0] if n_elements(factor0) EQ 0 then factor = 100. else factor = factor0[0] if n_elements(nprint0) EQ 0 then nprint = 1 else nprint = nprint0[0] if n_elements(iterproc0) EQ 0 then iterproc = 'DUSTEM_MPFIT_DEFITER' else iterproc = iterproc0[0] if n_elements(autoderiv0) EQ 0 then autoderiv = 1 else autoderiv = autoderiv0[0] if n_elements(fastnorm0) EQ 0 then fastnorm = 0 else fastnorm = fastnorm0[0] if n_elements(damp0) EQ 0 then damp = 0 else damp = damp0[0] ;; These are special configuration parameters that can't be easily ;; passed by MPFIT directly. ;; FASTNORM - decide on which sum-of-squares technique to use (1) ;; is fast, (0) is slower ;; PROC - user routine is a procedure (1) or function (0) ;; QANYTIED - set to 1 if any parameters are TIED, zero if none ;; PTIED - array of strings, one for each parameter common mpfit_config, mpconfig mpconfig = {fastnorm: keyword_set(fastnorm), proc: 0, nfev: 0L, damp: damp} common mpfit_machar, machvals iflag = 0L catch_msg = 'in DUSTEM_MPFIT' nfree = 0L npegged = 0L dof = 0L output_fjac = 0L ;; Set up a persistent fencepost that prevents recursive calls common mpfit_fencepost, mpfit_fencepost_active if n_elements(mpfit_fencepost_active) EQ 0 then mpfit_fencepost_active = 0 if mpfit_fencepost_active then begin errmsg = 'ERROR: recursion detected; you cannot run DUSTEM_MPFIT recursively' goto, TERMINATE endif ;; Only activate the fencepost if we are not in debugging mode if NOT keyword_set(nocatch) then mpfit_fencepost_active = 1 ;; Parameter damping doesn't work when user is providing their own ;; gradients. if damp NE 0 AND NOT keyword_set(autoderiv) then begin errmsg = 'ERROR: keywords DAMP and AUTODERIV are mutually exclusive' goto, TERMINATE endif ;; Process the ITERSTOP and ITERKEYSTOP keywords, and turn this into ;; a set of keywords to pass to MPFIT_DEFITER. if strupcase(iterproc) EQ 'DSUTEM_MPFIT_DEFITER' AND n_elements(iterargs) EQ 0 $ AND keyword_set(ss) then begin if n_elements(iterkeystop) GT 0 then begin sz = size(iterkeystop) tp = sz[sz[0]+1] if tp EQ 7 then begin ;; String - convert first char to byte iterkeybyte = (byte(iterkeystop[0]))[0] endif if (tp GE 1 AND tp LE 3) OR (tp GE 12 AND tp LE 15) then begin ;; Integer - convert to byte iterkeybyte = byte(iterkeystop[0]) endif if n_elements(iterkeybyte) EQ 0 then begin errmsg = 'ERROR: ITERKEYSTOP must be either a BYTE or STRING' goto, TERMINATE endif iterargs = {iterstop: 1, iterkeybyte: iterkeybyte} endif else begin iterargs = {iterstop: 1, iterkeybyte: 7b} endelse endif ;; Handle error conditions gracefully if NOT keyword_set(nocatch) then begin catch, catcherror if catcherror NE 0 then begin ;; An error occurred!!! catch, /cancel mpfit_fencepost_active = 0 err_string = ''+!error_state.msg message, /cont, 'Error detected while '+catch_msg+':' message, /cont, err_string message, /cont, 'Error condition detected. Returning to MAIN level.' if errmsg EQ '' then $ errmsg = 'Error detected while '+catch_msg+': '+err_string if info EQ 0 then info = -18 return, !values.d_nan endif endif mpconfig = create_struct(mpconfig, 'NOCATCH', keyword_set(nocatch)) ;; Parse FCN function name - be sure it is a scalar string sz = size(fcn) if sz[0] NE 0 then begin FCN_NAME: errmsg = 'ERROR: MYFUNCT must be a scalar string' goto, TERMINATE endif if sz[sz[0]+1] NE 7 then goto, FCN_NAME isext = 0 if fcn EQ '(EXTERNAL)' then begin if n_elements(efvec) EQ 0 OR n_elements(efjac) EQ 0 then begin errmsg = 'ERROR: when using EXTERNAL function, EXTERNAL_FVEC '+$ 'and EXTERNAL_FJAC must be defined' goto, TERMINATE endif nv = n_elements(efvec) nj = n_elements(efjac) if (nj MOD nv) NE 0 then begin errmsg = 'ERROR: the number of values in EXTERNAL_FJAC must be '+ $ 'a multiple of the number of values in EXTERNAL_FVEC' goto, TERMINATE endif isext = 1 endif ;; Parinfo: ;; --------------- STARTING/CONFIG INFO (passed in to routine, not changed) ;; .value - starting value for parameter ;; .fixed - parameter is fixed ;; .limited - a two-element array, if parameter is bounded on ;; lower/upper side ;; .limits - a two-element array, lower/upper parameter bounds, if ;; limited vale is set ;; .step - step size in Jacobian calc, if greater than zero catch_msg = 'parsing input parameters' ;; Parameters can either be stored in parinfo, or x. Parinfo takes ;; precedence if it exists. if n_elements(xall) EQ 0 AND n_elements(parinfo) EQ 0 then begin errmsg = 'ERROR: must pass parameters in P or PARINFO' goto, TERMINATE endif ;; Be sure that PARINFO is of the right type if n_elements(parinfo) GT 0 then begin ;; Make sure the array is 1-D parinfo = parinfo[*] parinfo_size = size(parinfo) if parinfo_size[parinfo_size[0]+1] NE 8 then begin errmsg = 'ERROR: PARINFO must be a structure.' goto, TERMINATE endif if n_elements(xall) GT 0 AND n_elements(xall) NE n_elements(parinfo) $ then begin errmsg = 'ERROR: number of elements in PARINFO and P must agree' goto, TERMINATE endif endif ;; If the parameters were not specified at the command line, then ;; extract them from PARINFO if n_elements(xall) EQ 0 then begin dustem_mpfit_parinfo, parinfo, tagnames, 'VALUE', xall, status=status if status EQ 0 then begin errmsg = 'ERROR: either P or PARINFO[*].VALUE must be supplied.' goto, TERMINATE endif sz = size(xall) ;; Convert to double if parameters are not float or double if sz[sz[0]+1] NE 4 AND sz[sz[0]+1] NE 5 then $ xall = double(xall) endif xall = xall[*] ;; Make sure the array is 1-D npar = n_elements(xall) zero = xall[0] * 0. one = zero + 1. fnorm = -one fnorm1 = -one ;; TIED parameters? dustem_mpfit_parinfo, parinfo, tagnames, 'TIED', ptied, default='', n=npar ptied = strtrim(ptied, 2) wh = where(ptied NE '', qanytied) qanytied = qanytied GT 0 mpconfig = create_struct(mpconfig, 'QANYTIED', qanytied, 'PTIED', ptied) ;; FIXED parameters ? dustem_mpfit_parinfo, parinfo, tagnames, 'FIXED', pfixed, default=0, n=npar pfixed = pfixed EQ 1 pfixed = pfixed OR (ptied NE '');; Tied parameters are also effectively fixed ;; Finite differencing step, absolute and relative, and sidedness of deriv. dustem_mpfit_parinfo, parinfo, tagnames, 'STEP', step, default=zero, n=npar dustem_mpfit_parinfo, parinfo, tagnames, 'RELSTEP', dstep, default=zero, n=npar dustem_mpfit_parinfo, parinfo, tagnames, 'MPSIDE', dside, default=0, n=npar ;; Debugging parameters dustem_mpfit_parinfo, parinfo, tagnames, 'MPDERIV_DEBUG', ddebug, default=0, n=npar dustem_mpfit_parinfo, parinfo, tagnames, 'MPDERIV_RELTOL', ddrtol, default=zero, n=npar dustem_mpfit_parinfo, parinfo, tagnames, 'MPDERIV_ABSTOL', ddatol, default=zero, n=npar ;; Maximum and minimum steps allowed to be taken in one iteration dustem_mpfit_parinfo, parinfo, tagnames, 'MPMAXSTEP', maxstep, default=zero, n=npar dustem_mpfit_parinfo, parinfo, tagnames, 'MPMINSTEP', minstep, default=zero, n=npar qmin = minstep * 0 ;; Remove minstep for now!! qmax = maxstep NE 0 wh = where(qmin AND qmax AND maxstep LT minstep, ct) if ct GT 0 then begin errmsg = 'ERROR: MPMINSTEP is greater than MPMAXSTEP' goto, TERMINATE endif ;; Finish up the free parameters ifree = where(pfixed NE 1, nfree) if nfree EQ 0 then begin errmsg = 'ERROR: no free parameters' goto, TERMINATE endif ;; An external Jacobian must be checked against the number of ;; parameters if isext then begin if (nj/nv) NE nfree then begin errmsg = string(nv, nfree, nfree, $ format=('("ERROR: EXTERNAL_FJAC must be a ",I0," x ",I0,' + $ '" array, where ",I0," is the number of free parameters")')) goto, TERMINATE endif endif ;; Compose only VARYING parameters xnew = xall ;; xnew is the set of parameters to be returned x = xnew[ifree] ;; x is the set of free parameters ; Same for min/max step diagnostics qmin = qmin[ifree] & minstep = minstep[ifree] qmax = qmax[ifree] & maxstep = maxstep[ifree] wh = where(qmin OR qmax, ct) qminmax = ct GT 0 ;; LIMITED parameters ? dustem_mpfit_parinfo, parinfo, tagnames, 'LIMITED', limited, status=st1 dustem_mpfit_parinfo, parinfo, tagnames, 'LIMITS', limits, status=st2 if st1 EQ 1 AND st2 EQ 1 then begin ;; Error checking on limits in parinfo wh = where((limited[0,*] AND xall LT limits[0,*]) OR $ (limited[1,*] AND xall GT limits[1,*]), ct) if ct GT 0 then begin errmsg = 'ERROR: parameters are not within PARINFO limits' goto, TERMINATE endif wh = where(limited[0,*] AND limited[1,*] AND $ limits[0,*] GE limits[1,*] AND $ pfixed EQ 0, ct) if ct GT 0 then begin errmsg = 'ERROR: PARINFO parameter limits are not consistent' goto, TERMINATE endif ;; Transfer structure values to local variables qulim = limited[1, ifree] ulim = limits [1, ifree] qllim = limited[0, ifree] llim = limits [0, ifree] wh = where(qulim OR qllim, ct) if ct GT 0 then qanylim = 1 else qanylim = 0 endif else begin ;; Fill in local variables with dummy values qulim = lonarr(nfree) ulim = x * 0. qllim = qulim llim = x * 0. qanylim = 0 endelse ;; Initialize the number of parameters pegged at a hard limit value wh = where((qulim AND (x EQ ulim)) OR (qllim AND (x EQ llim)), npegged) n = n_elements(x) if n_elements(maxiter) EQ 0 then maxiter = 200L ;; Check input parameters for errors if (n LE 0) OR (ftol LE 0) OR (xtol LE 0) OR (gtol LE 0) $ OR (maxiter LT 0) OR (factor LE 0) then begin errmsg = 'ERROR: input keywords are inconsistent' goto, TERMINATE endif if keyword_set(rescale) then begin errmsg = 'ERROR: DIAG parameter scales are inconsistent' if n_elements(diag) LT n then goto, TERMINATE wh = where(diag LE 0, ct) if ct GT 0 then goto, TERMINATE errmsg = '' endif if n_elements(state) NE 0 AND NOT keyword_set(extinit) then begin szst = size(state) if szst[szst[0]+1] NE 8 then begin errmsg = 'EXTERNAL_STATE keyword was not preserved' status = 0 goto, TERMINATE endif if nfree NE n_elements(state.ifree) then begin BAD_IFREE: errmsg = 'Number of free parameters must not change from one '+$ 'external iteration to the next' status = 0 goto, TERMINATE endif wh = where(ifree NE state.ifree, ct) if ct GT 0 then goto, BAD_IFREE tnames = tag_names(state) for i = 0L, n_elements(tnames)-1 do begin dummy = execute(tnames[i]+' = state.'+tnames[i]) endfor wa4 = reform(efvec, n_elements(efvec)) goto, RESUME_FIT endif common mpfit_error, mperr if NOT isext then begin mperr = 0 catch_msg = 'calling '+fcn fvec = dustem_mpfit_call(fcn, xnew, _EXTRA=fcnargs) iflag = mperr if iflag LT 0 then begin errmsg = 'ERROR: first call to "'+fcn+'" failed' goto, TERMINATE endif endif else begin fvec = reform(efvec, n_elements(efvec)) endelse catch_msg = 'calling DUSTEM_MPFIT_SETMACHAR' sz = size(fvec[0]) isdouble = (sz[sz[0]+1] EQ 5) dustem_mpfit_setmachar, double=isdouble common mpfit_profile, profvals ; prof_start = systime(1) MACHEP0 = machvals.machep DWARF = machvals.minnum szx = size(x) ;; The parameters and the squared deviations should have the same ;; type. Otherwise the MACHAR-based evaluation will fail. catch_msg = 'checking parameter data' tp = szx[szx[0]+1] if tp NE 4 AND tp NE 5 then begin if NOT keyword_set(quiet) then begin message, 'WARNING: input parameters must be at least FLOAT', /info message, ' (converting parameters to FLOAT)', /info endif x = float(x) xnew = float(x) szx = size(x) endif if isdouble AND tp NE 5 then begin if NOT keyword_set(quiet) then begin message, 'WARNING: data is DOUBLE but parameters are FLOAT', /info message, ' (converting parameters to DOUBLE)', /info endif x = double(x) xnew = double(xnew) endif m = n_elements(fvec) if (m LT n) then begin errmsg = 'ERROR: number of parameters must not exceed data' goto, TERMINATE endif fnorm = dustem_mpfit_enorm(fvec) ;; Initialize Levelberg-Marquardt parameter and iteration counter par = zero iter = 1L qtf = x * 0. ;; Beginning of the outer loop OUTER_LOOP: ;; If requested, call fcn to enable printing of iterates xnew[ifree] = x if qanytied then dustem_mpfit_tie, xnew, ptied dof = (n_elements(fvec) - nfree) > 1L if nprint GT 0 AND iterproc NE '' then begin catch_msg = 'calling '+iterproc iflag = 0L if (iter-1) MOD nprint EQ 0 then begin mperr = 0 xnew0 = xnew call_procedure, iterproc, fcn, xnew, iter, fnorm^2, $ FUNCTARGS=fcnargs, parinfo=parinfo, quiet=quiet, $ dof=dof, _EXTRA=iterargs iflag = mperr ;; Check for user termination if iflag LT 0 then begin errmsg = 'WARNING: premature termination by "'+iterproc+'"' goto, TERMINATE endif ;; If parameters were changed (grrr..) then re-tie if max(abs(xnew0-xnew)) GT 0 then begin if qanytied then dustem_mpfit_tie, xnew, ptied x = xnew[ifree] endif endif endif ;; Calculate the jacobian matrix iflag = 2 if NOT isext then begin catch_msg = 'calling DUSTEM_MPFIT_FDJAC2' ;; NOTE! If you change this call then change the one during ;; clean-up as well! fjac = dustem_mpfit_fdjac2(fcn, x, fvec, step, qulim, ulim, dside, $ iflag=iflag, epsfcn=epsfcn, $ autoderiv=autoderiv, dstep=dstep, $ FUNCTARGS=fcnargs, ifree=ifree, xall=xnew, $ deriv_debug=ddebug, deriv_reltol=ddrtol, deriv_abstol=ddatol) if iflag LT 0 then begin errmsg = 'WARNING: premature termination by FDJAC2' goto, TERMINATE endif endif else begin fjac = reform(efjac,n_elements(fvec),npar, /overwrite) endelse ;; Rescale the residuals and gradient, for use with "alternative" ;; statistics such as the Cash statistic. catch_msg = 'prescaling residuals and gradient' if n_elements(scalfcn) GT 0 then begin call_procedure, strtrim(scalfcn[0],2), fvec, fjac endif ;; Determine if any of the parameters are pegged at the limits npegged = 0L if qanylim then begin catch_msg = 'zeroing derivatives of pegged parameters' whlpeg = where(qllim AND (x EQ llim), nlpeg) whupeg = where(qulim AND (x EQ ulim), nupeg) npegged = nlpeg + nupeg ;; See if any "pegged" values should keep their derivatives if (nlpeg GT 0) then begin ;; Total derivative of sum wrt lower pegged parameters ;; Note: total(fvec*fjac[*,i]) is d(CHI^2)/dX[i] for i = 0L, nlpeg-1 do begin sum = total(fvec * fjac[*,whlpeg[i]]) if sum GT 0 then fjac[*,whlpeg[i]] = 0 endfor endif if (nupeg GT 0) then begin ;; Total derivative of sum wrt upper pegged parameters for i = 0L, nupeg-1 do begin sum = total(fvec * fjac[*,whupeg[i]]) if sum LT 0 then fjac[*,whupeg[i]] = 0 endfor endif endif ;; Save a copy of the Jacobian if the user requests it... if keyword_set(calc_fjac) then output_fjac = fjac ;; ===================== ;; Compute the QR factorization of the jacobian catch_msg = 'calling DUSTEM_MPFIT_QRFAC' ;; IN: Jacobian ;; OUT: Hh Vects Permutation RDIAG ACNORM dustem_mpfit_qrfac, fjac, ipvt, wa1, wa2, /pivot ;; Jacobian - jacobian matrix computed by mpfit_fdjac2 ;; Hh vects - house holder vectors from QR factorization & R matrix ;; Permutation - permutation vector for pivoting ;; RDIAG - diagonal elements of R matrix ;; ACNORM - norms of input Jacobian matrix before factoring ;; ===================== ;; On the first iteration if "diag" is unspecified, scale ;; according to the norms of the columns of the initial jacobian catch_msg = 'rescaling diagonal elements' if (iter EQ 1) then begin ;; Input: WA2 = root sum of squares of original Jacobian matrix ;; DIAG = user-requested diagonal (not documented!) ;; FACTOR = user-requested norm factor (not documented!) ;; Output: DIAG = Diagonal scaling values ;; XNORM = sum of squared scaled residuals ;; DELTA = rescaled XNORM if NOT keyword_set(rescale) OR (n_elements(diag) LT n) then begin diag = wa2 ;; Calculated from original Jacobian wh = where (diag EQ 0, ct) ;; Handle zero values if ct GT 0 then diag[wh] = one endif ;; On the first iteration, calculate the norm of the scaled x ;; and initialize the step bound delta wa3 = diag * x ;; WA3 is temp variable xnorm = dustem_mpfit_enorm(wa3) delta = factor*xnorm if delta EQ zero then delta = zero + factor endif ;; Form (q transpose)*fvec and store the first n components in qtf catch_msg = 'forming (q transpose)*fvec' wa4 = fvec for j=0L, n-1 do begin lj = ipvt[j] temp3 = fjac[j,lj] if temp3 NE 0 then begin fj = fjac[j:*,lj] wj = wa4[j:*] ;; *** optimization wa4(j:*) wa4[j] = wj - fj * total(fj*wj) / temp3 endif fjac[j,lj] = wa1[j] qtf[j] = wa4[j] endfor ;; From this point on, only the square matrix, consisting of the ;; triangle of R, is needed. fjac = fjac[0:n-1, 0:n-1] fjac = reform(fjac, n, n, /overwrite) fjac = fjac[*, ipvt] ;; Convert to permuted order fjac = reform(fjac, n, n, /overwrite) ;; Check for overflow. This should be a cheap test here since FJAC ;; has been reduced to a (small) square matrix, and the test is ;; O(N^2). wh = where(finite(fjac) EQ 0, ct) if ct GT 0 then goto, FAIL_OVERFLOW ;; Compute the norm of the scaled gradient catch_msg = 'computing the scaled gradient' gnorm = zero if fnorm NE 0 then begin for j=0L, n-1 do begin l = ipvt[j] if wa2[l] NE 0 then begin sum = total(fjac[0:j,j]*qtf[0:j])/fnorm gnorm = max([gnorm,abs(sum/wa2[l])]) endif endfor endif ;; Test for convergence of the gradient norm if gnorm LE gtol then info = 4 if info NE 0 then goto, TERMINATE if maxiter EQ 0 then begin info = 5 goto, TERMINATE endif ;; Rescale if necessary if NOT keyword_set(rescale) then $ diag = diag > wa2 ;; Beginning of the inner loop INNER_LOOP: ;; Determine the levenberg-marquardt parameter catch_msg = 'calculating LM parameter (DUSTEM_MPFIT_LMPAR)' par = dustem_mpfit_lmpar(fjac, ipvt, diag, qtf, delta, wa1, wa2, par=par) ;; Store the direction p and x+p. Calculate the norm of p wa1 = -wa1 if qanylim EQ 0 AND qminmax EQ 0 then begin ;; No parameter limits, so just move to new position WA2 alpha = one wa2 = x + wa1 endif else begin ;; Respect the limits. If a step were to go out of bounds, then ;; we should take a step in the same direction but shorter distance. ;; The step should take us right to the limit in that case. alpha = one if qanylim EQ 1 then begin ;; Do not allow any steps out of bounds catch_msg = 'checking for a step out of bounds' if nlpeg GT 0 then wa1[whlpeg] = wa1[whlpeg] > 0 if nupeg GT 0 then wa1[whupeg] = wa1[whupeg] < 0 dwa1 = abs(wa1) GT MACHEP0 whl = where(dwa1 AND qllim AND (x + wa1 LT llim), lct) if lct GT 0 then $ alpha = min([alpha, (llim[whl]-x[whl])/wa1[whl]]) whu = where(dwa1 AND qulim AND (x + wa1 GT ulim), uct) if uct GT 0 then $ alpha = min([alpha, (ulim[whu]-x[whu])/wa1[whu]]) endif ;; Obey any max step values. if qminmax EQ 1 then begin nwa1 = wa1 * alpha whmax = where(qmax AND maxstep GT 0, ct) if ct GT 0 then begin mrat = max(abs(nwa1[whmax])/abs(maxstep[whmax])) if mrat GT 1 then alpha = alpha / mrat endif endif ;; Scale the resulting vector wa1 = wa1 * alpha wa2 = x + wa1 ;; Adjust the final output values. If the step put us exactly ;; on a boundary, make sure we peg it there. sgnu = (ulim GE 0)*2d - 1d sgnl = (llim GE 0)*2d - 1d ;; Handles case of ;; ... nonzero *LIM ... ... zero *LIM ... ulim1 = ulim*(1-sgnu*MACHEP0) - (ulim EQ 0)*MACHEP0 llim1 = llim*(1+sgnl*MACHEP0) + (llim EQ 0)*MACHEP0 wh = where(qulim AND (wa2 GE ulim1), ct) if ct GT 0 then wa2[wh] = ulim[wh] wh = where(qllim AND (wa2 LE llim1), ct) if ct GT 0 then wa2[wh] = llim[wh] endelse wa3 = diag * wa1 pnorm = dustem_mpfit_enorm(wa3) ;; On the first iteration, adjust the initial step bound if iter EQ 1 then delta = min([delta,pnorm]) xnew[ifree] = wa2 if isext then goto, SAVE_STATE ;; Evaluate the function at x+p and calculate its norm mperr = 0 catch_msg = 'calling '+fcn wa4 = dustem_mpfit_call(fcn, xnew, _EXTRA=fcnargs) iflag = mperr if iflag LT 0 then begin errmsg = 'WARNING: premature termination by "'+fcn+'"' goto, TERMINATE endif RESUME_FIT: fnorm1 = dustem_mpfit_enorm(wa4) ;; Compute the scaled actual reduction catch_msg = 'computing convergence criteria' actred = -one if 0.1D * fnorm1 LT fnorm then actred = - (fnorm1/fnorm)^2 + 1. ;; Compute the scaled predicted reduction and the scaled directional ;; derivative for j = 0L, n-1 do begin wa3[j] = 0 wa3[0:j] = wa3[0:j] + fjac[0:j,j]*wa1[ipvt[j]] endfor ;; Remember, alpha is the fraction of the full LM step actually ;; taken temp1 = dustem_mpfit_enorm(alpha*wa3)/fnorm temp2 = (sqrt(alpha*par)*pnorm)/fnorm half = zero + 0.5 prered = temp1*temp1 + (temp2*temp2)/half dirder = -(temp1*temp1 + temp2*temp2) ;; Compute the ratio of the actual to the predicted reduction. ratio = zero tenth = zero + 0.1 if prered NE 0 then ratio = actred/prered ;; Update the step bound if ratio LE 0.25D then begin if actred GE 0 then temp = half $ else temp = half*dirder/(dirder + half*actred) if ((0.1D*fnorm1) GE fnorm) OR (temp LT 0.1D) then temp = tenth delta = temp*min([delta,pnorm/tenth]) par = par/temp endif else begin if (par EQ 0) OR (ratio GE 0.75) then begin delta = pnorm/half par = half*par endif endelse ;; Test for successful iteration if ratio GE 0.0001 then begin ;; Successful iteration. Update x, fvec, and their norms x = wa2 wa2 = diag * x fvec = wa4 xnorm = dustem_mpfit_enorm(wa2) fnorm = fnorm1 iter = iter + 1 ;DUSTEM SYSVARS ;#iterations !dustem_iter.act=iter ;#current parameter values (*!dustem_fit).current_param_values=ptr_new(x*(*(*!dustem_fit).param_init_values));So that parameters aren't just passed directly from MPFIT. ALSO allows to test with !dustem_current ; ;#Current parameter errors if n_elements(fnorm) GT 0 AND n_elements(fnorm1) GT 0 then begin fnorm = max([fnorm, fnorm1]) fnorm = fnorm^2. endif covar = !values.d_nan ;; (very carefully) set the covariance matrix COVAR if NOT keyword_set(nocovar) $ AND n_elements(n) GT 0 $ AND n_elements(fjac) GT 0 AND n_elements(ipvt) GT 0 then begin sz = size(fjac) if n GT 0 AND sz[0] GT 1 AND sz[1] GE n AND sz[2] GE n $ AND n_elements(ipvt) GE n then begin catch_msg = 'computing the covariance matrix' if n EQ 1 then $ cv = dustem_mpfit_covar(reform([fjac[0,0]],1,1), ipvt[0]) $ else $ cv = dustem_mpfit_covar(fjac[0:n-1,0:n-1], ipvt[0:n-1]) cv = reform(cv, n, n, /overwrite) nn = n_elements(xall) ;; Fill in actual covariance matrix, accounting for fixed ;; parameters. covar = replicate(zero, nn, nn) for i = 0L, n-1 do begin covar[ifree, ifree[i]] = cv[*,i] end ;; Compute errors in parameters catch_msg = 'computing parameter errors' i = lindgen(nn) perror = replicate(abs(covar[0])*0., nn) wh = where(covar[i,i] GE 0, ct) if ct GT 0 then $ perror[wh] = sqrt(covar[wh, wh]) (*!dustem_fit).current_param_errors=ptr_new(perror*(*(*!dustem_fit).param_init_values)) endif endif ; ;#current parameter errors ; fvec = mpfit_call(fcn, xnew, _EXTRA=fcnargs) ; ; fnorm = mpfit_enorm(fvec) ; ; if n_elements(fnorm) GT 0 AND n_elements(fnorm1) GT 0 then begin ; fnorm = max([fnorm, fnorm1]) ; fnorm = fnorm^2. ; endif ; stop ; covar = !values.d_nan ; ;; (very carefully) set the covariance matrix COVAR ;if info GT 0 AND NOT keyword_set(nocovar) AND n_elements(n) GT 0 AND n_elements(fjac) GT 0 AND n_elements(ipvt) GT 0 then begin ;if NOT keyword_set(nocovar) $ ; sz = size(fjac) ; if n GT 0 AND sz[0] GT 1 AND sz[1] GE n AND sz[2] GE n AND n_elements(ipvt) GE n then begin ; ; catch_msg = 'computing the covariance matrix' ; if n EQ 1 then $ ;cv = mpfit_covar(reform([fjac[0,0]],1,1), ipvt[0]) $ ;else $ ; cv = mpfit_covar(fjac[0:n-1,0:n-1], ipvt[0:n-1]) ; cv = reform(cv, n, n, /overwrite) ; nn = n_elements(xall) ; ; Fill in actual covariance matrix, accounting for fixed ; parameters. ; covar = replicate(zero, nn, nn) ; for i = 0L, n-1 do begin ; covar[ifree, ifree[i]] = cv[*,i] ; end ; Compute errors in parameters ; catch_msg = 'computing parameter errors' ; i = lindgen(nn) ; perror = replicate(abs(covar[0])*0., nn) ; wh = where(covar[i,i] GE 0, ct) ; ;stop ; if ct GT 0 then $ ; perror[wh] = sqrt(covar[wh, wh]) ; stop ; (*!dustem_fit).current_param_errors=ptr_new(perror*(*(*!dustem_fit).param_init_values)) ; endif; endif ;; Tests for convergence if (abs(actred) LE ftol) AND (prered LE ftol) $ AND (0.5D * ratio LE 1) then info = 1 if delta LE xtol*xnorm then info = 2 if (abs(actred) LE ftol) AND (prered LE ftol) $ AND (0.5D * ratio LE 1) AND (info EQ 2) then info = 3 if info NE 0 then goto, TERMINATE ;; Tests for termination and stringent tolerances if iter GE maxiter then info = 5 if (abs(actred) LE MACHEP0) AND (prered LE MACHEP0) $ AND (0.5*ratio LE 1) then info = 6 if delta LE MACHEP0*xnorm then info = 7 if gnorm LE MACHEP0 then info = 8 if info NE 0 then goto, TERMINATE ;; End of inner loop. Repeat if iteration unsuccessful if ratio LT 0.0001 then begin goto, INNER_LOOP endif ;; Check for over/underflow wh = where(finite(wa1) EQ 0 OR finite(wa2) EQ 0 OR finite(x) EQ 0, ct) if ct GT 0 OR finite(ratio) EQ 0 then begin FAIL_OVERFLOW: errmsg = ('ERROR: parameter or function value(s) have become '+$ 'infinite; check model function for over- '+$ 'and underflow') info = -16 goto, TERMINATE endif ;; End of outer loop. goto, OUTER_LOOP TERMINATE: catch_msg = 'in the termination phase' ;; Termination, either normal or user imposed. if iflag LT 0 then info = iflag iflag = 0 if n_elements(xnew) EQ 0 then goto, FINAL_RETURN if nfree EQ 0 then xnew = xall else xnew[ifree] = x if n_elements(qanytied) GT 0 then if qanytied then dustem_mpfit_tie, xnew, ptied dof = n_elements(fvec) - nfree ;; Call the ITERPROC at the end of the fit, if the fit status is ;; okay. Don't call it if the fit failed for some reason. if info GT 0 then begin mperr = 0 xnew0 = xnew call_procedure, iterproc, fcn, xnew, iter, fnorm^2, $ FUNCTARGS=fcnargs, parinfo=parinfo, quiet=quiet, $ dof=dof, _EXTRA=iterargs iflag = mperr if iflag LT 0 then begin errmsg = 'WARNING: premature termination by "'+iterproc+'"' endif else begin ;; If parameters were changed (grrr..) then re-tie if max(abs(xnew0-xnew)) GT 0 then begin if qanytied then dustem_mpfit_tie, xnew, ptied x = xnew[ifree] endif endelse endif ;; Initialize the number of parameters pegged at a hard limit value npegged = 0L if n_elements(qanylim) GT 0 then if qanylim then begin wh = where((qulim AND (x EQ ulim)) OR $ (qllim AND (x EQ llim)), npegged) endif ;; Calculate final function value (FNORM) and residuals (FVEC) if isext EQ 0 AND nprint GT 0 AND info GT 0 then begin catch_msg = 'calling '+fcn fvec = dustem_mpfit_call(fcn, xnew, _EXTRA=fcnargs) catch_msg = 'in the termination phase' fnorm = dustem_mpfit_enorm(fvec) endif if n_elements(fnorm) GT 0 AND n_elements(fnorm1) GT 0 then begin fnorm = max([fnorm, fnorm1]) fnorm = fnorm^2. endif covar = !values.d_nan ;; (very carefully) set the covariance matrix COVAR if info GT 0 AND NOT keyword_set(nocovar) $ AND n_elements(n) GT 0 $ AND n_elements(fjac) GT 0 AND n_elements(ipvt) GT 0 then begin sz = size(fjac) if n GT 0 AND sz[0] GT 1 AND sz[1] GE n AND sz[2] GE n $ AND n_elements(ipvt) GE n then begin catch_msg = 'computing the covariance matrix' if n EQ 1 then $ cv = dustem_mpfit_covar(reform([fjac[0,0]],1,1), ipvt[0]) $ else $ cv = dustem_mpfit_covar(fjac[0:n-1,0:n-1], ipvt[0:n-1]) cv = reform(cv, n, n, /overwrite) nn = n_elements(xall) ;; Fill in actual covariance matrix, accounting for fixed ;; parameters. covar = replicate(zero, nn, nn) for i = 0L, n-1 do begin covar[ifree, ifree[i]] = cv[*,i] end ;; Compute errors in parameters catch_msg = 'computing parameter errors' i = lindgen(nn) perror = replicate(abs(covar[0])*0., nn) wh = where(covar[i,i] GE 0, ct) if ct GT 0 then $ perror[wh] = sqrt(covar[wh, wh]) (*!dustem_fit).current_param_errors=ptr_new(perror*(*(*!dustem_fit).param_init_values)) endif endif ; catch_msg = 'returning the result' ; profvals.mpfit = profvals.mpfit + (systime(1) - prof_start) FINAL_RETURN: mpfit_fencepost_active = 0 nfev = mpconfig.nfev if n_elements(xnew) EQ 0 then return, !values.d_nan return, xnew ;; ------------------------------------------------------------------ ;; Alternate ending if the user supplies the function and gradients ;; externally ;; ------------------------------------------------------------------ SAVE_STATE: catch_msg = 'saving DUSTEM_MPFIT state' ;; Names of variables to save varlist = ['alpha', 'delta', 'diag', 'dwarf', 'factor', 'fnorm', $ 'fjac', 'gnorm', 'nfree', 'ifree', 'ipvt', 'iter', $ 'm', 'n', 'machvals', 'machep0', 'npegged', $ 'whlpeg', 'whupeg', 'nlpeg', 'nupeg', $ 'mpconfig', 'par', 'pnorm', 'qtf', $ 'wa1', 'wa2', 'wa3', 'xnorm', 'x', 'xnew'] cmd = '' ;; Construct an expression that will save them for i = 0L, n_elements(varlist)-1 do begin ival = 0 dummy = execute('ival = n_elements('+varlist[i]+')') if ival GT 0 then begin cmd = cmd + ',' + varlist[i]+':'+varlist[i] endif endfor cmd = 'state = create_struct({'+strmid(cmd,1)+'})' state = 0 if execute(cmd) NE 1 then $ message, 'ERROR: could not save DUSTEM_MPFIT state' ;; Set STATUS keyword to prepare for next iteration, and reset init ;; so we do not init the next time info = 9 extinit = 0 return, xnew end