function [xf,fval,exitflag,output] = fminbnd(funfcn,ax,bx,options,varargin)
%FMINBND Scalar bounded nonlinear function minimization.  
%   X = FMINBND(FUN,x1,x2) returns a value X that is a local minimizer
%   of the function that is described in FUN (usually an M-file)
%   in the interval x1 < X < x2. The function FUN should return a scalar 
%   function value F when called with feval: F=feval(FUN,X). 
%   See below for more options for FUN.
%
%   X = FMINBND(FUN,x1,x2,OPTIONS) minimizes with the default optimization
%   parameters replaced by values in the structure OPTIONS, created with
%   the OPTIMSET function.  See OPTIMSET for details.  FMINBND uses these options: 
%   Display, TolX, MaxFunEval and MaxIter. 
%
%   X = FMINBND(FUN,x1,x2,OPTIONS,P1,P2,...) provides for additional
%   arguments, which are passed to the objective function, FUN(X,P1,P2,...).
%   (Use OPTIONS = [] as a place holder if no options are set.)
%
%   [X,FVAL] = FMINBND(...) also returns the value of the objective function,
%   FVAL, computed in FUN, at X.
%
%   [X,FVAL,EXITFLAG] = FMINBND(...) also returns a string EXITFLAG that 
%   describes the exit condition of FMINBND.  
%   If EXITFLAG is:
%     1 then FMINBND converged with a solution X based on OPTIONS.TolFun
%     0 then the maximum number of function evaluations was reached
%   
%   [X,FVAL,EXITFLAG,OUTPUT] = FMINBND(...) also returns a structure
%   OUTPUT with the number of iterations taken in OUTPUT.iterations.
%
%   Examples
%       x = fminbnd('cos',3,4) computes pi to a few decimal places 
%         and gives a message upon termination.
%       [x,fval,exitflag] = fminbnd('cos',3,4,optimset('TolX',1e-12,'Display','off')) 
%         computes pi to about 12 decimal places, suppresses output, returns the  
%         function value at x, and returns an exitflag of 1.
%   
%   The argument FUN can be an inline function:
%      f = inline('sin(x)+3');
%      x = fminbnd(f,2,5);
%   Or create an M-file function afun.m on your path that contains
%     function f = afun(x)
%       f = sin(x) = 3;
%   then call fminbnd
%     x = fminbnd('afun',2,5);
%     
%
%   See also OPTIMSET, FMINSEARCH, FZERO.

%   Reference: "Computer Methods for Mathematical Computations",
%   Forsythe, Malcolm, and Moler, Prentice-Hall, 1976.

%   Original coding by Duane Hanselman, University of Maine.
%   Copyright (c) 1984-98 by The MathWorks, Inc.
%   $Revision: 1.6 $  $Date: 1998/08/17 19:31:57 $


defaultopt = optimset('display','final',...
      'maxfuneval',500,'maxiter',500,'TolX',1e-4);
% If just 'defaults' passed in, return the default options in X
if nargin==1 & nargout <= 1 & isequal(funfcn,'defaults')
   xf = defaultopt;
   return
end

if nargin < 3
   error('fminbnd requires three input arguments.')
end

% initialization
if nargin<4, 
   options = []; 
end

options = optimset(defaultopt,options);
printtype = optimget(options,'display');
tol = optimget(options,'tolx');
maxfun = optimget(options,'maxfuneval');
maxiter = optimget(options,'maxiter');
% fun evals are the same as iterations, so just check the minimum
maxcount = min(maxfun,maxiter);

switch printtype
case {'none','off'}
   print = 0;
case 'iter'
   print = 2;
case 'final'
   print = 1;
otherwise
   print = 1;
end
% checkbounds
if ax > bx
   exitflag = -1;
   output = []; xf=[]; fval = [];
   if print > 0
        msg=sprintf(['\nExiting due to infeasibility: the lower bound exceeds the' ...
         ' upper bound.\n']);
        disp(msg)
   end
   return
end

% Assume we'll converge
exitflag = 1;

header = ' Func-count     x          f(x)         Procedure';
step='       initial';

% Convert to inline function as needed.
funfcn = fcnchk(funfcn,length(varargin));

seps = sqrt(eps);
c = 0.5*(3.0 - sqrt(5.0));
a = ax; b = bx;
v = a + c*(b-a);
w = v; xf = v;
d = 0.0; e = 0.0;
x= xf; fx = feval(funfcn,x,varargin{:}); 
num = 1;
fmin_data = [ 1 xf fx ];
if print > 1
   disp(' ')
   disp(header)
   disp([sprintf('%5.0f   %12.6g %12.6g ',fmin_data), step]) 
end
fv = fx; fw = fx;
xm = 0.5*(a+b);
tol1 = seps*abs(xf) + tol/3.0;   
tol2 = 2.0*tol1;

% Main loop
while ( abs(xf-xm) > (tol2 - 0.5*(b-a)) ) 
   gs = 1;
   % Is a parabolic fit possible
   if abs(e) > tol1
      % Yes, so fit parabola
      gs = 0;
      r = (xf-w)*(fx-fv);
      q = (xf-v)*(fx-fw);
      p = (xf-v)*q-(xf-w)*r;
      q = 2.0*(q-r);
      if q > 0.0,  p = -p; end
      q = abs(q);
      r = e;  e = d;
      
      % Is the parabola acceptable
      if ( (abs(p)<abs(0.5*q*r)) & (p>q*(a-xf)) & (p<q*(b-xf)) )
         
         % Yes, parabolic interpolation step
         d = p/q;
         x = xf+d;
         step = '       parabolic';
         
         % f must not be evaluated too close to ax or bx
         if ((x-a) < tol2) | ((b-x) < tol2)
            si = sign(xm-xf) + ((xm-xf) == 0);
            d = tol1*si;
         end
      else
         % Not acceptable, must do a golden section step
         gs=1;
      end
   end
   if gs
      % A golden-section step is required
      if xf >= xm, e = a-xf;    else, e = b-xf;  end
      d = c*e;
      step = '       golden';
   end
   
   % The function must not be evaluated too close to xf
   si = sign(d) + (d == 0);
   x = xf + si * max( abs(d), tol1 );
   fu = feval(funfcn,x,varargin{:});  
   num = num+1;
   fmin_data = [num x fu];
   if print > 1
      disp([sprintf('%5.0f   %12.6g %12.6g ',fmin_data), step]) 
   end
   
   % Update a, b, v, w, x, xm, tol1, tol2
   if fu <= fx
      if x >= xf, a = xf; else, b = xf; end
      v = w; fv = fw;
      w = xf; fw = fx;
      xf = x; fx = fu;
   else % fu > fx
      if x < xf, a = x; else,b = x; end
      if ( (fu <= fw) | (w == xf) )
         v = w; fv = fw;
         w = x; fw = fu;
      elseif ( (fu <= fv) | (v == xf) | (v == w) )
         v = x; fv = fu;
      end
   end
   xm = 0.5*(a+b);
   tol1 = seps*abs(xf) + tol/3.0; tol2 = 2.0*tol1;
   if num > maxcount
      if print > 0     
         exitflag = 0;
         output.iterations=num;
         output.funcCount=num;
         fval = fx;
         terminate(x,exitflag,fval,maxfun,maxiter,options);
         return
      end
   end
end
fval = fx;
output.iterations = num;
output.funcCount = num;
output.algorithm = 'golden section search, parabolic interpolation';
if print > 0
   terminate(x,exitflag,fval,maxfun,maxiter,options);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function terminate(x,exitflag,finalf,maxfun,maxiter,options)

switch exitflag
case 1
   convmsg1 = sprintf([ ...
         '\nOptimization terminated successfully:\n',...
         ' the current x satisfies the termination criteria using OPTIONS.TolX of %e \n',...
          ],options.TolX);
   disp(convmsg1)
   
case 0
   if maxfun <= maxiter
      disp(' ')
      disp('Exiting: Maximum number of function evaluations has been exceeded')
      disp('         - increase MaxFunEvals option.')
      msg = sprintf('         Current function value: %f \n', finalf);
   else
      disp(' ')
      disp('Exiting: Maximum number of iterations has been exceeded')
      disp('         - increase MaxIter option.')
      msg = sprintf('         Current function value: %f \n', finalf);
      disp(msg)
   end
   
case -1
   ;
otherwise
   ;
end
