function [xmin,fmin] = nr_brent(f, ax, bx, cx, tol)
%NR_BRENT - Brent's Method for Minimization.
%
%  [ xmin, fmin ] = brent(f, ax, bx, cx, tol)
%
%  Given a function f, and given a bracketing triplet of abscissas ax,
%  bx and cx, and fx(b) is less than both f(ax) and f(cx), this routine
%  isolates the minimum to a fractional precision of about TOL using
%  Brent's method. The abscissa of the minimum is returned in xmin, and
%  the minimum function value is returned in fmin.
%
%  References:
%
%  Brent, R.P. 1973, Algorithms for Minimization without Derivatives
%     (Englewood Cliffs, NJ: Prentice Hall), Chapter 5
%  Press, W.H. et al, 1992,  Numerical recipes in C: the art of scientific
%     computing (Cambridge University Press), Chapter 10.2
%
%  Description of variables:
%
%    a, b   the minimum is bracketed between a and b
%    x      point with the least function value found so far
%    w      point with the second least function value found so far
%    v      previous value of w
%    u      point at which the function was evaluated most recently
%    xm     midpoint between a and b (the function is not evaluated there)
%
%    Note: these points are not necessarily all distinct.
%
%    e      movement from best current value in the last iteration step
%    etemp  movement from the best current value in the second last
%           iteration step
%
%  General principles:
%
%    - The parabola is fitted trough the points x, v and w
%    - To be acceptable the parabolic step must
%        (i)  fall within the bounding interval (a,b)
%        (ii) imply a movement from the best current value x that is
%             *less* than half the movement of the *step before last*
%    - The code never evaluates the function less than a distance tol1
%      from a point already evaluated or from a known bracketing point
% 
% Written by:
% -- 
% John L. Weatherwax                2004-12-11 
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----


VERBOSE = 1;                % Print steps taken for the solution.
CGOLD   = (3-sqrt(5))/2;    % golden ratio
ITMAX   = 100;              % max number of iterations
ZEPS    = 1e-10;            % absolute error tolerance

% initialization
e = 0;
if ax < cx
  a = ax; b = cx;
else
  a = cx; b = ax;
end;
v = bx; w = v; x = w;
fx = feval(f,x); fv = fx; fw = fv;

for iter = 1:ITMAX,
  if( VERBOSE ) 
    fprintf(1, 'k=%4d, |a-b|=%e\n', iter, abs(a-b));
  end
  xm = 0.5*(a+b);
  tol1 = tol*abs(x) + ZEPS;
  tol2 = 2*tol1;
  % Stopping criterion: equivalent to: max(x-a, b-x) <= tol2
  if abs(x-xm) <= tol2-0.5*(b-a)
    xmin = x;
    fmin = fx;
    return
  end
  if abs(e) > tol1
    %
    % The second last move was sufficently large:
    % let's construct the parabolic fit
    %
    r = (x-w)*(fx-fv);
    q = (x-v)*(fx-fw);
    p = (x-v)*q - (x-w)*r;
    q = 2.0*(q-r);
    if q > 0, p = -p; end
    q = abs(q);
    etemp = e;
    e = d;
    if abs(p) >= abs(0.5*q*etemp) | p <= q*(a-x) | p >= q*(b-x)
      %
      % The parabolic fit did not meet the reqirements above:
      % use a golden section refinement instead.
      %
      % Print an explanation of what condition failed:
      %
      if( VERBOSE ) 
        if abs(p) >= abs(0.5*q*etemp)
          disp('abs(p) >= abs(0.5*q*etemp)')
        end
        if p <= q*(a-x)
          disp('p <= q*(a-x)')
        end
        if p >= q*(b-x)
          disp('p >= q*(b-x)')
        end
      end
      %
      % Set e in such a way, that a parabolic fit is possible in the
      % next iteration.
      %
      if x >= xm
        e = a-x;
      else
        e = b-x;
      end
      d = CGOLD*e;
      if( VERBOSE ) 
        disp('golden section step');
      end
    else
      %
      % the parabolic fit meets the above requirements: use iteration
      %
      d = p/q;
      u = x+d;
      if u-a < tol2 | b-u < tol2
        %
        % the parabolic fit is too close to the boundaries of the
        % bracketing interval: put new point tol1 apart from the
        % midpoint
        %
        d = tol1*sign(xm-x);
        if( VERBOSE ) disp('too close to interval boundaries'); end
      else
        if( VERBOSE ) disp('parabolic step'); end
      end
    end
  else
    %
    % The second last move was not sufficently large: use golden section
    % set e in such a way, that a parabolic fit is possible in the next 
    % iteration.
    %
    if x >= xm
      e = a-x;
    else
      e = b-x;
    end
    d = CGOLD*e;
    if( VERBOSE ) 
      disp('abs(e) > tol1');
      disp('golden section step');
    end
  end
  if abs(d) >= tol1
    u = x+d;
  else
    %
    % u is too close to x: put u farer apart
    %
    u = x + tol1*sign(d);
    if( VERBOSE ) disp('u is too close to x'); end;
  end
  fu = feval(f,u);    % finally evaluate the function
  if fu <= fx
    %
    % u is better than x: x becomes new boundary
    %
    if u >= x
      a = x;
    else
      b = x;
    end
    v = w; w = x; x = u;
    fv = fw; fw = fx; fx = fu;
  else
    %
    % x is better than u: u becomes new boundary
    %
    if u < x
      a = u;
    else
      b = u;
    end
    if fu <= fw | w == x
      v = w; w = u;
      fv = fw; fw = fu;
    elseif fu <= fv | v == x | v == w
      v = u;
      fv = fu;
    end
  end
end
error('too many iterations in brent');