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