comparison toolboxes/FullBNT-1.0.7/netlab3.3/minbrack.m @ 0:e9a9cd732c1e tip

first hg version after svn

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