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