;+ ; NAME: ; TNMIN ; ; 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: ; Performs function minimization (Truncated-Newton Method) ; ; MAJOR TOPICS: ; Optimization and Minimization ; ; CALLING SEQUENCE: ; parms = TNMIN(MYFUNCT, X, FUNCTARGS=fcnargs, NFEV=nfev, ; MAXITER=maxiter, ERRMSG=errmsg, NPRINT=nprint, ; QUIET=quiet, XTOL=xtol, STATUS=status, ; FGUESS=fguess, PARINFO=parinfo, BESTMIN=bestmin, ; ITERPROC=iterproc, ITERARGS=iterargs, niter=niter) ; ; DESCRIPTION: ; ; TNMIN uses the Truncated-Newton method to minimize an arbitrary IDL ; function with respect to a given set of free parameters. The ; user-supplied function must compute the gradient with respect to ; each parameter. TNMIN is based on TN.F (TNBC) by Stephen Nash. ; ; If you want to solve a least-squares problem, to perform *curve* ; *fitting*, then you will probably want to use the routines MPFIT, ; MPFITFUN and MPFITEXPR. Those routines are specifically optimized ; for the least-squares problem. TNMIN is suitable for constrained ; and unconstrained optimization problems with a medium number of ; variables. Function *maximization* can be performed using the ; MAXIMIZE keyword. ; ; TNMIN is similar to MPFIT in that it allows parameters to be fixed, ; simple bounding limits to be placed on parameter values, and ; parameters to be tied to other parameters. One major difference ; between MPFIT and TNMIN is that TNMIN does not compute derivatives ; automatically by default. See PARINFO and AUTODERIVATIVE below for ; more discussion and examples. ; ; USER FUNCTION ; ; The user must define an IDL function which returns the desired ; value as a single scalar. The IDL function must accept a list of ; numerical parameters, P. Additionally, keyword parameters may be ; used to pass more data or information to the user function, via the ; FUNCTARGS keyword. ; ; The user function should be declared in the following way: ; ; FUNCTION MYFUNCT, p, dp [, keywords permitted ] ; ; Parameter values are passed in "p" ; f = .... ; Compute function value ; dp = .... ; Compute partial derivatives (optional) ; return, f ; END ; ; The function *must* accept at least one argument, the parameter ; list P. The vector P is implicitly assumed to be a one-dimensional ; array. Users may pass additional information to the function by ; composing and _EXTRA structure and passing it in the FUNCTARGS ; keyword. ; ; User functions may also indicate a fatal error condition using the ; ERROR_CODE common block variable, as described below under the ; TNMIN_ERROR common block definition (by setting ERROR_CODE to a ; number between -15 and -1). ; ; EXPLICIT vs. NUMERICAL DERIVATIVES ; ; By default, the user must compute the function gradient components ; explicitly using AUTODERIVATIVE=0. As explained below, numerical ; derivatives can also be calculated using AUTODERIVATIVE=1. ; ; For explicit derivatives, the user should call TNMIN using the ; default keyword value AUTODERIVATIVE=0. [ This is different ; behavior from MPFIT, where AUTODERIVATIVE=1 is the default. ] The ; IDL user routine should compute the gradient of the function as a ; one-dimensional array of values, one for each of the parameters. ; They are passed back to TNMIN via "dp" as shown above. ; ; 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. ; ; For numerical derivatives, a finite differencing approximation is ; used to estimate the gradient values. Users can activate this ; feature by passing the keyword AUTODERIVATIVE=1. Fine control over ; this behavior can be achieved using the STEP, RELSTEP and TNSIDE ; fields of the PARINFO structure. ; ; When finite differencing is used for computing derivatives (ie, ; when AUTODERIVATIVE=1), the parameter DP is not passed. Therefore ; functions can use N_PARAMS() to indicate whether they must compute ; the derivatives or not. However there is no penalty (other than ; computation time) for computing the gradient values and users may ; switch between AUTODERIVATIVE=0 or =1 in order to test both ; scenarios. ; ; CONSTRAINING PARAMETER VALUES WITH THE PARINFO KEYWORD ; ; The behavior of TNMIN can be modified with respect to each ; parameter to be optimized. 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 TNMIN. ; ; 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 ; TNMIN, 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 TNMIN does not use this tag in any ; way. ; ; .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. ; ; .TNSIDE - the sidedness of the finite difference when computing ; numerical derivatives. This field can take four ; values: ; ; 0 - one-sided derivative computed automatically ; 1 - one-sided derivative (f(x+h) - f(x) )/h ; -1 - one-sided derivative (f(x) - f(x-h))/h ; 2 - two-sided derivative (f(x+h) - f(x-h))/(2*h) ; ; Where H is the STEP parameter described above. The ; "automatic" one-sided derivative method will chose a ; direction for the finite difference which does not ; violate any constraints. The other methods do not ; perform this check. The two-sided method is in ; principle more precise, but requires twice as many ; function evaluations. Default: 0. ; ; .TIED - a string expression which "ties" the parameter to other ; free or fixed parameters. Any expression involving ; constants and the parameter array P are permitted. ; Example: if parameter 2 is always to be twice parameter ; 1 then use the following: parinfo(2).tied = '2 * P(1)'. ; Since they are totally constrained, tied parameters are ; considered to be fixed; no errors are computed for them. ; [ NOTE: the PARNAME can't be used in expressions. ] ; ; Future modifications to the PARINFO structure, if any, will involve ; adding structure tags beginning with the two letters "MP" or "TN". ; Therefore programmers are urged to avoid using tags starting with ; these two combinations of letters; otherwise they are free to ; include their own fields within the PARINFO structure, and they ; will be ignored. ; ; 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. ; ; ; INPUTS: ; ; MYFUNCT - a string variable containing the name of the function to ; be minimized (see USER FUNCTION above). The IDL routine ; should return the value of the function and optionally ; its gradients. ; ; X - An array of starting values for each of the parameters of the ; model. ; ; This parameter is optional if the PARINFO keyword is used. ; See above. The PARINFO keyword provides a mechanism to fix or ; constrain individual parameters. If both X and PARINFO are ; passed, then the starting *value* is taken from X, but the ; *constraints* are taken from PARINFO. ; ; ; RETURNS: ; ; Returns the array of best-fit parameters. ; ; ; 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: 0 (explicit derivatives required) ; ; BESTMIN - upon return, the best minimum function value that TNMIN ; could find. ; ; EPSABS - a nonnegative real variable. Termination occurs when the ; absolute error between consecutive iterates is at most ; EPSABS. Note that using EPSREL is normally preferable ; over EPSABS, except in cases where ABS(F) is much larger ; than 1 at the optimal point. A value of zero indicates ; the absolute error test is not applied. If EPSABS is ; specified, then both EPSREL and EPSABS tests are applied; ; if either succeeds then termination occurs. ; Default: 0 (i.e., only EPSREL is applied). ; ; EPSREL - a nonnegative input variable. Termination occurs when the ; relative error between two consecutive iterates is at ; most EPSREL. Therefore, EPSREL measures the relative ; error desired in the function. An additional, more ; lenient, stopping condition can be applied by specifying ; the EPSABS keyword. ; Default: 100 * Machine Precision ; ; ERRMSG - a string error or warning message is returned. ; ; FGUESS - provides the routine a guess to the minimum value. ; Default: 0 ; ; 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]} ; then the user supplied function should be declared ; like this: ; FUNCTION MYFUNCT, P, XVAL=x, YVAL=y ; ; By default, no extra parameters are passed to the ; user-supplied function. ; ; 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. ; ; ITERDERIV - Intended to print function gradient information. If ; set, then the ITERDERIV keyword is set in each call to ; ITERPROC. In the default ITERPROC, parameter values ; and gradient information are both printed when this ; keyword is set. ; ; ITERPROC - The name of a procedure to be called upon each NPRINT ; iteration of the TNMIN routine. It should be declared ; in the following way: ; ; PRO ITERPROC, MYFUNCT, p, iter, fnorm, FUNCTARGS=fcnargs, $ ; PARINFO=parinfo, QUIET=quiet, _EXTRA=extra ; ; perform custom iteration update ; END ; ; ITERPROC must accept the _EXTRA keyword, in case of ; future changes to the calling procedure. ; ; 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 is should be the function ; value. QUIET is set when no textual output should be ; printed. See below for documentation of PARINFO. ; ; In implementation, ITERPROC can perform updates to the ; terminal or graphical user interface, to provide ; feedback while the fit proceeds. If the fit is to be ; stopped for any reason, then ITERPROC should set the ; common block variable ERROR_CODE to negative value ; between -15 and -1 (see TNMIN_ERROR common block ; below). In principle, ITERPROC should probably not ; modify the parameter values, because it may interfere ; with the algorithm's stability. In practice it is ; allowed. ; ; Default: an internal routine is used to print the ; parameter values. ; ; MAXITER - The maximum number of iterations to perform. If the ; number is exceeded, then the STATUS value is set to 5 ; and TNMIN returns. ; Default: 200 iterations ; ; MAXIMIZE - If set, the function is maximized instead of ; minimized. ; ; MAXNFEV - The maximum number of function evaluations to perform. ; If the number is exceeded, then the STATUS value is set ; to -17 and TNMIN returns. A value of zero indicates no ; maximum. ; Default: 0 (no maximum) ; ; NFEV - upon return, the number of MYFUNCT function evaluations ; performed. ; ; NITER - upon return, number of iterations completed. ; ; 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. ; Default value: 1 ; ; PARINFO - 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 TNMIN. ; ; See description above for the structure of PARINFO. ; ; Default value: all parameters are free and unconstrained. ; ; QUIET - set this keyword when no textual output should be printed ; by TNMIN ; ; STATUS - an integer status code is returned. All values greater ; than zero can represent success (however STATUS EQ 5 may ; indicate failure to converge). Gaps in the numbering ; system below are to maintain compatibility with MPFIT. ; Upon return, STATUS can have one of the following values: ; ; -18 a fatal internal error occurred during optimization. ; ; -17 the maximum number of function evaluations has been ; reached without convergence. ; ; -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 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 TNMIN_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 TNMIN. ; ; 0 improper input parameters. ; ; 1 convergence was reached. ; ; 2-4 (RESERVED) ; ; 5 the maximum number of iterations has been reached ; ; 6-8 (RESERVED) ; ; ; EXAMPLE: ; ; FUNCTION F, X, DF, _EXTRA=extra ;; *** MUST ACCEPT KEYWORDS ; F = (X(0)-1)^2 + (X(1)+7)^2 ; DF = [ 2D * (X(0)-1), 2D * (X(1)+7) ] ; Gradient ; RETURN, F ; END ; ; P = TNMIN('F', [0D, 0D], BESTMIN=F0) ; Minimizes the function F(x0,x1) = (x0-1)^2 + (x1+7)^2. ; ; ; COMMON BLOCKS: ; ; COMMON TNMIN_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, TNMIN checks the value of this ; common block variable and exits immediately if the error ; condition has been set. By default the value of ERROR_CODE is ; zero, indicating a successful function/procedure call. ; ; ; REFERENCES: ; ; TRUNCATED-NEWTON METHOD, TN.F ; Stephen G. Nash, Operations Research and Applied Statistics ; Department ; http://www.netlib.org/opt/tn ; ; Nash, S. G. 1984, "Newton-Type Minimization via the Lanczos ; Method," SIAM J. Numerical Analysis, 21, p. 770-778 ; ; ; MODIFICATION HISTORY: ; Derived from TN.F by Stephen Nash with many changes and additions, ; to conform to MPFIT paradigm, 19 Jan 1999, CM ; Changed web address to COW, 25 Sep 1999, CM ; Alphabetized documented keyword parameters, 02 Oct 1999, CM ; Changed to ERROR_CODE for error condition, 28 Jan 2000, CM ; Continued and fairly major improvements (CM, 08 Jan 2001): ; - calling of user procedure is now concentrated in TNMIN_CALL, ; and arguments are reduced by storing a large number of them ; in common blocks; ; - finite differencing is done within TNMIN_CALL; added ; AUTODERIVATIVE=1 so that finite differencing can be enabled, ; both one and two sided; ; - a new procedure to parse PARINFO fields, borrowed from MPFIT; ; brought PARINFO keywords up to date with MPFIT; ; - go through and check for float vs. double discrepancies; ; - add explicit MAXIMIZE keyword, and support in TNMIN_CALL and ; TNMIN_DEFITER to print the correct values in that case; ; TNMIN_DEFITER now prints function gradient values if ; requested; ; - convert to common-based system of MPFIT for storing machine ; constants; revert TNMIN_ENORM to simple sum of squares, at ; least for now; ; - remove limit on number of function evaluations, at least for ; now, and until I can understand what happens when we do ; numerical derivatives. ; Further changes: more float vs double; disable TNMINSTEP for now; ; experimented with convergence test in case of function ; maximization, 11 Jan 2001, CM ; TNMINSTEP is parsed but not enabled, 13 Mar 2001 ; Major code cleanups; internal docs; reduced commons, CM, 08 Apr ; 2001 ; Continued code cleanups; documentation; the STATUS keyword ; actually means something, CM, 10 Apr 2001 ; Added reference to Nash paper, CM, 08 Feb 2002 ; Fixed 16-bit loop indices, D. Schelgel, 22 Apr 2003 ; Changed parens to square brackets because of conflicts with ; limits function. K. Tolbert, 23 Feb 2005 ; Some documentation clarifications, CM, 09 Nov 2007 ; Ensure that MY_FUNCT returns a scalar; make it more likely that ; error messages get back out to the user; fatal CATCH'd error ; now returns STATUS = -18, CM, 17 Sep 2008 ; Fix TNMIN_CALL when parameters are TIEd (thanks to Alfred de ; Wijn), CM, 22 Nov 2009 ; Remember to TIE the parameters before final return (thanks to ; Michael Smith), CM, 20 Jan 2010 ; ; TODO ; - scale derivatives semi-automatically; ; - ability to scale and offset parameters; ; ; $Id: tnmin.pro,v 1.20 2016/05/19 16:08:08 cmarkwar Exp $ ;- ; Copyright (C) 1998-2001,2002,2003,2007,2008,2009 Craig Markwardt ; This software is provided as is without any warranty whatsoever. ; Permission to use, copy and distribute unmodified copies for ; non-commercial purposes, and to modify and use for personal or ; internal use, is granted. All other rights are reserved. ;- ;%% TRUNCATED-NEWTON METHOD: SUBROUTINES ; FOR OTHER MACHINES, MODIFY ROUTINE MCHPR1 (MACHINE EPSILON) ; WRITTEN BY: STEPHEN G. NASH ; OPERATIONS RESEARCH AND APPLIED STATISTICS DEPT. ; GEORGE MASON UNIVERSITY ; FAIRFAX, VA 22030 ;****************************************************************** ;; Routine which declares functions and common blocks pro tnmin_dummy forward_function tnmin_enorm, tnmin_step1, tnmin forward_function tnmin_call, tnmin_autoder common tnmin_error, error_code common tnmin_machar, tnmin_machar_vals common tnmin_config, tnmin_tnconfig common tnmin_fcnargs, tnmin_tnfcnargs common tnmin_work, lsk, lyk, ldiagb, lsr, lyr a = 1 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 tnmin_setmachar, double=isdouble common tnmin_machar, tnmin_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 tnmin_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 ;; Procedure to parse the parameter values in PARINFO pro tnmin_parinfo, parinfo, tnames, tag, values, default=def, status=status, $ n_param=n 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)) endelse status = 1 return end ;; Procedure to tie one parameter to another. pro tnmin_tie, p, _ptied _wh = where(_ptied NE '', _ct) if _ct EQ 0 then return for _i = 0L, _ct-1 do begin _cmd = 'p('+strtrim(_wh(_i),2)+') = '+_ptied(_wh(_i)) _err = execute(_cmd) if _err EQ 0 then begin message, 'ERROR: Tied expression "'+_cmd+'" failed.' return endif endfor end function tnmin_autoder, fcn, x, dx, dside=dside common tnmin_machar, machvals common tnmin_config, tnconfig MACHEP0 = machvals.machep DWARF = machvals.minnum if n_elements(dside) NE n_elements(x) then dside = tnconfig.dside eps = sqrt(MACHEP0) h = eps * (1. + abs(x)) ;; if STEP is given, use that wh = where(tnconfig.step GT 0, ct) if ct GT 0 then h(wh) = tnconfig.step(wh) ;; if relative step is given, use that wh = where(tnconfig.dstep GT 0, ct) if ct GT 0 then h(wh) = abs(tnconfig.dstep(wh)*x(wh)) ;; 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. mask = (dside EQ -1 OR (tnconfig.ulimited AND (x GT tnconfig.ulimit-h))) wh = where(mask, ct) if ct GT 0 then h(wh) = -h(wh) dx = x * 0. f = tnmin_call(fcn, x) for i = 0L, n_elements(x)-1 do begin if tnconfig.pfixed(i) EQ 1 then goto, NEXT_PAR hh = h(i) RESTART_PAR: xp = x xp(i) = xp(i) + hh fp = tnmin_call(fcn, xp) if abs(dside(i)) LE 1 then begin ;; COMPUTE THE ONE-SIDED DERIVATIVE dx(i) = (fp-f)/hh endif else begin ;; COMPUTE THE TWO-SIDED DERIVATIVE xp(i) = x(i) - hh fm = tnmin_call(fcn, xp) dx(i) = (fp-fm)/(2*hh) endelse NEXT_PAR: endfor return, f end ;; Call user function or procedure, with _EXTRA or not, with ;; derivatives or not. function tnmin_call, fcn, x1, dx, fullparam_=xall ; on_error, 2 common tnmin_config, tnconfig common tnmin_fcnargs, fcnargs ifree = tnconfig.ifree ;; Following promotes the byte array to a floating point array so ;; that users who simply re-fill the array aren't surprised when ;; their gradient comes out as bytes. :-) dx = tnconfig.pfixed + x1(0)*0. if n_elements(xall) GT 0 then begin x = xall x(ifree) = x1 endif else begin x = x1 endelse ;; Enforce TIEd parameters if keyword_set(tnconfig.qanytied) then tnmin_tie, x, tnconfig.ptied ;; Decide whether we are calling a procedure or function if tnconfig.proc then proc = 1 else proc = 0 tnconfig.nfev = tnconfig.nfev + 1 if n_params() EQ 3 then begin if tnconfig.autoderiv then $ f = tnmin_autoder(fcn, x, dx) $ else if n_elements(fcnargs) GT 0 then $ f = call_function(fcn, x, dx, _EXTRA=fcnargs) $ else $ f = call_function(fcn, x, dx) dx = dx(ifree) if tnconfig.max then begin dx = -dx f = -f endif endif else begin if n_elements(fcnargs) GT 0 then $ f = call_function(fcn, x, _EXTRA=fcnargs) $ else $ f = call_function(fcn, x) if n_elements(f) NE 1 then begin message, 'ERROR: function "'+fcn+'" returned a vector when '+$ 'a scalar was expected.' endif endelse if n_elements(f) GT 1 then return, f $ else return, f(0) end function tnmin_enorm, vec common tnmin_config, tnconfig ;; Very simple-minded sum-of-squares if n_elements(tnconfig) GT 0 then if tnconfig.fastnorm then begin ans = sqrt(total(vec^2,1)) goto, TERMINATE endif common tnmin_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: return, ans end ; ; ROUTINES TO INITIALIZE PRECONDITIONER ; pro tnmin_initpc, diagb, emat, n, upd1, yksk, gsk, yrsr, lreset ;; Rename common variables as they appear in INITP3. Those ;; indicated in all caps are not used or renamed here. ; common tnmin_work, lsk, lyk, ldiagb, lsr, lyr common tnmin_work, sk, yk, LDIAGB, sr, yr ; I I I I ;; From INITP3 if keyword_set(upd1) then begin EMAT = DIAGB endif else if keyword_set(lreset) then begin BSK = DIAGB*SK SDS = TOTAL(SK*BSK) EMAT = DIAGB - DIAGB*DIAGB*SK*SK/SDS + YK*YK/YKSK endif else begin BSK = DIAGB * SR SDS = TOTAL(SR*BSK) SRDS = TOTAL(SK*BSK) YRSK = TOTAL(YR*SK) BSK = DIAGB*SK - BSK*SRDS/SDS+YR*YRSK/YRSR EMAT = DIAGB-DIAGB*DIAGB*SR*SR/SDS+YR*YR/YRSR SDS = TOTAL(SK*BSK) EMAT = EMAT - BSK*BSK/SDS+YK*YK/YKSK endelse return end pro tnmin_ssbfgs, n, gamma, sj, yj, hjv, hjyj, yjsj, yjhyj, $ vsj, vhyj, hjp1v ; ; SELF-SCALED BFGS ; DELTA = (1. + GAMMA*YJHYJ/YJSJ)*VSJ/YJSJ - GAMMA*VHYJ/YJSJ BETA = -GAMMA*VSJ/YJSJ HJP1V = GAMMA*HJV + DELTA*SJ + BETA*HJYJ RETURN end ; ; THIS ROUTINE ACTS AS A PRECONDITIONING STEP FOR THE ; LINEAR CONJUGATE-GRADIENT ROUTINE. IT IS ALSO THE ; METHOD OF COMPUTING THE SEARCH DIRECTION FROM THE ; GRADIENT FOR THE NON-LINEAR CONJUGATE-GRADIENT CODE. ; IT REPRESENTS A TWO-STEP SELF-SCALED BFGS FORMULA. ; pro tnmin_msolve, g, y, n, upd1, yksk, gsk, yrsr, lreset, first, $ hyr, hyk, ykhyk, yrhyr ;; Rename common variables as they appear in MSLV ; common tnmin_work, lsk, lyk, ldiagb, lsr, lyr common tnmin_work, sk, yk, diagb, sr, yr ; I I I I I ;; From MSLV if keyword_set(UPD1) then begin Y = G / DIAGB RETURN endif ONE = G(0)*0 + 1. GSK = TOTAL(G*SK) if keyword_set(lreset) then begin ; ; COMPUTE GH AND HY WHERE H IS THE INVERSE OF THE DIAGONALS ; HG = G/DIAGB if keyword_set(FIRST) then begin HYK = YK/DIAGB YKHYK = TOTAL(YK*HYK) endif GHYK = TOTAL(G*HYK) TNMIN_SSBFGS,N,ONE,SK,YK,HG,HYK,YKSK, YKHYK,GSK,GHYK,Y RETURN endif ; ; COMPUTE HG AND HY WHERE H IS THE INVERSE OF THE DIAGONALS ; HG = G/DIAGB if keyword_set(FIRST) then begin HYK = YK/DIAGB HYR = YR/DIAGB YKSR = TOTAL(YK*SR) YKHYR = TOTAL(YK*HYR) endif GSR = TOTAL(G*SR) GHYR = TOTAL(G*HYR) if keyword_set(FIRST) then begin YRHYR = TOTAL(YR*HYR) endif TNMIN_SSBFGS,N,ONE,SR,YR,HG,HYR,YRSR, YRHYR,GSR,GHYR,HG if keyword_set(FIRST) then begin TNMIN_SSBFGS,N,ONE,SR,YR,HYK,HYR,YRSR, YRHYR,YKSR,YKHYR,HYK endif YKHYK = TOTAL(HYK*YK) GHYK = TOTAL(HYK*G) TNMIN_SSBFGS,N,ONE,SK,YK,HG,HYK,YKSK, YKHYK,GSK,GHYK,Y RETURN end ; ; THIS ROUTINE COMPUTES THE PRODUCT OF THE MATRIX G TIMES THE VECTOR ; V AND STORES THE RESULT IN THE VECTOR GV (FINITE-DIFFERENCE VERSION) ; pro tnmin_gtims, v, gv, n, x, g, fcn, first, delta, accrcy, xnorm, $ xnew IF keyword_set(FIRST) THEN BEGIN ;; Extra factor of ten is to avoid clashing with the finite ;; difference scheme which computes the derivatives DELTA = SQRT(100*ACCRCY)*(1.+XNORM) ;; XXX diff than TN.F ; DELTA = SQRT(ACCRCY)*(1.+XNORM) FIRST = 0 ENDIF DINV = 1. /DELTA F = tnmin_call(FCN, X + DELTA*V, GV, fullparam_=xnew) GV = (GV-G)*DINV return end ; ; UPDATE THE PRECONDITIOING MATRIX BASED ON A DIAGONAL VERSION ; OF THE BFGS QUASI-NEWTON UPDATE. ; pro tnmin_ndia3, n, e, v, gv, r, vgv VR = TOTAL(V*R) E = E - R*R/VR + GV*GV/VGV wh = where(e LE 1D-6, ct) if ct GT 0 then e(wh) = 1 return end pro tnmin_fix, whlpeg, whupeg, z if whlpeg(0) NE -1 then z(whlpeg) = 0 if whupeg(0) NE -1 then z(whupeg) = 0 end ; ; THIS ROUTINE PERFORMS A PRECONDITIONED CONJUGATE-GRADIENT ; ITERATION IN ORDER TO SOLVE THE NEWTON EQUATIONS FOR A SEARCH ; DIRECTION FOR A TRUNCATED-NEWTON ALGORITHM. WHEN THE VALUE OF THE ; QUADRATIC MODEL IS SUFFICIENTLY REDUCED, ; THE ITERATION IS TERMINATED. ; ; PARAMETERS ; ; ZSOL - COMPUTED SEARCH DIRECTION ; G - CURRENT GRADIENT ; GV,GZ1,V - SCRATCH VECTORS ; R - RESIDUAL ; DIAGB,EMAT - DIAGONAL PRECONDITONING MATRIX ; NITER - NONLINEAR ITERATION # ; FEVAL - VALUE OF QUADRATIC FUNCTION pro tnmin_modlnp, zsol, gv, r, v, diagb, emat, $ x, g, zk, n, niter, maxit, nmodif, nlincg, $ upd1, yksk, gsk, yrsr, lreset, fcn, whlpeg, whupeg, $ accrcy, gtp, gnorm, xnorm, xnew ; ; GENERAL INITIALIZATION ; zero = x(0)* 0. one = zero + 1 IF (MAXIT EQ 0) THEN RETURN FIRST = 1 RHSNRM = GNORM TOL = zero + 1.E-12 QOLD = zero ; ; INITIALIZATION FOR PRECONDITIONED CONJUGATE-GRADIENT ALGORITHM ; tnmin_initpc, diagb, emat, n, upd1, yksk, gsk, yrsr, lreset R = -G V = G*0. ZSOL = V ; ; ************************************************************ ; MAIN ITERATION ; ************************************************************ ; FOR K = 1L, MAXIT DO BEGIN NLINCG = NLINCG + 1 ; ; CG ITERATION TO SOLVE SYSTEM OF EQUATIONS ; tnmin_fix, whlpeg, whupeg, r TNMIN_MSOLVE, R, ZK, N, UPD1, YKSK, GSK, YRSR, LRESET, FIRST, $ HYR, HYK, YKHYK, YRHYR tnmin_fix, whlpeg, whupeg, zk RZ = TOTAL(R*ZK) IF (RZ/RHSNRM LT TOL) THEN GOTO, MODLNP_80 IF (K EQ 1) THEN BETA = ZERO $ ELSE BETA = RZ/RZOLD V = ZK + BETA*V tnmin_fix, whlpeg, whupeg, v TNMIN_GTIMS, V, GV, N, X, G, FCN, FIRST, DELTA, ACCRCY, XNORM, XNEW tnmin_fix, whlpeg, whupeg, gv VGV = TOTAL(V*GV) IF (VGV/RHSNRM LT TOL) THEN GOTO, MODLNP_50 TNMIN_NDIA3, N,EMAT,V,GV,R,VGV ; ; COMPUTE LINEAR STEP LENGTH ; ALPHA = RZ / VGV ; ; COMPUTE CURRENT SOLUTION AND RELATED VECTORS ; ZSOL = ZSOL + ALPHA*V R = R - ALPHA*GV ; ; TEST FOR CONVERGENCE ; GTP = TOTAL(ZSOL*G) PR = TOTAL(R*ZSOL) QNEW = 5.E-1 * (GTP + PR) QTEST = K * (1.E0 - QOLD/QNEW) IF (QTEST LT 0.D0) THEN GOTO, MODLNP_70 QOLD = QNEW IF (QTEST LE 5.D-1) THEN GOTO, MODLNP_70 ; ; PERFORM CAUTIONARY TEST ; IF (GTP GT 0) THEN GOTO, MODLNP_40 RZOLD = RZ ENDFOR ; ; TERMINATE ALGORITHM ; K = K-1 goto, MODLNP_70 MODLNP_40: ZSOL = ZSOL - ALPHA*V GTP = TOTAL(ZSOL*G) goto, MODLNP_90 MODLNP_50: ;; printed output MODLNP_60: IF (K GT 1) THEN GOTO, MODLNP_70 TNMIN_MSOLVE,G,ZSOL,N,UPD1,YKSK,GSK,YRSR,LRESET,FIRST, $ HYR, HYK, YKHYK, YRHYR ZSOL = -ZSOL tnmin_fix, whlpeg, whupeg, zsol GTP = TOTAL(ZSOL*G) MODLNP_70: goto, MODLNP_90 MODLNP_80: IF (K GT 1) THEN GOTO, MODLNP_70 ZSOL = -G tnmin_fix, whlpeg, whupeg, zsol GTP = TOTAL(ZSOL*G) goto, MODLNP_70 ; ; STORE (OR RESTORE) DIAGONAL PRECONDITIONING ; MODLNP_90: diagb = emat return end function tnmin_step1, fnew, fm, gtp, smax, epsmch ; ******************************************************** ; STEP1 RETURNS THE LENGTH OF THE INITIAL STEP TO BE TAKEN ALONG THE ; VECTOR P IN THE NEXT LINEAR SEARCH. ; ******************************************************** D = ABS(FNEW-FM) ALPHA = FNEW(0)*0 + 1. IF (2.D0*D LE (-GTP) AND D GE EPSMCH) THEN $ ALPHA = -2.*D/GTP IF (ALPHA GE SMAX) THEN ALPHA = SMAX return, alpha end ; ; ************************************************************ ; GETPTC, AN ALGORITHM FOR FINDING A STEPLENGTH, CALLED REPEATEDLY BY ; ROUTINES WHICH REQUIRE A STEP LENGTH TO BE COMPUTED USING CUBIC ; INTERPOLATION. THE PARAMETERS CONTAIN INFORMATION ABOUT THE INTERVAL ; IN WHICH A LOWER POINT IS TO BE FOUND AND FROM THIS GETPTC COMPUTES A ; POINT AT WHICH THE FUNCTION CAN BE EVALUATED BY THE CALLING PROGRAM. ; THE VALUE OF THE INTEGER PARAMETERS IENTRY DETERMINES THE PATH TAKEN ; THROUGH THE CODE. ; ************************************************************ pro tnmin_getptc, big, small, rtsmll, reltol, abstol, tnytol, $ fpresn, eta, rmu, xbnd, u, fu, gu, xmin, fmin, gmin, $ xw, fw, gw, a, b, oldf, b1, scxbnd, e, step, factor, $ braktd, gtest1, gtest2, tol, ientry, itest ;; This chicanery is so that we get the data types right ZERO = fu(0)* 0. ; a1 = zero & scale = zero & chordm = zero ; chordu = zero & d1 = zero & d2 = zero ; denom = zero POINT1 = ZERO + 0.1 HALF = ZERO + 0.5 ONE = ZERO + 1 THREE = ZERO + 3 FIVE = ZERO + 5 ELEVEN = ZERO + 11 if ientry EQ 1 then begin ;; else clause = 20 (OK) ; ; IENTRY=1 ; CHECK INPUT PARAMETERS ; ;; GETPTC_10: ITEST = 2 IF (U LE ZERO OR XBND LE TNYTOL OR GU GT ZERO) THEN RETURN ITEST = 1 IF (XBND LT ABSTOL) THEN ABSTOL = XBND TOL = ABSTOL TWOTOL = TOL + TOL ; ; A AND B DEFINE THE INTERVAL OF UNCERTAINTY, X AND XW ARE POINTS ; WITH LOWEST AND SECOND LOWEST FUNCTION VALUES SO FAR OBTAINED. ; INITIALIZE A,SMIN,XW AT ORIGIN AND CORRESPONDING VALUES OF ; FUNCTION AND PROJECTION OF THE GRADIENT ALONG DIRECTION OF SEARCH ; AT VALUES FOR LATEST ESTIMATE AT MINIMUM. ; A = ZERO XW = ZERO XMIN = ZERO OLDF = FU FMIN = FU FW = FU GW = GU GMIN = GU STEP = U FACTOR = FIVE ; ; THE MINIMUM HAS NOT YET BEEN BRACKETED. ; BRAKTD = 0 ; ; SET UP XBND AS A BOUND ON THE STEP TO BE TAKEN. (XBND IS NOT COMPUTED ; EXPLICITLY BUT SCXBND IS ITS SCALED VALUE.) SET THE UPPER BOUND ; ON THE INTERVAL OF UNCERTAINTY INITIALLY TO XBND + TOL(XBND). ; SCXBND = XBND B = SCXBND + RELTOL*ABS(SCXBND) + ABSTOL E = B + B B1 = B ; ; COMPUTE THE CONSTANTS REQUIRED FOR THE TWO CONVERGENCE CRITERIA. ; GTEST1 = -RMU*GU GTEST2 = -ETA*GU ; ; SET IENTRY TO INDICATE THAT THIS IS THE FIRST ITERATION ; IENTRY = 2 goto, GETPTC_210 endif ; ; IENTRY = 2 ; ; UPDATE A,B,XW, AND XMIN ; ;; GETPTC_20: IF (FU GT FMIN) THEN GOTO, GETPTC_60 ; ; IF FUNCTION VALUE NOT INCREASED, NEW POINT BECOMES NEXT ; ORIGIN AND OTHER POINTS ARE SCALED ACCORDINGLY. ; CHORDU = OLDF - (XMIN + U)*GTEST1 if NOT (FU LE CHORDU) then begin ; ; THE NEW FUNCTION VALUE DOES NOT SATISFY THE SUFFICIENT DECREASE ; CRITERION. PREPARE TO MOVE THE UPPER BOUND TO THIS POINT AND ; FORCE THE INTERPOLATION SCHEME TO EITHER BISECT THE INTERVAL OF ; UNCERTAINTY OR TAKE THE LINEAR INTERPOLATION STEP WHICH ESTIMATES ; THE ROOT OF F(ALPHA)=CHORD(ALPHA). ; CHORDM = OLDF - XMIN*GTEST1 GU = -GMIN DENOM = CHORDM-FMIN IF (ABS(DENOM) LT 1.D-15) THEN BEGIN DENOM = ZERO + 1.E-15 IF (CHORDM-FMIN LT 0.D0) THEN DENOM = -DENOM ENDIF IF (XMIN NE ZERO) THEN GU = GMIN*(CHORDU-FU)/DENOM FU = (HALF*U*(GMIN+GU) + FMIN) > FMIN ; ; IF FUNCTION VALUE INCREASED, ORIGIN REMAINS UNCHANGED ; BUT NEW POINT MAY NOW QUALIFY AS W. ; GETPTC_60: IF (U GE ZERO) THEN BEGIN B = U BRAKTD = 1 ENDIF ELSE BEGIN A = U ENDELSE XW = U FW = FU GW = GU endif else begin ;; GETPTC_30: FW = FMIN FMIN = FU GW = GMIN GMIN = GU XMIN = XMIN + U A = A-U B = B-U XW = -U SCXBND = SCXBND - U IF (GU GT ZERO) THEN BEGIN B = ZERO BRAKTD = 1 ENDIF ELSE BEGIN A = ZERO ENDELSE TOL = ABS(XMIN)*RELTOL + ABSTOL endelse TWOTOL = TOL + TOL XMIDPT = HALF*(A + B) ; ; CHECK TERMINATION CRITERIA ; CONVRG = ABS(XMIDPT) LE TWOTOL - HALF*(B-A) OR $ ABS(GMIN) LE GTEST2 AND FMIN LT OLDF AND $ (ABS(XMIN - XBND) GT TOL OR NOT BRAKTD) IF CONVRG THEN BEGIN ITEST = 0 IF (XMIN NE ZERO) THEN RETURN ; ; IF THE FUNCTION HAS NOT BEEN REDUCED, CHECK TO SEE THAT THE RELATIVE ; CHANGE IN F(X) IS CONSISTENT WITH THE ESTIMATE OF THE DELTA- ; UNIMODALITY CONSTANT, TOL. IF THE CHANGE IN F(X) IS LARGER THAN ; EXPECTED, REDUCE THE VALUE OF TOL. ; ITEST = 3 IF (ABS(OLDF-FW) LE FPRESN*(ONE + ABS(OLDF))) THEN RETURN TOL = POINT1*TOL IF (TOL LT TNYTOL) THEN RETURN RELTOL = POINT1*RELTOL ABSTOL = POINT1*ABSTOL TWOTOL = POINT1*TWOTOL endif ; ; CONTINUE WITH THE COMPUTATION OF A TRIAL STEP LENGTH ; ;; GETPTC_100: R = ZERO Q = ZERO S = ZERO IF (ABS(E) GT TOL) THEN BEGIN ; ; FIT CUBIC THROUGH XMIN AND XW ; R = THREE*(FMIN-FW)/XW + GMIN + GW ABSR = ABS(R) Q = ABSR IF (GW EQ ZERO OR GMIN EQ ZERO) EQ 0 THEN BEGIN ;; else clause = 140 (OK) ; ; COMPUTE THE SQUARE ROOT OF (R*R - GMIN*GW) IN A WAY ; WHICH AVOIDS UNDERFLOW AND OVERFLOW. ; ABGW = ABS(GW) ABGMIN = ABS(GMIN) S = SQRT(ABGMIN)*SQRT(ABGW) IF ((GW/ABGW)*GMIN LE ZERO) THEN BEGIN ; ; COMPUTE THE SQUARE ROOT OF R*R + S*S. ; SUMSQ = ONE P = ZERO IF (ABSR LT S) THEN BEGIN ;; else clause = 110 (OK) ; ; THERE IS A POSSIBILITY OF OVERFLOW. ; IF (S GT RTSMLL) THEN P = S*RTSMLL IF (ABSR GE P) THEN SUMSQ = ONE +(ABSR/S)^2 SCALE = S endif else begin ;; flow to 120 (OK) ; ; THERE IS A POSSIBILITY OF UNDERFLOW. ; ;; GETPTC_110: IF (ABSR GT RTSMLL) THEN P = ABSR*RTSMLL IF (S GE P) THEN SUMSQ = ONE + (S/ABSR)^2 SCALE = ABSR ENDELSE ;; flow to 120 (OK) ;; GETPTC_120: SUMSQ = SQRT(SUMSQ) Q = BIG IF (SCALE LT BIG/SUMSQ) THEN Q = SCALE*SUMSQ endif else begin ;; flow to 140 ; ; COMPUTE THE SQUARE ROOT OF R*R - S*S ; ;; GETPTC_130: Q = SQRT(ABS(R+S))*SQRT(ABS(R-S)) IF (R GE S OR R LE (-S)) EQ 0 THEN BEGIN R = ZERO Q = ZERO goto, GETPTC_150 endif endelse endif ; ; COMPUTE THE MINIMUM OF FITTED CUBIC ; ;; GETPTC_140: IF (XW LT ZERO) THEN Q = -Q S = XW*(GMIN - R - Q) Q = GW - GMIN + Q + Q IF (Q GT ZERO) THEN S = -S IF (Q LE ZERO) THEN Q = -Q R = E IF (B1 NE STEP OR BRAKTD) THEN E = STEP endif ; ; CONSTRUCT AN ARTIFICIAL BOUND ON THE ESTIMATED STEPLENGTH ; GETPTC_150: A1 = A B1 = B STEP = XMIDPT IF (BRAKTD) EQ 0 THEN BEGIN ;; else flow to 160 (OK) STEP = -FACTOR*XW IF (STEP GT SCXBND) THEN STEP = SCXBND IF (STEP NE SCXBND) THEN FACTOR = FIVE*FACTOR ;; flow to 170 (OK) endif else begin ; ; IF THE MINIMUM IS BRACKETED BY 0 AND XW THE STEP MUST LIE ; WITHIN (A,B). ; ;; GETPTC_160: if (a NE zero OR xw GE zero) AND (b NE zero OR xw LE zero) then $ goto, GETPTC_180 ; ; IF THE MINIMUM IS NOT BRACKETED BY 0 AND XW THE STEP MUST LIE ; WITHIN (A1,B1). ; D1 = XW D2 = A IF (A EQ ZERO) THEN D2 = B ; THIS LINE MIGHT BE ; IF (A EQ ZERO) THEN D2 = E U = - D1/D2 STEP = FIVE*D2*(POINT1 + ONE/U)/ELEVEN IF (U LT ONE) THEN STEP = HALF*D2*SQRT(U) endelse ;; GETPTC_170: IF (STEP LE ZERO) THEN A1 = STEP IF (STEP GT ZERO) THEN B1 = STEP ; ; REJECT THE STEP OBTAINED BY INTERPOLATION IF IT LIES OUTSIDE THE ; REQUIRED INTERVAL OR IT IS GREATER THAN HALF THE STEP OBTAINED ; DURING THE LAST-BUT-ONE ITERATION. ; GETPTC_180: if NOT (abs(s) LE abs(half*q*r) OR s LE q*a1 OR s GE q*b1) then begin ;; else clause = 200 (OK) ; ; A CUBIC INTERPOLATION STEP ; STEP = S/Q ; ; THE FUNCTION MUST NOT BE EVALUTATED TOO CLOSE TO A OR B. ; if NOT (step - a GE twotol AND b - step GE twotol) then begin ;; else clause = 210 (OK) IF (XMIDPT LE ZERO) THEN STEP = -TOL ELSE STEP = TOL endif ;; flow to 210 (OK) endif else begin ;; GETPTC_200: E = B-A endelse ; ; IF THE STEP IS TOO LARGE, REPLACE BY THE SCALED BOUND (SO AS TO ; COMPUTE THE NEW POINT ON THE BOUNDARY). ; GETPTC_210: if (step GE scxbnd) then begin ;; else clause = 220 (OK) STEP = SCXBND ; ; MOVE SXBD TO THE LEFT SO THAT SBND + TOL(XBND) = XBND. ; SCXBND = SCXBND - (RELTOL*ABS(XBND)+ABSTOL)/(ONE + RELTOL) endif ;; GETPTC_220: U = STEP IF (ABS(STEP) LT TOL AND STEP LT ZERO) THEN U = -TOL IF (ABS(STEP) LT TOL AND STEP GE ZERO) THEN U = TOL ITEST = 1 RETURN end ; ; LINE SEARCH ALGORITHMS OF GILL AND MURRAY ; pro tnmin_linder, n, fcn, small, epsmch, reltol, abstol, $ tnytol, eta, sftbnd, xbnd, p, gtp, x, f, alpha, g, $ iflag, xnew zero = f(0) * 0. one = zero + 1. LSPRNT = 0L NPRNT = 10000L RTSMLL = SQRT(SMALL) BIG = 1./SMALL ITCNT = 0L ; ; SET THE ESTIMATED RELATIVE PRECISION IN F(X). ; FPRESN = 10.*EPSMCH U = ALPHA FU = F FMIN = F GU = GTP RMU = zero + 1E-4 ; ; FIRST ENTRY SETS UP THE INITIAL INTERVAL OF UNCERTAINTY. ; IENTRY = 1L LINDER_10: ; ; TEST FOR TOO MANY ITERATIONS ; ITCNT = ITCNT + 1 IF (ITCNT GT 30) THEN BEGIN ;; deviation from Nash: allow optimization to continue in outer ;; loop even if we fail to converge, if IFLAG EQ 0. A value of ;; 1 indicates failure. I believe that I tried IFLAG=0 once and ;; there was some problem, but I forget what it was. IFLAG = 1 F = FMIN ALPHA = XMIN X = X + ALPHA*P RETURN ENDIF IFLAG = 0 TNMIN_GETPTC,BIG,SMALL,RTSMLL,RELTOL,ABSTOL,TNYTOL, $ FPRESN,ETA,RMU,XBND,U,FU,GU,XMIN,FMIN,GMIN, $ XW,FW,GW,A,B,OLDF,B1,SCXBND,E,STEP,FACTOR, $ BRAKTD,GTEST1,GTEST2,TOL,IENTRY,ITEST ; ; IF ITEST=1, THE ALGORITHM REQUIRES THE FUNCTION VALUE TO BE ; CALCULATED. ; IF (ITEST EQ 1) THEN BEGIN UALPHA = XMIN + U FU = TNMIN_CALL(FCN, X + UALPHA*P, LG, fullparam_=xnew) GU = TOTAL(LG*P) ; ; THE GRADIENT VECTOR CORRESPONDING TO THE BEST POINT IS ; OVERWRITTEN IF FU IS LESS THAN FMIN AND FU IS SUFFICIENTLY ; LOWER THAN F AT THE ORIGIN. ; IF (FU LE FMIN AND FU LE OLDF-UALPHA*GTEST1) THEN $ G = LG ; print, 'fu = ', fu GOTO, LINDER_10 ENDIF ; ; IF ITEST=2 OR 3 A LOWER POINT COULD NOT BE FOUND ; IFLAG = 1 IF (ITEST NE 0) THEN RETURN ; ; IF ITEST=0 A SUCCESSFUL SEARCH HAS BEEN MADE ; ; print, 'itcnt = ', itcnt IFLAG = 0 F = FMIN ALPHA = XMIN X = X + ALPHA*P RETURN END pro tnmin_defiter, fcn, x, iter, fnorm, fmt=fmt, FUNCTARGS=fcnargs, $ quiet=quiet, deriv=df, dprint=dprint, pfixed=pfixed, $ maximize=maximize, _EXTRA=iterargs if keyword_set(quiet) then return if n_params() EQ 3 then begin fnorm = tnmin_call(fcn, x, df) endif if keyword_set(maximize) then f = -fnorm else f = fnorm print, iter, f, format='("Iter ",I6," FUNCTION = ",G20.8)' if n_elements(fmt) GT 0 then begin print, x, format=fmt endif else begin n = n_elements(x) ii = lindgen(n) p = ' P('+strtrim(ii,2)+') = '+string(x,format='(G)') if keyword_set(dprint) then begin p1 = strarr(n) wh = where(pfixed EQ 0, ct) if ct GT 0 AND n_elements(df) GE ct then begin if keyword_set(maximize) then df1 = -df else df1 = df p1(wh) = string(df1, format='(G)') endif wh = where(pfixed EQ 1, ct) if ct GT 0 then $ p1(wh) = ' (FIXED)' p = p + ' : DF/DP('+strtrim(ii,2)+') = '+p1 endif print, p, format='(A)' endelse return end ; SUBROUTINE TNBC (IERROR, N, X, F, G, W, LW, SFUN, LOW, UP, IPIVOT) ; IMPLICIT DOUBLE PRECISION (A-H,O-Z) ; INTEGER IERROR, N, LW, IPIVOT(N) ; DOUBLE PRECISION X(N), G(N), F, W(LW), LOW(N), UP(N) ; ; THIS ROUTINE SOLVES THE OPTIMIZATION PROBLEM ; ; MINIMIZE F(X) ; X ; SUBJECT TO LOW <= X <= UP ; ; WHERE X IS A VECTOR OF N REAL VARIABLES. THE METHOD USED IS ; A TRUNCATED-NEWTON ALGORITHM (SEE "NEWTON-TYPE MINIMIZATION VIA ; THE LANCZOS ALGORITHM" BY S.G. NASH (TECHNICAL REPORT 378, MATH. ; THE LANCZOS METHOD" BY S.G. NASH (SIAM J. NUMER. ANAL. 21 (1984), ; PP. 770-778). THIS ALGORITHM FINDS A LOCAL MINIMUM OF F(X). IT DOES ; NOT ASSUME THAT THE FUNCTION F IS CONVEX (AND SO CANNOT GUARANTEE A ; GLOBAL SOLUTION), BUT DOES ASSUME THAT THE FUNCTION IS BOUNDED BELOW. ; IT CAN SOLVE PROBLEMS HAVING ANY NUMBER OF VARIABLES, BUT IT IS ; ESPECIALLY USEFUL WHEN THE NUMBER OF VARIABLES (N) IS LARGE. ; ; SUBROUTINE PARAMETERS: ; ; IERROR - (INTEGER) ERROR CODE ; ( 0 => NORMAL RETURN ; ( 2 => MORE THAN MAXFUN EVALUATIONS ; ( 3 => LINE SEARCH FAILED TO FIND LOWER ; ( POINT (MAY NOT BE SERIOUS) ; (-1 => ERROR IN INPUT PARAMETERS ; N - (INTEGER) NUMBER OF VARIABLES ; X - (REAL*8) VECTOR OF LENGTH AT LEAST N; ON INPUT, AN INITIAL ; ESTIMATE OF THE SOLUTION; ON OUTPUT, THE COMPUTED SOLUTION. ; G - (REAL*8) VECTOR OF LENGTH AT LEAST N; ON OUTPUT, THE FINAL ; VALUE OF THE GRADIENT ; F - (REAL*8) ON INPUT, A ROUGH ESTIMATE OF THE VALUE OF THE ; OBJECTIVE FUNCTION AT THE SOLUTION; ON OUTPUT, THE VALUE ; OF THE OBJECTIVE FUNCTION AT THE SOLUTION ; W - (REAL*8) WORK VECTOR OF LENGTH AT LEAST 14*N ; LW - (INTEGER) THE DECLARED DIMENSION OF W ; SFUN - A USER-SPECIFIED SUBROUTINE THAT COMPUTES THE FUNCTION ; AND GRADIENT OF THE OBJECTIVE FUNCTION. IT MUST HAVE ; THE CALLING SEQUENCE ; SUBROUTINE SFUN (N, X, F, G) ; INTEGER N ; DOUBLE PRECISION X(N), G(N), F ; LOW, UP - (REAL*8) VECTORS OF LENGTH AT LEAST N CONTAINING ; THE LOWER AND UPPER BOUNDS ON THE VARIABLES. IF ; THERE ARE NO BOUNDS ON A PARTICULAR VARIABLE, SET ; THE BOUNDS TO -1.D38 AND 1.D38, RESPECTIVELY. ; IPIVOT - (INTEGER) WORK VECTOR OF LENGTH AT LEAST N, USED ; TO RECORD WHICH VARIABLES ARE AT THEIR BOUNDS. ; ; THIS IS AN EASY-TO-USE DRIVER FOR THE MAIN OPTIMIZATION ROUTINE ; LMQNBC. MORE EXPERIENCED USERS WHO WISH TO CUSTOMIZE PERFORMANCE ; OF THIS ALGORITHM SHOULD CALL LMQBC DIRECTLY. ; ;---------------------------------------------------------------------- ; THIS ROUTINE SETS UP ALL THE PARAMETERS FOR THE TRUNCATED-NEWTON ; ALGORITHM. THE PARAMETERS ARE: ; ; ETA - SEVERITY OF THE LINESEARCH ; MAXFUN - MAXIMUM ALLOWABLE NUMBER OF FUNCTION EVALUATIONS ; XTOL - DESIRED ACCURACY FOR THE SOLUTION X* ; STEPMX - MAXIMUM ALLOWABLE STEP IN THE LINESEARCH ; ACCRCY - ACCURACY OF COMPUTED FUNCTION VALUES ; MSGLVL - CONTROLS QUANTITY OF PRINTED OUTPUT ; 0 = NONE, 1 = ONE LINE PER MAJOR ITERATION. ; MAXIT - MAXIMUM NUMBER OF INNER ITERATIONS PER STEP ; function tnmin, fcn, xall, fguess=fguess, functargs=fcnargs, parinfo=parinfo, $ epsrel=epsrel0, epsabs=epsabs0, fastnorm=fastnorm, $ nfev=nfev, maxiter=maxiter0, maxnfev=maxfun0, maximize=fmax, $ errmsg=errmsg, nprint=nprint, status=status, nocatch=nocatch, $ iterproc=iterproc, iterargs=iterargs, niter=niter,quiet=quiet,$ autoderivative=autoderiv, iterderiv=iterderiv, bestmin=f if n_elements(nprint) EQ 0 then nprint = 1 if n_elements(iterproc) EQ 0 then iterproc = 'TNMIN_DEFITER' if n_elements(autoderiv) EQ 0 then autoderiv = 0 if n_elements(msglvl) EQ 0 then msglvl = 0 if n_params() EQ 0 then begin message, "USAGE: PARMS = TNMIN('MYFUNCT', START_PARAMS, ... )", /info return, !values.d_nan endif iterd = keyword_set(iterderiv) maximize = keyword_set(fmax) status = 0L nfev = 0L errmsg = '' catch_msg = 'in TNMIN' common tnmin_config, tnconfig tnconfig = {fastnorm: keyword_set(fastnorm), proc: 0, nfev: 0L, $ autoderiv: keyword_set(autoderiv), max: maximize} ;; Handle error conditions gracefully if NOT keyword_set(nocatch) then begin catch, catcherror if catcherror NE 0 then begin catch, /cancel 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 status EQ 0 then status = -18 return, !values.d_nan endif 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 X or PARINFO' goto, TERMINATE endif ;; Be sure that PARINFO is of the right type if n_elements(parinfo) GT 0 then begin 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 X 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 tnmin_parinfo, parinfo, tagnames, 'VALUE', xall, status=stx if stx EQ 0 then begin errmsg = 'ERROR: either X 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 npar = n_elements(xall) zero = xall(0) * 0. one = zero + 1 ten = zero + 10 twothird = (zero+2)/(zero+3) quarter = zero + 0.25 half = zero + 0.5 ;; Extract machine parameters sz = size(xall) tp = sz(sz(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 xall = float(xall) sz = size(xall) endif isdouble = (sz(sz(0)+1) EQ 5) common tnmin_machar, machvals tnmin_setmachar, double=isdouble MCHPR1 = machvals.machep ;; TIED parameters? tnmin_parinfo, parinfo, tagnames, 'TIED', ptied, default='', n=npar ptied = strtrim(ptied, 2) wh = where(ptied NE '', qanytied) qanytied = qanytied GT 0 tnconfig = create_struct(tnconfig, 'QANYTIED', qanytied, 'PTIED', ptied) ;; FIXED parameters ? tnmin_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 derivative tnmin_parinfo, parinfo, tagnames, 'STEP', step, default=zero, n=npar tnmin_parinfo, parinfo, tagnames, 'RELSTEP', dstep, default=zero, n=npar tnmin_parinfo, parinfo, tagnames, 'TNSIDE', dside, default=2, n=npar ;; Maximum and minimum steps allowed to be taken in one iteration tnmin_parinfo, parinfo, tagnames, 'TNMAXSTEP', maxstep, default=zero, n=npar tnmin_parinfo, parinfo, tagnames, 'TNMINSTEP', minstep, default=zero, n=npar qmin = minstep * 0 ;; Disable 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: TNMINSTEP is greater than TNMAXSTEP' goto, TERMINATE endif wh = where(qmin AND qmax, ct) qminmax = ct GT 0 ;; Finish up the free parameters ifree = where(pfixed NE 1, ct) if ct EQ 0 then begin errmsg = 'ERROR: no free parameters' goto, TERMINATE 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 ;; LIMITED parameters ? tnmin_parinfo, parinfo, tagnames, 'LIMITED', limited, status=st1 tnmin_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(n_elements(ifree)) ulim = x * 0. qllim = qulim llim = x * 0. qanylim = 0 endelse tnconfig = create_struct(tnconfig, $ 'PFIXED', pfixed, 'IFREE', ifree, $ 'STEP', step, 'DSTEP', dstep, 'DSIDE', dside, $ 'ULIMITED', qulim, 'ULIMIT', ulim) common tnmin_fcnargs, tnfcnargs tnfcnargs = 0 & dummy = temporary(tnfcnargs) if n_elements(fcnargs) GT 0 then tnfcnargs = fcnargs ;; SET UP CUSTOMIZING PARAMETERS ;; ETA - SEVERITY OF THE LINESEARCH ;; MAXFUN - MAXIMUM ALLOWABLE NUMBER OF FUNCTION EVALUATIONS ;; XTOL - DESIRED ACCURACY FOR THE SOLUTION X* ;; STEPMX - MAXIMUM ALLOWABLE STEP IN THE LINESEARCH ;; ACCRCY - ACCURACY OF COMPUTED FUNCTION VALUES ;; MSGLVL - DETERMINES QUANTITY OF PRINTED OUTPUT ;; 0 = NONE, 1 = ONE LINE PER MAJOR ITERATION. ;; MAXIT - MAXIMUM NUMBER OF INNER ITERATIONS PER STEP n = n_elements(x) if n_elements(maxit) EQ 0 then begin maxit = (n/2) < 50 > 2 ;; XXX diff than TN.F endif if n_elements(maxfun0) EQ 0 then $ maxfun = 0L $ else $ maxfun = floor(maxfun0(0)) > 1 ; maxfun = 150L*n ; if keyword_set(autoderiv) then $ ; maxfun = maxfun*(1L + round(total(abs(dside)> 1 < 2))) eta = zero + 0.25 stepmx = zero + 10 if n_elements(maxiter0) EQ 0 then $ maxiter = 200L $ else $ maxiter = floor(maxiter0(0)) > 1 g = replicate(x(0)* 0., n) ;; call minimizer ; ; THIS ROUTINE IS A BOUNDS-CONSTRAINED TRUNCATED-NEWTON METHOD. ; THE TRUNCATED-NEWTON METHOD IS PRECONDITIONED BY A LIMITED-MEMORY ; QUASI-NEWTON METHOD (THIS PRECONDITIONING STRATEGY IS DEVELOPED ; IN THIS ROUTINE) WITH A FURTHER DIAGONAL SCALING (SEE ROUTINE NDIA3). ; FOR FURTHER DETAILS ON THE PARAMETERS, SEE ROUTINE TNBC. ; ; ; initialize variables ; common tnmin_work, lsk, lyk, ldiagb, lsr, lyr ; I/O I/O I/O I/O I/O lsk = 0 & lyk = 0 & ldiagb = 0 & lsr = 0 & lyr = 0 zero = x(0)* 0. one = zero + 1 if n_elements(fguess) EQ 0 then fguess = one if maximize then f = -fguess else f = fguess conv = 0 & lreset = 0 & upd1 = 0 & newcon = 0 gsk = zero & yksk = zero & gtp = zero & gtpnew = zero & yrsr = zero upd1 = 1 ireset = 0L nmodif = 0L nlincg = 0L fstop = f conv = 0 nm1 = n - 1 ;; From CHKUCP ; ; CHECKS PARAMETERS AND SETS CONSTANTS WHICH ARE COMMON TO BOTH ; DERIVATIVE AND NON-DERIVATIVE ALGORITHMS ; EPSMCH = MCHPR1 SMALL = EPSMCH*EPSMCH TINY = SMALL NWHY = -1L ; ; SET CONSTANTS FOR LATER ; ;; Some of these constants have been moved around for clarity (!) if n_elements(epsrel0) EQ 0 then epsrel = 100*MCHPR1 $ else epsrel = epsrel0(0)+0. if n_elements(epsabs0) EQ 0 then epsabs = zero $ else epsabs = abs(epsabs0(0))+0. ACCRCY = epsrel XTOL = sqrt(ACCRCY) RTEPS = SQRT(EPSMCH) RTOL = XTOL IF (ABS(RTOL) LT ACCRCY) THEN RTOL = 10. *RTEPS FTOL2 = 0 FTOL1 = RTOL^2 + EPSMCH ;; For func chg convergence test (U1a) if epsabs NE 0 then $ FTOL2 = EPSABS + EPSMCH ;; For absolute func convergence test (U1b) PTOL = RTOL + RTEPS ;; For parm chg convergence test (U2) GTOL1 = ACCRCY^TWOTHIRD ;; For gradient convergence test (U3, squared) GTOL2 = (1D-2*XTOL)^2 ;; For gradient convergence test (U4, squared) ; ; CHECK FOR ERRORS IN THE INPUT PARAMETERS ; IF (ETA LT 0.D0 OR STEPMX LT RTOL) THEN BEGIN errmsg = 'ERROR: input keywords are inconsistent' goto, TERMINATE endif ;; Check input parameters for errors if (n LE 0) OR (xtol LE 0) OR (maxit LE 0) then begin errmsg = 'ERROR: input keywords are inconsistent' goto, TERMINATE endif NWHY = 0L XNORM = TNMIN_ENORM(X) ALPHA = zero TEST = zero ; From SETUCR ; ; CHECK INPUT PARAMETERS, COMPUTE THE INITIAL FUNCTION VALUE, SET ; CONSTANTS FOR THE SUBSEQUENT MINIMIZATION ; fm = f ; ; COMPUTE THE INITIAL FUNCTION VALUE ; catch_msg = 'calling TNMIN_CALL' fnew = tnmin_call(fcn, x, g, fullparam_=xnew) ; ; SET CONSTANTS FOR LATER ; NITER = 0L OLDF = FNEW GTG = TOTAL(G*G) common tnmin_error, tnerr if nprint GT 0 AND iterproc NE '' then begin iflag = 0L if (niter-1) MOD nprint EQ 0 then begin catch_msg = 'calling '+iterproc tnerr = 0 call_procedure, iterproc, fcn, xnew, niter, fnew, $ FUNCTARGS=fcnargs, parinfo=parinfo, quiet=quiet, $ dprint=iterd, deriv=g, pfixed=pfixed, maximize=maximize, $ _EXTRA=iterargs iflag = tnerr if iflag LT 0 then begin errmsg = 'WARNING: premature termination by "'+iterproc+'"' nwhy = 4L goto, CLEANUP endif endif endif fold = fnew flast = fnew ; From LMQNBC ; ; TEST THE LAGRANGE MULTIPLIERS TO SEE IF THEY ARE NON-NEGATIVE. ; BECAUSE THE CONSTRAINTS ARE ONLY LOWER BOUNDS, THE COMPONENTS ; OF THE GRADIENT CORRESPONDING TO THE ACTIVE CONSTRAINTS ARE THE ; LAGRANGE MULTIPLIERS. AFTERWORDS, THE PROJECTED GRADIENT IS FORMED. ; catch_msg = 'zeroing derivatives of pegged parameters' lmask = qllim AND (x EQ llim) AND (g GE 0) umask = qulim AND (x EQ ulim) AND (g LE 0) whlpeg = where(lmask, nlpeg) whupeg = where(umask, nupeg) tnmin_fix, whlpeg, whupeg, g GTG = TOTAL(G*G) ; ; CHECK IF THE INITIAL POINT IS A LOCAL MINIMUM. ; FTEST = ONE + ABS(FNEW) IF (GTG LT GTOL2*FTEST*FTEST) THEN GOTO, CLEANUP ; ; SET INITIAL VALUES TO OTHER PARAMETERS ; ICYCLE = NM1 GNORM = SQRT(GTG) DIFNEW = ZERO EPSRED = HALF/TEN FKEEP = FNEW ; ; SET THE DIAGONAL OF THE APPROXIMATE HESSIAN TO UNITY. ; LDIAGB = replicate(one, n) ; ; ..................START OF MAIN ITERATIVE LOOP.......... ; ; COMPUTE THE NEW SEARCH DIRECTION ; catch_msg = 'calling TNMIN_MODLNP' tnmin_modlnp, lpk, lgv, lz1, lv, ldiagb, lemat, $ x, g, lzk, n, niter, maxit, nmodif, nlincg, upd1, yksk, $ gsk, yrsr, lreset, fcn, whlpeg, whupeg, accrcy, gtpnew, gnorm, xnorm, $ xnew ITER_LOOP: catch_msg = 'computing step length' LOLDG = G PNORM = tnmin_enorm(LPK) OLDF = FNEW OLDGTP = GTPNEW ; ; PREPARE TO COMPUTE THE STEP LENGTH ; PE = PNORM + EPSMCH ; ; COMPUTE THE ABSOLUTE AND RELATIVE TOLERANCES FOR THE LINEAR SEARCH ; RELTOL = RTEPS*(XNORM + ONE)/PE ABSTOL = - EPSMCH*FTEST/(OLDGTP - EPSMCH) ; ; COMPUTE THE SMALLEST ALLOWABLE SPACING BETWEEN POINTS IN ; THE LINEAR SEARCH ; TNYTOL = EPSMCH*(XNORM + ONE)/PE ;; From STPMAX SPE = STEPMX/PE mmask = (NOT lmask AND NOT umask) wh = where(mmask AND (lpk GT 0) AND qulim AND (ulim - x LT spe*lpk), ct) if ct GT 0 then begin spe = min( (ulim(wh)-x(wh)) / lpk(wh)) endif wh = where(mmask AND (lpk LT 0) AND qllim AND (llim - x GT spe*lpk), ct) if ct GT 0 then begin spe = min( (llim(wh)-x(wh)) / lpk(wh)) endif ;; From LMQNBC ; ; SET THE INITIAL STEP LENGTH. ; ALPHA = TNMIN_STEP1(FNEW,FM,OLDGTP,SPE, epsmch) ; ; PERFORM THE LINEAR SEARCH ; catch_msg = 'performing linear search' tnmin_linder, n, fcn, small, epsmch, reltol, abstol, tnytol, $ eta, zero, spe, lpk, oldgtp, x, fnew, alpha, g, nwhy, xnew NEWCON = 0 IF (ABS(ALPHA-SPE) GT 1.D1*EPSMCH) EQ 0 THEN BEGIN NEWCON = 1 NWHY = 0L ;; From MODZ mmask = (NOT lmask AND NOT umask) wh = where(mmask AND (lpk LT 0) AND qllim $ AND (x-llim LE 10*epsmch*(abs(llim)+one)),ct) if ct GT 0 then begin flast = fnew lmask(wh) = 1 x(wh) = llim(wh) whlpeg = where(lmask, nlpeg) endif wh = where(mmask AND (lpk GT 0) AND qulim $ AND (ulim-x LE 10*epsmch*(abs(ulim)+one)),ct) if ct GT 0 then begin flast = fnew umask(wh) = 1 x(wh) = ulim(wh) whupeg = where(umask, nupeg) endif xnew(ifree) = x ;; From LMQNBC FLAST = FNEW endif FOLD = FNEW NITER = NITER + 1 ; ; IF REQUIRED, PRINT THE DETAILS OF THIS ITERATION ; if nprint GT 0 AND iterproc NE '' then begin iflag = 0L xx = xnew xx(ifree) = x if (niter-1) MOD nprint EQ 0 then begin catch_msg = 'calling '+iterproc tnerr = 0 call_procedure, iterproc, fcn, xx, niter, fnew, $ FUNCTARGS=fcnargs, parinfo=parinfo, quiet=quiet, $ dprint=iterd, deriv=g, pfixed=pfixed, maximize=maximize, $ _EXTRA=iterargs iflag = tnerr if iflag LT 0 then begin errmsg = 'WARNING: premature termination by "'+iterproc+'"' nwhy = 4L goto, CLEANUP endif endif endif catch_msg = 'testing for convergence' IF (NWHY LT 0) THEN BEGIN NWHY = -2L goto, CLEANUP ENDIF IF (NWHY NE 0 AND NWHY NE 2) THEN BEGIN ;; THE LINEAR SEARCH HAS FAILED TO FIND A LOWER POINT NWHY = 3L goto, CLEANUP endif if nwhy GT 1 then begin fnew = tnmin_call(fcn, x, g, fullparam_=xnew) endif wh = where(finite(x) EQ 0, ct) if ct GT 0 OR finite(fnew) EQ 0 then begin nwhy = -3L goto, CLEANUP endif ; ; TERMINATE IF MORE THAN MAXFUN EVALUATIONS HAVE BEEN MADE ; NWHY = 2L if maxfun GT 0 AND tnconfig.nfev GT maxfun then goto, CLEANUP if niter GT maxiter then goto, CLEANUP NWHY = 0L ; ; SET UP PARAMETERS USED IN CONVERGENCE AND RESETTING TESTS ; DIFOLD = DIFNEW DIFNEW = OLDF - FNEW ; ; IF THIS IS THE FIRST ITERATION OF A NEW CYCLE, COMPUTE THE ; PERCENTAGE REDUCTION FACTOR FOR THE RESETTING TEST. ; IF (ICYCLE EQ 1) THEN BEGIN IF (DIFNEW GT 2.D0*DIFOLD) THEN EPSRED = EPSRED + EPSRED IF (DIFNEW LT 5.0D-1*DIFOLD) THEN EPSRED = HALF*EPSRED ENDIF LGV = G tnmin_fix, whlpeg, whupeg, lgv GTG = TOTAL(LGV*LGV) GNORM = SQRT(GTG) FTEST = ONE + ABS(FNEW) XNORM = tnmin_enorm(X) ;; From CNVTST LTEST = (FLAST - FNEW) LE (-5.D-1*GTPNEW) wh = where((lmask AND g LT 0) OR (umask AND g GT 0), ct) if ct GT 0 then begin conv = 0 if NOT ltest then begin mx = max(abs(g(wh)), wh2) lmask(wh(wh2)) = 0 & umask(wh(wh2)) = 0 FLAST = FNEW goto, CNVTST_DONE endif endif ;; Gill Murray and Wright tests are listed to the right. ;; Modifications due to absolute function value test are done here. fconv = abs(DIFNEW) LT FTOL1*FTEST ;; U1a if ftol2 EQ 0 then begin pconv = ALPHA*PNORM LT PTOL*(ONE + XNORM) ;; U2 gconv = GTG LT GTOL1*FTEST*FTEST ;; U3 endif else begin ;; Absolute tolerance implies a loser constraint on parameters fconv = fconv OR (abs(difnew) LT ftol2) ;; U1b acc2 = (FTOL2/FTEST + EPSMCH) pconv = ALPHA*PNORM LT SQRT(acc2)*(ONE + XNORM) ;; U2 gconv = GTG LT (acc2^twothird)*FTEST*FTEST ;; U3 endelse IF ((PCONV AND FCONV AND GCONV) $ ;; U1 + U2 + U3 OR (GTG LT GTOL2*FTEST*FTEST)) THEN BEGIN ;; U4 CONV = 1 ENDIF ELSE BEGIN ;; Convergence failed CONV = 0 ENDELSE ; ; FOR DETAILS, SEE GILL, MURRAY, AND WRIGHT (1981, P. 308) AND ; FLETCHER (1981, P. 116). THE MULTIPLIER TESTS (HERE, TESTING ; THE SIGN OF THE COMPONENTS OF THE GRADIENT) MAY STILL NEED TO ; MODIFIED TO INCORPORATE TOLERANCES FOR ZERO. ; CNVTST_DONE: ;; From LMQNBC IF (CONV) THEN GOTO, CLEANUP tnmin_fix, whlpeg, whupeg, g ; ; COMPUTE THE CHANGE IN THE ITERATES AND THE CORRESPONDING CHANGE ; IN THE GRADIENTS ; IF NEWCON EQ 0 THEN BEGIN LYK = G - LOLDG LSK = ALPHA*LPK ; ; SET UP PARAMETERS USED IN UPDATING THE PRECONDITIONING STRATEGY. ; YKSK = TOTAL(LYK*LSK) LRESET = 0 IF (ICYCLE EQ NM1 OR DIFNEW LT EPSRED*(FKEEP-FNEW)) THEN LRESET = 1 IF (LRESET EQ 0) THEN BEGIN YRSR = TOTAL(LYR*LSR) IF (YRSR LE ZERO) THEN LRESET = 1 ENDIF UPD1 = 0 ENDIF ; ; COMPUTE THE NEW SEARCH DIRECTION ; ;; TNMIN_90: catch_msg = 'calling TNMIN_MODLNP' tnmin_modlnp, lpk, lgv, lz1, lv, ldiagb, lemat, $ x, g, lzk, n, niter, maxit, nmodif, nlincg, upd1, yksk, $ gsk, yrsr, lreset, fcn, whlpeg, whupeg, accrcy, gtpnew, gnorm, xnorm, $ xnew IF (NEWCON) THEN GOTO, ITER_LOOP ; IF (NOT LRESET) OR ICYCLE EQ 1 AND n_elements(LSR) GT 0 THEN BEGIN ;; For testing IF (LRESET EQ 0) THEN BEGIN ; ; COMPUTE THE ACCUMULATED STEP AND ITS CORRESPONDING ; GRADIENT DIFFERENCE. ; LSR = LSR + LSK LYR = LYR + LYK ICYCLE = ICYCLE + 1 goto, ITER_LOOP ENDIF ; ; RESET ; ;; TNMIN_110: IRESET = IRESET + 1 ; ; INITIALIZE THE SUM OF ALL THE CHANGES IN X. ; LSR = LSK LYR = LYK FKEEP = FNEW ICYCLE = 1L goto, ITER_LOOP ; ; ...............END OF MAIN ITERATION....................... ; CLEANUP: nfev = tnconfig.nfev tnfcnargs = 0 catch, /cancel case NWHY of -3: begin ;; INDEFINITE VALUE status = -16L if errmsg EQ '' then $ errmsg = ('ERROR: parameter or function value(s) have become '+$ 'infinite; check model function for over- '+$ 'and underflow') return, !values.d_nan end -2: begin ;; INTERNAL ERROR IN LINE SEARCH status = -18L if errmsg EQ '' then $ errmsg = 'ERROR: Fatal error during line search' return, !values.d_nan end -1: begin TERMINATE: ;; FATAL TERMINATION status = 0L if errmsg EQ '' then errmsg = 'ERROR: Invalid inputs' return, !values.d_nan end 0: begin CONVERGED: status = 1L end 2: begin ;; MAXIMUM NUMBER of FUNC EVALS or ITERATIONS REACHED if maxfun GT 0 AND nfev GT maxfun then begin status = -17L if errmsg EQ '' then $ errmsg = ('WARNING: no convergence within maximum '+$ 'number of function calls') endif else begin status = 5L if errmsg EQ '' then $ errmsg = ('WARNING: no convergence within maximum '+$ 'number of iterations') endelse FNEW = OLDF end 3: begin status = -18L if errmsg EQ '' then errmsg = 'ERROR: Line search failed to converge' end 4: begin ;; ABNORMAL TERMINATION BY USER ROUTINE status = iflag end endcase ;; Successful return F = FNEW xnew(ifree) = x if keyword_set(tnconfig.qanytied) then tnmin_tie, xnew, tnconfig.ptied return, xnew end