wolffd@0: function [x, options] = linemin(f, pt, dir, fpt, options, ...
wolffd@0: 	varargin)
wolffd@0: %LINEMIN One dimensional minimization.
wolffd@0: %
wolffd@0: %	Description
wolffd@0: %	[X, OPTIONS] = LINEMIN(F, PT, DIR, FPT, OPTIONS) uses Brent's
wolffd@0: %	algorithm to find the minimum of the function F(X) along the line DIR
wolffd@0: %	through the point PT.  The function value at the starting point is
wolffd@0: %	FPT.  The point at which F has a local minimum is returned as X.  The
wolffd@0: %	function value at that point is returned in OPTIONS(8).
wolffd@0: %
wolffd@0: %	LINEMIN(F, PT, DIR, FPT, OPTIONS, P1, P2, ...) allows  additional
wolffd@0: %	arguments to be passed to F().
wolffd@0: %
wolffd@0: %	The optional parameters have the following interpretations.
wolffd@0: %
wolffd@0: %	OPTIONS(1) is set to 1 to display error values.
wolffd@0: %
wolffd@0: %	OPTIONS(2) is a measure of the absolute precision required for the
wolffd@0: %	value of X at the solution.
wolffd@0: %
wolffd@0: %	OPTIONS(3) is a measure of the precision required of the objective
wolffd@0: %	function at the solution.  Both this and the previous condition must
wolffd@0: %	be satisfied for termination.
wolffd@0: %
wolffd@0: %	OPTIONS(14) is the maximum number of iterations; default 100.
wolffd@0: %
wolffd@0: %	See also
wolffd@0: %	CONJGRAD, MINBRACK, QUASINEW
wolffd@0: %
wolffd@0: 
wolffd@0: %	Copyright (c) Ian T Nabney (1996-2001)
wolffd@0: 
wolffd@0: % Set up the options.
wolffd@0: if(options(14))
wolffd@0:   niters = options(14);
wolffd@0: else
wolffd@0:   niters = 100;
wolffd@0: end
wolffd@0: options(10) = 0; % Initialise count of function evaluations
wolffd@0: 
wolffd@0: display = options(1);
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: cphi = 1 - 1/phi;
wolffd@0: TOL = sqrt(eps);	% Maximal fractional precision
wolffd@0: TINY = 1.0e-10;         % Can't use fractional precision when minimum is at 0
wolffd@0: 
wolffd@0: % Bracket the minimum
wolffd@0: [br_min, br_mid, br_max, num_evals] = feval('minbrack', 'linef', ...
wolffd@0:   0.0, 1.0, fpt, f, pt, dir, varargin{:});
wolffd@0: options(10) = options(10) + num_evals;  % Increment number of fn. evals
wolffd@0: 					% No gradient evals in minbrack
wolffd@0: 
wolffd@0: % Use Brent's algorithm to find minimum
wolffd@0: % Initialise the points and function values
wolffd@0: w = br_mid;   	% Where second from minimum is
wolffd@0: v = br_mid;   	% Previous value of w
wolffd@0: x = v;   	% Where current minimum is
wolffd@0: e = 0.0; 	% Distance moved on step before last
wolffd@0: fx = feval('linef', x, f, pt, dir, varargin{:});
wolffd@0: options(10) = options(10) + 1;
wolffd@0: fv = fx; fw = fx;
wolffd@0: 
wolffd@0: for n = 1:niters
wolffd@0:   xm = 0.5.*(br_min+br_max);  % Middle of bracket
wolffd@0:   % Make sure that tolerance is big enough
wolffd@0:   tol1 = TOL * (max(abs(x))) + TINY;
wolffd@0:   % Decide termination on absolute precision required by options(2)
wolffd@0:   if (max(abs(x - xm)) <= options(2) & br_max-br_min < 4*options(2))
wolffd@0:     options(8) = fx;
wolffd@0:     return;
wolffd@0:   end
wolffd@0:   % Check if step before last was big enough to try a parabolic step.
wolffd@0:   % Note that this will fail on first iteration, which must be a golden
wolffd@0:   % section step.
wolffd@0:   if (max(abs(e)) > tol1)
wolffd@0:     % Construct a trial parabolic fit through x, v and w
wolffd@0:     r = (fx - fv) .* (x - w);
wolffd@0:     q = (fx - fw) .* (x - v);
wolffd@0:     p = (x - v).*q - (x - w).*r;
wolffd@0:     q = 2.0 .* (q - r);
wolffd@0:     if (q > 0.0) p = -p; end
wolffd@0:     q = abs(q);
wolffd@0:     % Test if the parabolic fit is OK
wolffd@0:     if (abs(p) >= abs(0.5*q*e) | p <= q*(br_min-x) | p >= q*(br_max-x))
wolffd@0:       % No it isn't, so take a golden section step
wolffd@0:       if (x >= xm)
wolffd@0:         e = br_min-x;
wolffd@0:       else
wolffd@0:         e = br_max-x;
wolffd@0:       end
wolffd@0:       d = cphi*e;
wolffd@0:     else
wolffd@0:       % Yes it is, so take the parabolic step
wolffd@0:       e = d;
wolffd@0:       d = p/q;
wolffd@0:       u = x+d;
wolffd@0:       if (u-br_min < 2*tol1 | br_max-u < 2*tol1)
wolffd@0:         d = sign(xm-x)*tol1;
wolffd@0:       end
wolffd@0:     end
wolffd@0:   else
wolffd@0:     % Step before last not big enough, so take a golden section step
wolffd@0:     if (x >= xm)
wolffd@0:       e = br_min - x;
wolffd@0:     else
wolffd@0:       e = br_max - x;
wolffd@0:     end
wolffd@0:     d = cphi*e;
wolffd@0:   end
wolffd@0:   % Make sure that step is big enough
wolffd@0:   if (abs(d) >= tol1)
wolffd@0:     u = x+d;
wolffd@0:   else
wolffd@0:     u = x + sign(d)*tol1;
wolffd@0:   end
wolffd@0:   % Evaluate function at u
wolffd@0:   fu = feval('linef', u, f, pt, dir, varargin{:});
wolffd@0:   options(10) = options(10) + 1;
wolffd@0:   % Reorganise bracket
wolffd@0:   if (fu <= fx)
wolffd@0:     if (u >= x)
wolffd@0:       br_min = x;
wolffd@0:     else
wolffd@0:       br_max = x;
wolffd@0:     end
wolffd@0:     v = w; w = x; x = u;
wolffd@0:     fv = fw; fw = fx; fx = fu;
wolffd@0:   else
wolffd@0:     if (u < x)
wolffd@0:       br_min = u;   
wolffd@0:     else
wolffd@0:       br_max = u;
wolffd@0:     end
wolffd@0:     if (fu <= fw | w == x)
wolffd@0:       v = w; w = u;
wolffd@0:       fv = fw; fw = fu;
wolffd@0:     elseif (fu <= fv | v == x | v == w)
wolffd@0:       v = u;
wolffd@0:       fv = fu;
wolffd@0:     end
wolffd@0:   end
wolffd@0:   if (display == 1)
wolffd@0:     fprintf(1, 'Cycle %4d  Error %11.6f\n', n, fx);
wolffd@0:   end
wolffd@0: end
wolffd@0: options(8) = fx;