%  
% Perform Golden Section search on the function fun,
% over the interval [a,b].  The function fun must be UNIMODAL
% over [a,b].  Search stops when error <= tol(b-a).
%
% This is an example of how you would call this function
% [xmin,fxmin,newa,newb]=golden('f',a,b,tol)
%
% golden.m returns the estimated minimum point,
%                  the estimated objective,
%                  the final estimated interval.
%
%
function [xmin,fxmin,newa,newb]=golden(fun,a,b,tol)
%
tic
TAU = (-1.0 + sqrt(5.0))/2.0;

%     the current search interval will be [ak, bk],
%     with function evaluations at points left, right where
%            temp = TAU*(bk - ak)
%            right = ak + TAU*(bk - ak)
%            left = bk - temp 
%      or equivalently 
%            left = ak + bk - right;    
%     i.e. ak < left < right < bk.


%
disp('Start Golden Section Search:');
disp('----------------------------');

if ((b-a < 0)| (tol<=0))
   error('Bad input data: TERMINATE');
else

  disp('The search is terminated when the final interval has');
  disp('length no greater than ');
    acc = tol*(b-a) 

%   initialize ... 
    ak = a; bk = b;
    fl = feval(fun,ak) ;
    fr = feval(fun,bk) ;
    numEval = 2;
    right = a + TAU*(b-a);
    fright = feval(fun,right); 

    left = ak + bk - right;
    fleft = feval(fun,left); 

    disp('   ak            bk           f(Left)       f(Right)');
  
    while ((bk-ak)>acc)
  
      s=sprintf('%7.5f        %7.5f       %7.5f       %7.5f', ak,bk,fl,fr);
      disp(s)
      numEval=numEval+1;
%      decide which subinterval to discard
      if (fleft < fright)
        bk = right; fr = fright;
        right = left; fright=fleft;
        left = ak + bk - right;
        fleft= feval(fun,left);

      else if (fleft >= fright)
        ak = left; fl =fleft;
        left = right; fleft=fright;
        right = ak + bk - left;
        fright = feval(fun,right);
      end;
    end;
end;

s=sprintf('\nThe final interval has endpoints %8.6f, %8.6f.', ak, bk);
disp(s);
s=sprintf('Predicted no. of function evaluations = %d.', 2 + ceil( log(tol)/log(TAU)) );
disp(s);
s=sprintf('Required  no. of function evaluations = %d.', numEval);
disp(s);

xmin=(ak+bk)/2;
fxmin=feval(fun,xmin);
newa=ak;
newb=bk;
%
% End of the function.  
%

toc
end;