; function GOLDEN - golden-section optimization
;
; This function finds the independent variable which optimizes the given function
;
; Input:
;   ax,bx,cx - Bracket. Values of the independent variable such that
;                1) bx is between ax and cx
;                2) f(bx) is more optimum both f(ax) and f(cx)
;   f -        User function. String with the name of the function to call
;   tol -      fractional precision. Default is 1e-3 if not specified
;   /max -     If set, maximizes the function, otherwise minimizes it
;   _extra -   All other parameters are passed to the user function
; Output:
;   xopt -     value of independent variable which optimizes the function
; Return:
;   The value of the function at its optimum point
;
;Given a function f, and given a bracketing triplet of abscissas ax, bx, cx (such that bx is
;between ax and cx, and f(bx) is more optimum than both f(ax) and f(cx)), this routine performs a
;golden section search for the optimum, isolating it to a fractional precision of about tol. The
;abscissa of the minimum is returned as xopt, and the minimum function value is returned as
;golden, the returned function value.

function golden,ax,bx,cx,f,tol,xopt=xopt,max=max,_extra=extra
  if n_elements(tol) eq 0 then tol=1d-3
  R=(sqrt(5d)-1)/2 ;The golden ratios.
  C=1d -R
  ;#define SHFT2(a,b,c) (a)=(b);(b)=(c);
  ;#define SHFT3(a,b,c,d) (a)=(b);(b)=(c);(c)=(d);

  ;From Numerical Recipes, section 10.1
  f1=0d
  f2=0d

  x0=ax         ; At any given time we will keep track of four points,
  x3=cx         ; x0,x1,x2,x3.
  if (abs(cx-bx) gt abs(bx-ax)) then begin ;Make x0 to x1 the smaller segment,
    x1=bx
    x2=bx+C*(cx-bx); and fill in the new point to be tried.
  end else begin
    x2=bx;
    x1=bx-C*(bx-ax);
  end

  f1=call_function(f,x1,_extra=extra) ; The initial function evaluations. Note that
                                      ; we never need to evaluate the function
                                      ; at the original endpoints.
  if keyword_set(max) then f1=-f1
  f2=call_function(f,x2,_extra=extra);
  if keyword_set(max) then f2=-f2
  while (abs(x3-x0) gt tol*(abs(x1)+abs(x2))) do begin
    if (f2 lt f1) then begin ;One possible outcome,
      ;SHFT3(x0,x1,x2,R*x1+C*x3) its housekeeping,
      x0=x1
      x1=x2
      x2=R*x1+C*x3
      ;SHFT2(f1,f2,(*f)(x2)) and a new function evaluation.
      f1=f2
      f2=call_function(f,x2,_extra=extra)
      if keyword_set(max) then f2=-f2
    end else begin ; The other outcome,
      ;SHFT3(x3,x2,x1,R*x2+C*x0)
      x3=x2
      x2=x1
      x1=R*x2+C*x0
      ;SHFT2(f2,f1,(*f)(x1)) and its new function evaluation.
      f2=f1
      f1=call_function(f,x1,_extra=extra)
      if keyword_set(max) then f1=-f1
    end
  end ; Back to see if we are done.
  if (f1 lt f2) then begin ; We are done. Output the best of the two
    xopt=x1                ; current values.
    if keyword_set(max) then f1=-f1
    return,f1
  end else begin
    xopt=x2;
    if keyword_set(max) then f2=-f2
    return,f2
  end
end