function [x,y] = localmin(func, a, b, opts, varargin)
%
% LOCALMIN - constrained minimisation of a univariate function
%
%    X = LOCALMIN(FUNC,A,B) finds a local minumum of a univariate function, 
%    via a function handle FUNC, constrained to lie in the range (A, B).  The
%    local minimum is located via an algorithm analogous to Dekker's method
%    [1], which does not require any gradient information.
%
%    [X,Y] = LOCALMIN(FUNC,A,B) also returns the value of the function, Y, at
%    the local minimum, X.
%
%    X = LOCALMIN(FUNC,A,B,OPTS) allows the optimisation parameters to be
%    specified via a structure, OPTS, with members
%
%       opts.MaxIter     - Maximum number of optimisation steps.
%       opts.MaxFunEvals - Maximum number of function evaluations.
%       opts.TolFun      - Stopping criterion based on the relative change in
%                          value of the function in each step.
%       opts.TolX        - Stopping criterion based on the change in the
%                          minimiser in each step.
%
%    OPTS = LOCALMIN() returns a structure containing the default optimisation
%    parameters, with the following values:
%
%       opts.MaxIter     = 200
%       opts.MaxFunEvals = 1000
%       opts.TolFun      = 1e-3
%       opts.TolX        = 1e-3
%
%    X = LOCALMIN(FUNC,A,B,OPTS, P1, P2, ...) allows addinal parameters to be
%    passed to the function to be minimised.
%
%    References:
%
%       [1] Richard P. Brent, "Algorithms for minimization without
%           derivatives", Dover, 2002.  ISBN 0-486-41998-3
%
%    See also: SIMPLEX

%
% File        : localmin.m
%
% Date        : Friday 3rd July 2009
%
% Author      : Dr Gavin C. Cawley
%
% Description : Simple implementation of an algorithm for bounded minimisation
%               of a scalar function without gradient information 
%               [1, chapter 5],  Similar to the fminbnd routine from the
%               MATLAB optimisation toolbox.
%
% References  : [1] Richard P. Brent, "Algorithms for minimization without
%                   derivatives", Dover, 2002.  ISBN 0-486-41998-3
%
% History     : 03/07/2009 - v1.00
%
% Copyright   : (c) Dr Gavin C. Cawley, July 2009.
%
%    This program is free software; you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation; either version 2 of the License, or
%    (at your option) any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program; if not, write to the Free Software
%    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
%

% use default optimisation parameters if none given

if nargin < 4

   opts.MaxIter     = 20;
   opts.MaxFunEvals = 1000;
   opts.TolFun      = 1e-9;
   opts.TolX        = 1e-9;
   opts.Display     = 'iter';

end

% return structure containing default optimisation parameters

if nargin == 0

   x = opts;

   return

end

% initialisation

c      = 0.38196601125011;   % c = (3 - sqrt(5))/2;
x      = a + c*(b - a);
w      = x;
v      = x;
e      = 0;
fx     = func(x, varargin{:}); 
fw     = fx;
fv     = fx;
fcount = 1;
format = '    %d\t\t%f\t%f\t%s\n';

fprintf(1, '\n Func-count        x               f(x)        Procedure\n');
fprintf(1, format, fcount, x, fx, 'initial');

% get on with it!

for i=1:opts.MaxIter

   m    = 0.5*(a + b);
   tol  = eps*abs(x) + opts.TolFun;
   t2   = 2*tol;

   % check stopping criterion

   if abs(x - m) < (t2 - 0.5*(b-a))

      break;

   end

   if abs(e) > tol

      % fit parabola

      r = (x - w)*(fx- fv);
      q = (x - v)*(fx - fw);
      p = (x - v)*q - (x - w)*r;
      q = 2*(q - r);

      if q > 0

         p = -p;

      else

         q = -q;

      end

      r = e;
      e = d;

   else

      p = 0;
      q = 0;
      r = 0;

   end

   if abs(p) < abs(0.5*q*r) && p > q*(a - x) && p < q*(b - x) 

      % perform parabolic interpolation step

      method = 'parabolic';
      d      = p/q;
      u      = x + d;
      
      % f must not be evaluated too close to a or b

      if u - a < t2 || b - u < t2

         if x < m

            d = tol; 

         else

            d = -tol;

         end

      end

   else

      % perform  golden section step

      method = 'golden';

      if x < m

         e = b - x;

      else

         e = a - x;

      end

      d = c*e;

   end

   % f must not be evaluated too close to x

   if abs(d) >= tol

      u = x + d;

   elseif d > 0

      u = x + tol;

   else

      u = x - tol;

   end

   fu     = func(u, varargin{:});
   fcount = fcount + 1;

   % update a, b, v, w and x

   if fu < fx

      if u < x

         b = x;
         
      else

         a = x;

      end

      v  = w;
      fv = fw;
      w  = x;
      fw = fx;
      x  = u; 
      fx = fu;

   else

      if u < x

         a = u;

      else

         b = u;

      end

      if fu <= fw || w == x

         v  = w;
         fv = fw;
         w  = u;
         fw = fu;

      elseif fu <= fv || v == x ||v == w

         v  = u;
         fv = fu;

      end

   end

   % display progress

   fprintf(1, format, fcount, x, fx, method);

end

fprintf(1, '\n');

% bye bye...

