brent.pro (Chris Beaumont's IDL Library)

;+
; PURPOSE:
;  This function minimizes a function of one variable. The function
;  requires 3 input abscissas which bracket at least one local minimum
;  (see, e.g. bracket.pro)
;
;  The algorithm is adapted from Numerical Recipes, but is
;  substantially modified in an attempt to simplify the
;  procedure. This may make the program slower, but I don't
;  think it should be any less stable.
;
;  The program attempts to improve upon goldenmin by using parabolic
;  extrapolation in the vicinity of the minimum. If the function is
;  smooth, then this should quickly converge. If it seems that
;  parabolic extrapolation isn't behaving well, then it
;  defaults to a goldenmin style partitioning of (abc) into golden
;  sections. This ensures convergence.
;
;  The  biggest change from the routine in Numerical Recipes is the
;  test for when to use parabolic extrapolation. The default state is
;  to use parabolic extrapolation if the extrapolated minimum is
;  interior to (a,c). However, if the range of (a,c) decreases by less
;  than a factor of 1.5 in any iteration, then the algorithm is
;  'punished', and parabolic extrapolation is not used for a
;  round. The value of 1.5 is used because the goldenmin style of
;  partitioning shrinks (a,c) by about 1.67 per iteration. For smooth
;  functions, typical shrinkages seem to be >2 per iteration. They can
;  be as high as several thousand for nearly parabolic functions.
;
; CATEGORY:
;  Numerical Recipes
;
; CALLING SEQUENCE:
;  result = brent(func, ax, bx, cx, 
;                 fax, fbx, fcx, [tol tol, 
;                 fmin = fmin, /verbose, _extra = extra, /plot, 
;                 golden = golden)
;
; INPUTS:
;  func: The name of a user written function to minimize. The function
;        must have a calling sequence like result = func(x, _extra =
;        extra). It may declare extra keywords, which are supplied to
;        brent. It must return a scalar.
;
;   ax: The first point bracketing the minimum
;   bx: The second point bracketing the minimum
;   cx: The third point bracketing the minimum
;   fax: func(ax)
;   fbx: func(bx)
;   fcx: func(cx)
;
; KEYWORD PARAMETERS:
;  tol: The requested fractional precision of xmin. Defaults to 1d-3
;  verbose: Set to produce textual output
;  fmin: Set to a variable to hold f(xmin)
;  _extra: extra kewyords to pass along to calls of func
;  plot: Set to plot the points as it plugs along.
;  golden: If set, turn off parabolic interpolation (i.e. basically
;  golden search)
;
; OUTPUTS:
;  xmin: The approximate location of a minimum within (abc)
;
;- PSEUDOCODE:
;   xmin = brent(func,  // name of function
;                a,b,c, // points bracketing minimum
;                tol)   // fractional accuarcy of xmin
;  
;     step1 = step2 = 0
;     While (TRUE) do:
;       if (c - a) lt abs(b * tol) then return b
;       u  <-- chooseCandidate(a,b,c, step2). // new candidate min
;       fu <-- func(u)
;       temp <-- b               
;       updateBracket(a,b,c,u)   // narrow the bracket
;     endwhile
;   END
;
; SEE ALSO:
;  bracket, goldenmin
;
; MODIFICATION HISTORY:
;  June 2009: Adapted from NR by Chris Beaumont
;-
function brent, func, ax, bx, cx, $
                fax, fbx, fcx, $
                tol = tol, $
                verbose = verbose, $
                fmin = fmin, $
                _extra = extra, $
                plot = plot, $
                golden = golden

  compile_opt idl2
;- check inputs;
  if n_params() ne 7 then begin
     print, 'brent calling sequence:'
     print, ' xmin = brent(func, ax, bx, cx, fax, fbx, fcx, '
     print, '              [tol = tol, /verbose, fmin = fmin, '
     print, '               golden = golden, _extra = extra, /plot])'
     return, !values.f_nan
  endif
  
  if size(func, /type) ne 7 then $
     message, 'func must be a string'
  
  if (bx - ax) * (cx - bx) le 0 then  $
     message, 'b must be between a and c'
  
  if fbx ge fax || fbx ge fcx then $
     message, 'fbx must be < fax, fcx'
  
  if ~keyword_set(tol) then tol = 1d-3
  
  a = (ax lt cx) ? ax : cx
  c = (ax gt cx) ? ax : cx
  fa = (ax lt cx) ? fax : fcx
  fc = (ax gt cx) ? fax : fcx
  b = bx
  fb = double(fbx)
  
  shrinkage = 0D
  punishParabola = 0

  ITMAX = 100
  GOLD  = .3819660
  EPS   = 1d-7

  if keyword_set(plot) then begin
     loadct, 13, /silent
     plot, [a,b,c], [fa, fb, fc], psym = symcat(16), color = fsc_color('white')
  endif
 
  for i = 0, ITMAX - 1, 1 do begin
     if keyword_set(plot) then oplot, [a,b,c],[fa,fb,fc], sym = 8, $
                                     color = (5 * i) mod 255
     
     t1 = (tol * abs(b) + EPS)
     
     ;- output
     verbiage, string(i,a,b,c,format = '("iter: ", i3, "'+$
                      '   x:  ", f, ", ", f, ", ", f)'), 3, verbose
     verbiage, string(fa, fb, fc, $
                      format='("         f(x):  ", f, ", ", f, ", ", f)'), $
               3, verbose
     verbiage, string(abs(c - a) / (2 * t1), format='("       bracket / tol: ", e9.1)'), $
               3, verbose
     
     ;- test for convergence
     if abs(c - a) lt 2 * t1 then begin
        verbiage, 'Converged in '+strtrim(i,2)+' iterations', 2, verbose
        verbiage, string(b, fb, format='("xmin: ", f, " f(xmin): ", f)'), 2, verbose
        fmin = fb
        return, b
     endif
     
    assert, (b gt a) && (c gt b)

     ;- calculate two trial points
     u1 = (c - b) gt (b - a) ? $
          b + GOLD * (c - b) : $
          b - GOLD * (b - a)
     par = parabola(a, b, c, fa, fb, fc)
     u2 = - par[1] / (2 * par[2])

          ;- special case: u2 is very close to b
     ;- this is probably a good thing - b is very close to the right answer.
     if abs(u2 - b) lt t1 then begin
        ;- to the right of b?
        if u2 gt b then begin
           if (c - b) lt t1 then u2 = b - t1 else u2 = b + t1
        endif else begin
           if (b - a) lt t1 then u2 = b + t1 else u2 = b - t1
        endelse
     endif
     
     verbiage, string(u2, format='("   Parabola x:", f)'), 3, verbose
     

     ;- is u2 acceptable?
     if ~keyword_set(golden) && ~punishParabola && $
        (u2 gt a) && (u2 lt c) then begin
        verbiage, '     parabola accepted.', 3, verbose
        u = u2 
     endif else begin
        u = u1
        punishParabola = 0
     endelse

     ;- evaluate function at u and update bracket
     fu = call_function(func, u, _extra = extra)
     
     shrinkage = c - a
     if u gt b then begin
        ; buc bracket
        if fu lt fb then begin
           fa = fb
           fb = fu
           a = b
           b = u
        ; abu bracket
        endif else begin
           fc = fu
           c = u
        endelse
     endif else begin
      ;-aub bracket
        if (fu lt fb) then begin
           fc = fb
           fb = fu
           c = b
           b = u
        ;-ubc bracket
        endif else begin
           fa = fu
           a = u
        endelse
     endelse
     
     shrinkage = shrinkage / (c - a)
     verbiage, string(shrinkage, format='("    shrinkage:", e9.1)'), 3, verbose
     ;- golden shrinkage is about 1.618. 
     ;- punish parabola for not shrinking enough.
     if shrinkage lt 1.5 then punishParabola = 1
     assert, shrinkage gt 1
  endfor 

  ;- failure
  verbiage, 'Did not converge to a minimum within '+strtrim(MAXITER, 2)+' iterations', 1, verbose
  return, !values.f_nan
  
end