; Procedure minF_parabol_D, ; first, a utility function which gets derivative in 1-D: ;------------------------------------------------------------------------------ function call_func_deriv, func_name, x, deriv, POINT_NDIM=pn, DIRECTION=dirn f = call_function( func_name, pn + x * dirn, grad ) deriv = total( [grad * dirn] ) return, f end ;------------------------------------------------------------------------------ pro minF_parabol_D, xa,xb,xc, xmin, fmin, FUNC_NAME=func_name, $ MAX_ITERATIONS=maxit, $ TOLERANCE=TOL, $ POINT_NDIM=pn, DIRECTION=dirn ;+ ; NAME: ; MINF_PARABOL_D ; PURPOSE: ; Minimize a function using a modified Brent's method with derivatives ; EXPLANATION: ; Based on the procedure DBRENT in Numerical Recipes by Press et al. ; Finds a local minimum of a 1-D function up to specified tolerance, ; using the first derivative of function in the algorithm. ; This routine assumes that the function has a minimum nearby. ; (recommend first calling minF_bracket, xa,xb,xc, to bracket minimum). ; Routine can also be applied to a scalar function of many variables, ; for such case the local minimum in a specified direction is found, ; This routine is called by minF_conj_grad, to locate minimum in the ; direction of the conjugate gradient of function of many variables. ; ; CALLING EXAMPLES: ; minF_parabol_D, xa,xb,xc, xmin, fmin, FUNC_NAME="name" ;for 1-D func. ; or: ; minF_parabol_D, xa,xb,xc, xmin, fmin, FUNC="name", $ ; POINT=[0,1,1], $ ; DIRECTION=[2,1,1] ;for 3-D func. ; INPUTS: ; xa,xb,xc = scalars, 3 points which bracket location of minimum, ; that is, f(xb) < f(xa) and f(xb) < f(xc), so minimum exists. ; When working with function of N variables ; (xa,xb,xc) are then relative distances from POINT_NDIM, ; in the direction specified by keyword DIRECTION, ; with scale factor given by magnitude of DIRECTION. ; KEYWORDS: ; FUNC_NAME = function name (string) ; Calling mechanism should be: F = func_name( px, gradient ) ; where: ; px = scalar or vector of independent variables, input. ; F = scalar value of function at px. ; gradient = derivative of function, a scalar if 1-D, ; a gradient vector if N-D, ; (should only be computed if arg. is present). ; ; POINT_NDIM = when working with function of N variables, ; use this keyword to specify the starting point in N-dim space. ; Default = 0, which assumes function is 1-D. ; DIRECTION = when working with function of N variables, ; use this keyword to specify the direction in N-dim space ; along which to bracket the local minimum, (default=1 for 1-D). ; (xa, xb, xc, x_min are then relative distances from POINT_NDIM) ; MAX_ITER = maximum allowed number iterations, default=100. ; TOLERANCE = desired accuracy of minimum location, default=sqrt(1.e-7). ; ; OUTPUTS: ; xmin = estimated location of minimum. ; When working with function of N variables, ; xmin is the relative distance from POINT_NDIM, ; in the direction specified by keyword DIRECTION, ; with scale factor given by magnitude of DIRECTION, ; so that min. Loc. Pmin = Point_Ndim + xmin * Direction. ; fmin = value of function at xmin (or Pmin). ; PROCEDURE: ; Brent's method to minimize a function by using parabolic interpolation ; and using first derivative of function, ; from Numerical Recipes (by Press, et al.), sec.10.3 (p.287), ; MODIFICATION HISTORY: ; Written, Frank Varosi NASA/GSFC 1992. ;- zeps = 1.e-7 ;machine epsilon, smallest addition. if N_elements( TOL ) NE 1 then TOL = sqrt( zeps ) if N_elements( maxit ) NE 1 then maxit = 100 if N_elements( pn ) LE 0 then begin pn = 0 dirn = 1 endif xLo = xa < xc xHi = xa > xc xmin = xb fmin = call_func_deriv( func_name, xmin, dx, POINT=pn, DIR=dirn ) xv = xmin & xw = xmin fv = fmin & fw = fmin dv = dx & dw = dx es = 0. for iter = 1,maxit do begin xm = (xLo + xHi)/2. TOL1 = TOL * abs(xmin) + zeps TOL2 = 2*TOL1 if ( abs( xmin - xm ) LE ( TOL2 - (xHi-xLo)/2. ) ) then return if (abs( es ) GT TOL1) then begin d1 = 2*(xHi-xLo) d2 = d1 if (dw NE dx) then d1 = (xw-xmin)*dx/(dx-dw) if (dv NE dx) then d2 = (xv-xmin)*dx/(dx-dv) u1 = xmin + d1 u2 = xmin + d2 ok1 = ((xLo-u1)*(u1-xHi) GT 0) AND (dx*d1 LE 0) ok2 = ((xLo-u2)*(u2-xHi) GT 0) AND (dx*d2 LE 0) olde = es es = ds if NOT (ok1 OR ok2) then goto,BISECT if (ok1 AND ok2) then begin if (abs( d1 ) LT abs( d2 )) then ds=d1 else ds=d2 endif else if (ok1) then ds=d1 else ds=d2 if (abs( ds ) LE abs( olde/2 )) then begin xu = xmin + ds if ((xu-xLo) LT TOL2) OR $ ((xHi-xu) LT TOL2) then $ ds = TOL1 * (1-2*((xm-xmin) LT 0)) goto,STEP endif endif BISECT: if (dx GE 0) then es = xLo-xmin else es = xHi-xmin ds = es/2 STEP: sign = 1 - 2*(ds LT 0) xu = xmin + sign * ( abs( ds ) > TOL1 ) fu = call_func_deriv( func_name, xu, du, POINT=pn, DIR=dirn ) if (fu GT fmin) AND (abs( ds ) LT TOL1) then return if (fu LE fmin) then begin if (xu GE xmin) then xLo=xmin else xHi=xmin xv = xw & fv = fw & dv = dw xw = xmin & fw = fmin & dw = dx xmin = xu & fmin = fu & dx = du endif else begin if (xu LT xmin) then xLo=xu else xHi=xu if (fu LE fw) OR (xw EQ xmin) then begin xv = xw & fv = fw & dv = dw xw = xu & fw = fu & dw = du endif else if (fu LE fv) OR (xv EQ xmin) $ OR (xv EQ xw) then begin xv = xu & fv = fu & dv = du endif endelse endfor message,"exceeded maximum number of iterations: "+strtrim(iter,2),/INFO return end