pro minF_parabolic, xa,xb,xc, xmin, fmin, FUNC_NAME=func_name, \$ MAX_ITERATIONS=maxit, \$ TOLERANCE=TOL, \$ POINT_NDIM=pn, DIRECTION=dirn ;+ ; NAME: ; MINF_PARABOLIC ; PURPOSE: ; Minimize a function using Brent's method with parabolic interpolation ; EXPLANATION: ; Find a local minimum of a 1-D function up to specified tolerance. ; 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_parabolic, xa,xb,xc, xmin, fmin, FUNC_NAME="name" ;for 1-D func. ; or: ; minF_parabolic, 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. ; INPUT KEYWORDS: ; FUNC_NAME = function name (string) ; Calling mechanism should be: F = func_name( px ) ; where: ; px = scalar or vector of independent variables, input. ; F = scalar value of function at px. ; ; 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. ; Based on function BRENT in Numerical Recipes in FORTRAN (Press et al. ; 1992), sec.10.2 (p. 397). ; MODIFICATION HISTORY: ; Written, Frank Varosi NASA/GSFC 1992. ; Converted to IDL V5.0 W. Landsman September 1997 ;- zeps = 1.e-7 ;machine epsilon, smallest addition. goldc = 1 - (sqrt(5)-1)/2 ;complement of golden mean. 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_function( func_name, pn + xmin * dirn ) xv = xmin & xw = xmin fv = fmin & fw = fmin es = 0. for iter = 1,maxit do begin goldstep = 1 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 r = (xmin-xw) * (fmin-fv) q = (xmin-xv) * (fmin-fw) p = (xmin-xv) * q + (xmin-xw) * r q = 2 * (q-r) if (q GT 0) then p = -p q = abs( q ) etemp = es es = ds if (p GT q*(xLo-xmin)) AND \$ (p LT q*(xHi-xmin)) AND \$ (abs( p ) LT abs( q*etemp/2 )) then begin ds = p/q xu = xmin + ds if (xu-xLo LT TOL2) OR (xHi-xu LT TOL2) then \$ ds = TOL1 * (1-2*((xm-xmin) LT 0)) goldstep = 0 endif endif if (goldstep) then begin if (xmin GE xm) then es = xLo-xmin else es = xHi-xmin ds = goldc * es endif xu = xmin + (1-2*(ds LT 0)) * ( abs( ds ) > TOL1 ) fu = call_function( func_name, pn + xu * dirn ) if (fu LE fmin) then begin if (xu GE xmin) then xLo=xmin else xHi=xmin xv = xw & fv = fw xw = xmin & fw = fmin xmin = xu & fmin = fu 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 xw = xu & fw = fu endif else if (fu LE fv) OR (xv EQ xmin) \$ OR (xv EQ xw) then begin xv = xu & fv = fu endif endelse endfor message,"exceeded maximum number of iterations: "+strtrim(iter,2),/INFO return end