1;+ 2; NAME: 3; MPCURVEFIT 4; 5; AUTHOR: 6; Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770 7; craigm@lheamail.gsfc.nasa.gov 8; UPDATED VERSIONs can be found on my WEB PAGE: 9; http://cow.physics.wisc.edu/~craigm/idl/idl.html 10; 11; PURPOSE: 12; Perform Levenberg-Marquardt least-squares fit (replaces CURVEFIT) 13; 14; MAJOR TOPICS: 15; Curve and Surface Fitting 16; 17; CALLING SEQUENCE: 18; YFIT = MPCURVEFIT(X, Y, WEIGHTS, P, [SIGMA,] FUNCTION_NAME=FUNC, 19; ITER=iter, ITMAX=itmax, 20; CHISQ=chisq, NFREE=nfree, DOF=dof, 21; NFEV=nfev, COVAR=covar, [/NOCOVAR, ] [/NODERIVATIVE, ] 22; FUNCTARGS=functargs, PARINFO=parinfo, 23; FTOL=ftol, XTOL=xtol, GTOL=gtol, TOL=tol, 24; ITERPROC=iterproc, ITERARGS=iterargs, 25; NPRINT=nprint, QUIET=quiet, 26; ERRMSG=errmsg, STATUS=status) 27; 28; DESCRIPTION: 29; 30; MPCURVEFIT fits a user-supplied model -- in the form of an IDL 31; function -- to a set of user-supplied data. MPCURVEFIT calls 32; MPFIT, the MINPACK-1 least-squares minimizer, to do the main 33; work. 34; 35; Given the data and their uncertainties, MPCURVEFIT finds the best 36; set of model parameters which match the data (in a least-squares 37; sense) and returns them in the parameter P. 38; 39; MPCURVEFIT returns the best fit function. 40; 41; The user must supply the following items: 42; - An array of independent variable values ("X"). 43; - An array of "measured" *dependent* variable values ("Y"). 44; - An array of weighting values ("WEIGHTS"). 45; - The name of an IDL function which computes Y given X ("FUNC"). 46; - Starting guesses for all of the parameters ("P"). 47; 48; There are very few restrictions placed on X, Y or FUNCT. Simply 49; put, FUNCT must map the "X" values into "Y" values given the 50; model parameters. The "X" values may represent any independent 51; variable (not just Cartesian X), and indeed may be multidimensional 52; themselves. For example, in the application of image fitting, X 53; may be a 2xN array of image positions. 54; 55; MPCURVEFIT carefully avoids passing large arrays where possible to 56; improve performance. 57; 58; See below for an example of usage. 59; 60; USER FUNCTION 61; 62; The user must define a function which returns the model value. For 63; applications which use finite-difference derivatives -- the default 64; -- the user function should be declared in the following way: 65; 66; ; MYFUNCT - example user function 67; ; X - input independent variable (vector same size as data) 68; ; P - input parameter values (N-element array) 69; ; YMOD - upon return, user function values 70; ; DP - upon return, the user function must return 71; ; an ARRAY(M,N) of derivatives in this parameter 72; ; 73; PRO MYFUNCT, x, p, ymod, dp 74; ymod = F(x, p) ;; Model function 75; 76; if n_params() GE 4 then begin 77; ; Create derivative and compute derivative array 78; dp = make_array(n_elements(x), n_elements(p), value=x[0]*0) 79; 80; ; Compute derivative if requested by caller 81; for i = 0, n_elements(p)-1 do dp(*,i) = FGRAD(x, p, i) 82; endif 83; END 84; 85; where FGRAD(x, p, i) is a model function which computes the 86; derivative of the model F(x,p) with respect to parameter P(i) at X. 87; The returned array YMOD must have the same dimensions and type as 88; the "measured" Y values. The returned array DP[i,j] is the 89; derivative of the ith function value with respect to the jth 90; parameter. 91; 92; User functions may also indicate a fatal error condition 93; using the ERROR_CODE common block variable, as described 94; below under the MPFIT_ERROR common block definition. 95; 96; If NODERIVATIVE=1, then MPCURVEFIT will never request explicit 97; derivatives from the user function, and instead will user numerical 98; estimates (i.e. by calling the user function multiple times). 99; 100; CONSTRAINING PARAMETER VALUES WITH THE PARINFO KEYWORD 101; 102; The behavior of MPFIT can be modified with respect to each 103; parameter to be fitted. A parameter value can be fixed; simple 104; boundary constraints can be imposed; limitations on the parameter 105; changes can be imposed; properties of the automatic derivative can 106; be modified; and parameters can be tied to one another. 107; 108; These properties are governed by the PARINFO structure, which is 109; passed as a keyword parameter to MPFIT. 110; 111; PARINFO should be an array of structures, one for each parameter. 112; Each parameter is associated with one element of the array, in 113; numerical order. The structure can have the following entries 114; (none are required): 115; 116; .VALUE - the starting parameter value (but see the START_PARAMS 117; parameter for more information). 118; 119; .FIXED - a boolean value, whether the parameter is to be held 120; fixed or not. Fixed parameters are not varied by 121; MPFIT, but are passed on to MYFUNCT for evaluation. 122; 123; .LIMITED - a two-element boolean array. If the first/second 124; element is set, then the parameter is bounded on the 125; lower/upper side. A parameter can be bounded on both 126; sides. Both LIMITED and LIMITS must be given 127; together. 128; 129; .LIMITS - a two-element float or double array. Gives the 130; parameter limits on the lower and upper sides, 131; respectively. Zero, one or two of these values can be 132; set, depending on the values of LIMITED. Both LIMITED 133; and LIMITS must be given together. 134; 135; .PARNAME - a string, giving the name of the parameter. The 136; fitting code of MPFIT does not use this tag in any 137; way. However, the default ITERPROC will print the 138; parameter name if available. 139; 140; .STEP - the step size to be used in calculating the numerical 141; derivatives. If set to zero, then the step size is 142; computed automatically. Ignored when AUTODERIVATIVE=0. 143; This value is superceded by the RELSTEP value. 144; 145; .RELSTEP - the *relative* step size to be used in calculating 146; the numerical derivatives. This number is the 147; fractional size of the step, compared to the 148; parameter value. This value supercedes the STEP 149; setting. If the parameter is zero, then a default 150; step size is chosen. 151; 152; .MPSIDE - the sidedness of the finite difference when computing 153; numerical derivatives. This field can take four 154; values: 155; 156; 0 - one-sided derivative computed automatically 157; 1 - one-sided derivative (f(x+h) - f(x) )/h 158; -1 - one-sided derivative (f(x) - f(x-h))/h 159; 2 - two-sided derivative (f(x+h) - f(x-h))/(2*h) 160; 161; Where H is the STEP parameter described above. The 162; "automatic" one-sided derivative method will chose a 163; direction for the finite difference which does not 164; violate any constraints. The other methods do not 165; perform this check. The two-sided method is in 166; principle more precise, but requires twice as many 167; function evaluations. Default: 0. 168; 169; .MPMAXSTEP - the maximum change to be made in the parameter 170; value. During the fitting process, the parameter 171; will never be changed by more than this value in 172; one iteration. 173; 174; A value of 0 indicates no maximum. Default: 0. 175; 176; .TIED - a string expression which "ties" the parameter to other 177; free or fixed parameters. Any expression involving 178; constants and the parameter array P are permitted. 179; Example: if parameter 2 is always to be twice parameter 180; 1 then use the following: parinfo(2).tied = '2 * P(1)'. 181; Since they are totally constrained, tied parameters are 182; considered to be fixed; no errors are computed for them. 183; [ NOTE: the PARNAME can't be used in expressions. ] 184; 185; .MPPRINT - if set to 1, then the default ITERPROC will print the 186; parameter value. If set to 0, the parameter value 187; will not be printed. This tag can be used to 188; selectively print only a few parameter values out of 189; many. Default: 1 (all parameters printed) 190; 191; 192; Future modifications to the PARINFO structure, if any, will involve 193; adding structure tags beginning with the two letters "MP". 194; Therefore programmers are urged to avoid using tags starting with 195; the same letters; otherwise they are free to include their own 196; fields within the PARINFO structure, and they will be ignored. 197; 198; PARINFO Example: 199; parinfo = replicate({value:0.D, fixed:0, limited:[0,0], $ 200; limits:[0.D,0]}, 5) 201; parinfo(0).fixed = 1 202; parinfo(4).limited(0) = 1 203; parinfo(4).limits(0) = 50.D 204; parinfo(*).value = [5.7D, 2.2, 500., 1.5, 2000.] 205; 206; A total of 5 parameters, with starting values of 5.7, 207; 2.2, 500, 1.5, and 2000 are given. The first parameter 208; is fixed at a value of 5.7, and the last parameter is 209; constrained to be above 50. 210; 211; INPUTS: 212; X - Array of independent variable values. 213; 214; Y - Array of "measured" dependent variable values. Y should have 215; the same data type as X. The function FUNCT should map 216; X->Y. 217; 218; WEIGHTS - Array of weights to be used in calculating the 219; chi-squared value. If WEIGHTS is specified then the ERR 220; parameter is ignored. The chi-squared value is computed 221; as follows: 222; 223; CHISQ = TOTAL( (Y-FUNCT(X,P))^2 * ABS(WEIGHTS) ) 224; 225; Here are common values of WEIGHTS: 226; 227; 1D/ERR^2 - Normal weighting (ERR is the measurement error) 228; 1D/Y - Poisson weighting (counting statistics) 229; 1D - Unweighted 230; 231; P - An array of starting values for each of the parameters of the 232; model. The number of parameters should be fewer than the 233; number of measurements. Also, the parameters should have the 234; same data type as the measurements (double is preferred). 235; 236; Upon successful completion the new parameter values are 237; returned in P. 238; 239; If both START_PARAMS and PARINFO are passed, then the starting 240; *value* is taken from START_PARAMS, but the *constraints* are 241; taken from PARINFO. 242; 243; SIGMA - The formal 1-sigma errors in each parameter, computed from 244; the covariance matrix. If a parameter is held fixed, or 245; if it touches a boundary, then the error is reported as 246; zero. 247; 248; If the fit is unweighted (i.e. no errors were given, or 249; the weights were uniformly set to unity), then SIGMA will 250; probably not represent the true parameter uncertainties. 251; 252; *If* you can assume that the true reduced chi-squared 253; value is unity -- meaning that the fit is implicitly 254; assumed to be of good quality -- then the estimated 255; parameter uncertainties can be computed by scaling SIGMA 256; by the measured chi-squared value. 257; 258; DOF = N_ELEMENTS(X) - N_ELEMENTS(P) ; deg of freedom 259; CSIGMA = SIGMA * SQRT(CHISQ / DOF) ; scaled uncertainties 260; 261; RETURNS: 262; 263; Returns the array containing the best-fitting function. 264; 265; KEYWORD PARAMETERS: 266; 267; CHISQ - the value of the summed, squared, weighted residuals for 268; the returned parameter values, i.e. the chi-square value. 269; 270; COVAR - the covariance matrix for the set of parameters returned 271; by MPFIT. The matrix is NxN where N is the number of 272; parameters. The square root of the diagonal elements 273; gives the formal 1-sigma statistical errors on the 274; parameters IF errors were treated "properly" in MYFUNC. 275; Parameter errors are also returned in PERROR. 276; 277; To compute the correlation matrix, PCOR, use this: 278; IDL> PCOR = COV * 0 279; IDL> FOR i = 0, n-1 DO FOR j = 0, n-1 DO $ 280; PCOR(i,j) = COV(i,j)/sqrt(COV(i,i)*COV(j,j)) 281; 282; If NOCOVAR is set or MPFIT terminated abnormally, then 283; COVAR is set to a scalar with value !VALUES.D_NAN. 284; 285; DOF - number of degrees of freedom, computed as 286; DOF = N_ELEMENTS(DEVIATES) - NFREE 287; Note that this doesn't account for pegged parameters (see 288; NPEGGED). 289; 290; ERRMSG - a string error or warning message is returned. 291; 292; FTOL - a nonnegative input variable. Termination occurs when both 293; the actual and predicted relative reductions in the sum of 294; squares are at most FTOL (and STATUS is accordingly set to 295; 1 or 3). Therefore, FTOL measures the relative error 296; desired in the sum of squares. Default: 1D-10 297; 298; FUNCTION_NAME - a scalar string containing the name of an IDL 299; procedure to compute the user model values, as 300; described above in the "USER MODEL" section. 301; 302; FUNCTARGS - A structure which contains the parameters to be passed 303; to the user-supplied function specified by FUNCT via 304; the _EXTRA mechanism. This is the way you can pass 305; additional data to your user-supplied function without 306; using common blocks. 307; 308; By default, no extra parameters are passed to the 309; user-supplied function. 310; 311; GTOL - a nonnegative input variable. Termination occurs when the 312; cosine of the angle between fvec and any column of the 313; jacobian is at most GTOL in absolute value (and STATUS is 314; accordingly set to 4). Therefore, GTOL measures the 315; orthogonality desired between the function vector and the 316; columns of the jacobian. Default: 1D-10 317; 318; ITER - the number of iterations completed. 319; 320; ITERARGS - The keyword arguments to be passed to ITERPROC via the 321; _EXTRA mechanism. This should be a structure, and is 322; similar in operation to FUNCTARGS. 323; Default: no arguments are passed. 324; 325; ITERPROC - The name of a procedure to be called upon each NPRINT 326; iteration of the MPFIT routine. It should be declared 327; in the following way: 328; 329; PRO ITERPROC, FUNCT, p, iter, fnorm, FUNCTARGS=fcnargs, $ 330; PARINFO=parinfo, QUIET=quiet, ... 331; ; perform custom iteration update 332; END 333; 334; ITERPROC must either accept all three keyword 335; parameters (FUNCTARGS, PARINFO and QUIET), or at least 336; accept them via the _EXTRA keyword. 337; 338; FUNCT is the user-supplied function to be minimized, 339; P is the current set of model parameters, ITER is the 340; iteration number, and FUNCTARGS are the arguments to be 341; passed to FUNCT. FNORM should be the 342; chi-squared value. QUIET is set when no textual output 343; should be printed. See below for documentation of 344; PARINFO. 345; 346; In implementation, ITERPROC can perform updates to the 347; terminal or graphical user interface, to provide 348; feedback while the fit proceeds. If the fit is to be 349; stopped for any reason, then ITERPROC should set the 350; common block variable ERROR_CODE to negative value (see 351; MPFIT_ERROR common block below). In principle, 352; ITERPROC should probably not modify the parameter 353; values, because it may interfere with the algorithm's 354; stability. In practice it is allowed. 355; 356; Default: an internal routine is used to print the 357; parameter values. 358; 359; ITMAX - The maximum number of iterations to perform. If the 360; number is exceeded, then the STATUS value is set to 5 361; and MPFIT returns. 362; Default: 200 iterations 363; 364; NFEV - the number of FUNCT function evaluations performed. 365; 366; NFREE - the number of free parameters in the fit. This includes 367; parameters which are not FIXED and not TIED, but it does 368; include parameters which are pegged at LIMITS. 369; 370; NOCOVAR - set this keyword to prevent the calculation of the 371; covariance matrix before returning (see COVAR) 372; 373; NODERIVATIVE - if set, then the user function will not be queried 374; for analytical derivatives, and instead the 375; derivatives will be computed by finite differences 376; (and according to the PARINFO derivative settings; 377; see above for a description). 378; 379; NPRINT - The frequency with which ITERPROC is called. A value of 380; 1 indicates that ITERPROC is called with every iteration, 381; while 2 indicates every other iteration, etc. Note that 382; several Levenberg-Marquardt attempts can be made in a 383; single iteration. 384; Default value: 1 385; 386; PARINFO - Provides a mechanism for more sophisticated constraints 387; to be placed on parameter values. When PARINFO is not 388; passed, then it is assumed that all parameters are free 389; and unconstrained. Values in PARINFO are never 390; modified during a call to MPFIT. 391; 392; See description above for the structure of PARINFO. 393; 394; Default value: all parameters are free and unconstrained. 395; 396; QUIET - set this keyword when no textual output should be printed 397; by MPFIT 398; 399; STATUS - an integer status code is returned. All values greater 400; than zero can represent success (however STATUS EQ 5 may 401; indicate failure to converge). It can have one of the 402; following values: 403; 404; 0 improper input parameters. 405; 406; 1 both actual and predicted relative reductions 407; in the sum of squares are at most FTOL. 408; 409; 2 relative error between two consecutive iterates 410; is at most XTOL 411; 412; 3 conditions for STATUS = 1 and STATUS = 2 both hold. 413; 414; 4 the cosine of the angle between fvec and any 415; column of the jacobian is at most GTOL in 416; absolute value. 417; 418; 5 the maximum number of iterations has been reached 419; 420; 6 FTOL is too small. no further reduction in 421; the sum of squares is possible. 422; 423; 7 XTOL is too small. no further improvement in 424; the approximate solution x is possible. 425; 426; 8 GTOL is too small. fvec is orthogonal to the 427; columns of the jacobian to machine precision. 428; 429; TOL - synonym for FTOL. Use FTOL instead. 430; 431; XTOL - a nonnegative input variable. Termination occurs when the 432; relative error between two consecutive iterates is at most 433; XTOL (and STATUS is accordingly set to 2 or 3). Therefore, 434; XTOL measures the relative error desired in the approximate 435; solution. Default: 1D-10 436; 437; YERROR - upon return, the root-mean-square variance of the 438; residuals. 439; 440; 441; EXAMPLE: 442; 443; ; First, generate some synthetic data 444; npts = 200 445; x = dindgen(npts) * 0.1 - 10. ; Independent variable 446; yi = gauss1(x, [2.2D, 1.4, 3000.]) ; "Ideal" Y variable 447; y = yi + randomn(seed, npts) * sqrt(1000. + yi); Measured, w/ noise 448; sy = sqrt(1000.D + y) ; Poisson errors 449; 450; ; Now fit a Gaussian to see how well we can recover 451; p0 = [1.D, 1., 1000.] ; Initial guess 452; yfit = mpcurvefit(x, y, 1/sy^2, p0, $ ; Fit a function 453; FUNCTION_NAME='GAUSS1P',/autoderivative) 454; print, p 455; 456; Generates a synthetic data set with a Gaussian peak, and Poisson 457; statistical uncertainty. Then the same function is fitted to the 458; data to see how close we can get. GAUSS1 and GAUSS1P are 459; available from the same web page. 460; 461; 462; COMMON BLOCKS: 463; 464; COMMON MPFIT_ERROR, ERROR_CODE 465; 466; User routines may stop the fitting process at any time by 467; setting an error condition. This condition may be set in either 468; the user's model computation routine (MYFUNCT), or in the 469; iteration procedure (ITERPROC). 470; 471; To stop the fitting, the above common block must be declared, 472; and ERROR_CODE must be set to a negative number. After the user 473; procedure or function returns, MPFIT checks the value of this 474; common block variable and exits immediately if the error 475; condition has been set. By default the value of ERROR_CODE is 476; zero, indicating a successful function/procedure call. 477; 478; REFERENCES: 479; 480; MINPACK-1, Jorge More', available from netlib (www.netlib.org). 481; "Optimization Software Guide," Jorge More' and Stephen Wright, 482; SIAM, *Frontiers in Applied Mathematics*, Number 14. 483; 484; MODIFICATION HISTORY: 485; Translated from MPFITFUN, 25 Sep 1999, CM 486; Alphabetized documented keywords, 02 Oct 1999, CM 487; Added QUERY keyword and query checking of MPFIT, 29 Oct 1999, CM 488; Check to be sure that X and Y are present, 02 Nov 1999, CM 489; Documented SIGMA for unweighted fits, 03 Nov 1999, CM 490; Changed to ERROR_CODE for error condition, 28 Jan 2000, CM 491; Copying permission terms have been liberalized, 26 Mar 2000, CM 492; Propagated improvements from MPFIT, 17 Dec 2000, CM 493; Corrected behavior of NODERIVATIVE, 13 May 2002, CM 494; Documented RELSTEP field of PARINFO (!!), CM, 25 Oct 2002 495; Make more consistent with comparable IDL routines, 30 Jun 2003, CM 496; Minor documentation adjustment, 03 Feb 2004, CM 497; Fix error in documentation, 26 Aug 2005, CM 498; Convert to IDL 5 array syntax (!), 16 Jul 2006, CM 499; Move STRICTARR compile option inside each function/procedure, 9 Oct 2006 500; Fix bug in handling of explicit derivatives with errors/weights 501; (the weights were not being applied), CM, 2012-07-22 502; Add more documentation on calling interface for user function and 503; parameter derivatives, CM, 2012-07-22 504; Better documentation for STATUS, CM, 2016-04-29 505; 506; $Id: mpcurvefit.pro,v 1.12 2016/05/19 16:08:49 cmarkwar Exp $ 507;- 508; Copyright (C) 1997-2000, 2002, 2003, 2004, 2005, 2012, 2016 Craig Markwardt 509; This software is provided as is without any warranty whatsoever. 510; Permission to use, copy, modify, and distribute modified or 511; unmodified copies is granted, provided this copyright and disclaimer 512; are included unchanged. 513;- 514 515FORWARD_FUNCTION mpcurvefit_eval, mpcurvefit, mpfit 516 517; This is the call-back function for MPFIT. It evaluates the 518; function, subtracts the data, and returns the residuals. 519function mpcurvefit_eval, p, dp, _EXTRA=extra 520 521 COMPILE_OPT strictarr 522 common mpcurvefit_common, fcn, x, y, wts, f, fcnargs 523 524 ;; The function is evaluated here. There are four choices, 525 ;; depending on whether (a) FUNCTARGS was passed to MPCURVEFIT, which 526 ;; is passed to this function as "hf"; or (b) the derivative 527 ;; parameter "dp" is passed, meaning that derivatives should be 528 ;; calculated analytically by the function itself. 529 if n_elements(fcnargs) GT 0 then begin 530 if n_params() GT 1 then call_procedure, fcn, x, p, f, dp,_EXTRA=fcnargs $ 531 else call_procedure, fcn, x, p, f, _EXTRA=fcnargs 532 endif else begin 533 if n_params() GT 1 then call_procedure, fcn, x, p, f, dp $ 534 else call_procedure, fcn, x, p, f 535 endelse 536 537 ;; Compute the deviates, applying the weights 538 result = (y-f)*wts 539 540 ;; Apply weights to derivative quantities 541 if n_params() GT 1 then begin 542 np = n_elements(p) 543 nf = n_elements(f) 544 for j = 0L, np-1 do dp[j*nf] = dp[j*nf:j*nf+nf-1] * wts 545 endif 546 547 ;; Make sure the returned result is one-dimensional. 548 result = reform(result, n_elements(result), /overwrite) 549 return, result 550 551end 552 553function mpcurvefit, x, y, wts, p, perror, function_name=fcn, $ 554 iter=iter, itmax=maxiter, $ 555 chisq=bestnorm, nfree=nfree, dof=dof, $ 556 nfev=nfev, covar=covar, nocovar=nocovar, yerror=yerror, $ 557 noderivative=noderivative, tol=tol, ftol=ftol, $ 558 FUNCTARGS=fa, parinfo=parinfo, $ 559 errmsg=errmsg, STATUS=status, QUIET=quiet, $ 560 query=query, _EXTRA=extra 561 562 COMPILE_OPT strictarr 563 status = 0L 564 errmsg = '' 565 566 ;; Detect MPFIT and crash if it was not found 567 catch, catcherror 568 if catcherror NE 0 then begin 569 MPFIT_NOTFOUND: 570 catch, /cancel 571 message, 'ERROR: the required function MPFIT must be in your IDL path', /info 572 return, !values.d_nan 573 endif 574 if mpfit(/query) NE 1 then goto, MPFIT_NOTFOUND 575 catch, /cancel 576 if keyword_set(query) then return, 1 577 578 if n_params() EQ 0 then begin 579 message, "USAGE: YFIT = MPCURVEFIT(X, Y, WTS, P, DP)", /info 580 return, !values.d_nan 581 endif 582 if n_elements(x) EQ 0 OR n_elements(y) EQ 0 then begin 583 message, 'ERROR: X and Y must be defined', /info 584 return, !values.d_nan 585 endif 586 if n_elements(fcn) EQ 0 then fcn = 'funct' 587 if n_elements(noderivative) EQ 0 then noderivative = 0 588 589 common mpcurvefit_common, fc, xc, yc, wc, mc, ac 590 fc = fcn & xc = x & yc = y & wc = sqrt(abs(wts)) & mc = 0L 591 ac = 0 & dummy = size(temporary(ac)) 592 if n_elements(fa) GT 0 then ac = fa 593 594 if n_elements(tol) GT 0 then ftol = tol 595 596 result = mpfit('mpcurvefit_eval', p, maxiter=maxiter, $ 597 autoderivative=noderivative, ftol=ftol, $ 598 parinfo=parinfo, STATUS=status, nfev=nfev, BESTNORM=bestnorm,$ 599 covar=covar, perror=perror, niter=iter, nfree=nfree, dof=dof,$ 600 ERRMSG=errmsg, quiet=quiet, _EXTRA=extra) 601 602 ;; Retrieve the fit value 603 yfit = temporary(mc) 604 ;; Now do some clean-up 605 xc = 0 & yc = 0 & wc = 0 & mc = 0 & ac = 0 606 607 if NOT keyword_set(quiet) AND errmsg NE '' then $ 608 message, errmsg, /info $ 609 else $ 610 p = result 611 612 yerror = p[0]*0 613 if n_elements(dof) GT 0 AND dof[0] GT 0 then begin 614 yerror[0] = sqrt( total( (y-yfit)^2 ) / dof[0] ) 615 endif 616 617 return, yfit 618end 619