function xopt = gold_sect_search (fhandle, bracket, accuracy, debug_flag)

%GOLD_SECTION_SEARCH  optimization using golden section search
%   XOPT = GOLD_SECT_SEARCH (FHANDLE, BRACKET, ACCURACY, DEBUG_FLAG) uses
%   golden section search to find a minimum of a one-dimensional function
%   FHANDLE within the interval BRACKET with an accuracy of at least
%   ACCURACY.  The X value at which the function reaches a minimum is
%   returned as XOPT.
%
%   Inputs:
%      FHANDLE:    The function handle to the function to be minimized
%      BRACKET:    Bracket(1) is the lower bound, Bracket(2) the upper
%                  bound of the interval within which the minimum is
%                  searched.
%      ACCURACY:   Bound on the absolute error on XOPT.
%      DEBUG_FLAG: (optional) when non-zero, the function generates one
%                  line of output for every iteration.
%
%   Outputs:
%      XOPT:       The X value at which the function reaches a minimum
%
%   Assumptions:
%            - the function FHANDLE is unimodal
%            - the minimum is reached in at most 2000 iterations
%
%% ME 2016
%% Fall 2007
%%
%% AUTHOR:  Chris Paredis
%%

% NOTE: this funciton implements MINIMIZATION!!!
maxiter = 2000;
ratio = 0.5*(sqrt(5)-1);

% initialize the bracket and two interior points
xl = bracket(1);
xu = bracket(2);
size = ratio*(xu - xl);
x1 = xl + size;
x2 = xu - size;  % Note: x2<x1
% the corresponding function values
fl = feval(fhandle,xl);
fu = feval(fhandle,xu);
f1 = feval(fhandle,x1);
f2 = feval(fhandle,x2);

for i=5:maxiter
    % check whether we have achieved the desired accuracy
    if 0.5*(xu - xl) <= accuracy
        % the best guess would be the middle of the current bracket
        xopt = 0.5*(xu + xl);
        return;
    end
    
    % update the size of the next bracket
    size = ratio * size;
    
    % find the new bracket
    if f1 > f2
        % [xl, x1] is the new bracket
        xu = x1;
        x1 = x2;
        x2 = xu - size;
        fu = f1;
        f1 = f2;
        f2 = feval(fhandle, x2);
        if debug_flag
            fprintf('%5d | %12g | %12g\n',i,x2,f2);
        end
    else
        % [x2, xu] is the new bracket
        xl = x2;
        x2 = x1;
        x1 = xl + size;
        fl = f2;
        f2 = f1;
        f1 = feval(fhandle, x1);
        if debug_flag
            fprintf('%5d | %12g | %12g\n',i,x1,f1);
        end
    end
end
warning('maximum number of iterations has been exceeded');