wolffd@0: function  [br_min, br_mid, br_max, num_evals] = minbrack(f, a, b, fa,  ...
wolffd@0: 			 varargin)
wolffd@0: %MINBRACK Bracket a minimum of a function of one variable.
wolffd@0: %
wolffd@0: %	Description
wolffd@0: %	BRMIN, BRMID, BRMAX, NUMEVALS] = MINBRACK(F, A, B, FA) finds a
wolffd@0: %	bracket of three points around a local minimum of F.  The function F
wolffd@0: %	must have a one dimensional domain. A < B is an initial guess at the
wolffd@0: %	minimum and maximum points of a bracket, but MINBRACK will search
wolffd@0: %	outside this interval if necessary. The bracket consists of three
wolffd@0: %	points (in increasing order) such that F(BRMID) < F(BRMIN) and
wolffd@0: %	F(BRMID) < F(BRMAX). FA is the value of the function at A: it is
wolffd@0: %	included to avoid unnecessary function evaluations in the
wolffd@0: %	optimization routines. The return value NUMEVALS is the number of
wolffd@0: %	function evaluations in MINBRACK.
wolffd@0: %
wolffd@0: %	MINBRACK(F, A, B, FA, P1, P2, ...) allows additional arguments to be
wolffd@0: %	passed to F
wolffd@0: %
wolffd@0: %	See also
wolffd@0: %	LINEMIN, LINEF
wolffd@0: %
wolffd@0: 
wolffd@0: %	Copyright (c) Ian T Nabney (1996-2001)
wolffd@0: 
wolffd@0: % Check function string
wolffd@0: f = fcnchk(f, length(varargin));
wolffd@0: 
wolffd@0: % Value of golden section (1 + sqrt(5))/2.0
wolffd@0: phi = 1.6180339887499;
wolffd@0: 
wolffd@0: % Initialise count of number of function evaluations
wolffd@0: num_evals = 0;
wolffd@0: 
wolffd@0: % A small non-zero number to avoid dividing by zero in quadratic interpolation
wolffd@0: TINY = 1.e-10;
wolffd@0: 
wolffd@0: % Maximal proportional step to take: don't want to make this too big
wolffd@0: % as then spend a lot of time finding the minimum inside the bracket
wolffd@0: max_step = 10.0;
wolffd@0: 
wolffd@0: fb = feval(f, b, varargin{:});
wolffd@0: num_evals = num_evals + 1;
wolffd@0: 
wolffd@0: % Assume that we know going from a to b is downhill initially 
wolffd@0: % (usually because gradf(a) < 0).
wolffd@0: if (fb > fa)
wolffd@0:   % Minimum must lie between a and b: do golden section until we find point
wolffd@0:   % low enough to be middle of bracket
wolffd@0:   c = b;
wolffd@0:   b = a + (c-a)/phi;
wolffd@0:   fb = feval(f, b, varargin{:});
wolffd@0:   num_evals = num_evals + 1;
wolffd@0:   while (fb > fa)
wolffd@0:     c = b;
wolffd@0:     b = a + (c-a)/phi;
wolffd@0:     fb = feval(f, b, varargin{:});
wolffd@0:     num_evals = num_evals + 1;
wolffd@0:   end
wolffd@0: else  
wolffd@0:   % There is a valid bracket upper bound greater than b
wolffd@0:   c = b + phi*(b-a);
wolffd@0:   fc = feval(f, c, varargin{:});
wolffd@0:   num_evals = num_evals + 1;
wolffd@0:   bracket_found = 0;
wolffd@0:   
wolffd@0:   while (fb > fc)
wolffd@0:     % Do a quadratic interpolation (i.e. to minimum of quadratic)
wolffd@0:     r = (b-a).*(fb-fc);
wolffd@0:     q = (b-c).*(fb-fa);
wolffd@0:     u = b - ((b-c)*q - (b-a)*r)/(2.0*(sign(q-r)*max([abs(q-r), TINY])));
wolffd@0:     ulimit = b + max_step*(c-b);
wolffd@0:     
wolffd@0:     if ((b-u)'*(u-c) > 0.0)
wolffd@0:       % Interpolant lies between b and c
wolffd@0:       fu = feval(f, u, varargin{:});
wolffd@0:       num_evals = num_evals + 1;
wolffd@0:       if (fu < fc)
wolffd@0: 	% Have a minimum between b and c
wolffd@0: 	br_min = b;
wolffd@0: 	br_mid = u;
wolffd@0: 	br_max = c;
wolffd@0: 	return;
wolffd@0:       elseif (fu > fb)
wolffd@0: 	% Have a minimum between a and u
wolffd@0: 	br_min = a;
wolffd@0: 	br_mid = c;
wolffd@0: 	br_max = u;
wolffd@0: 	return;
wolffd@0:       end
wolffd@0:       % Quadratic interpolation didn't give a bracket, so take a golden step
wolffd@0:       u = c + phi*(c-b);
wolffd@0:     elseif ((c-u)'*(u-ulimit) > 0.0)
wolffd@0:       % Interpolant lies between c and limit
wolffd@0:       fu = feval(f, u, varargin{:});
wolffd@0:       num_evals = num_evals + 1;
wolffd@0:       if (fu < fc)
wolffd@0: 	% Move bracket along, and then take a golden section step
wolffd@0: 	b = c;
wolffd@0: 	c = u;
wolffd@0: 	u = c + phi*(c-b);
wolffd@0:       else
wolffd@0: 	bracket_found = 1;
wolffd@0:       end
wolffd@0:     elseif ((u-ulimit)'*(ulimit-c) >= 0.0)
wolffd@0:       % Limit parabolic u to maximum value
wolffd@0:       u = ulimit;
wolffd@0:     else
wolffd@0:       % Reject parabolic u and use golden section step
wolffd@0:       u = c + phi*(c-b);
wolffd@0:     end
wolffd@0:     if ~bracket_found
wolffd@0:       fu = feval(f, u, varargin{:});
wolffd@0:       num_evals = num_evals + 1;
wolffd@0:     end
wolffd@0:     a = b; b = c; c = u;
wolffd@0:     fa = fb; fb = fc; fc = fu;
wolffd@0:   end % while loop
wolffd@0: end   % bracket found
wolffd@0: br_mid = b;
wolffd@0: if (a < c)
wolffd@0:   br_min = a;
wolffd@0:   br_max = c;
wolffd@0: else
wolffd@0:   br_min = c; 
wolffd@0:   br_max = a;
wolffd@0: end