function [x_min,f_min,b_iter,s_iter] = line_search(f,x1,lb,ub,h,tol,max_b_iter,max_s_iter)
% This function finds the minimum of f(x) within lower and upper bounds.
% f = f(x).
% x1 = starting value of x.
% h = initial step size.
% lb = lower bound.
% ub = upper bound.
% tol = error tolerance on the minimum of f(x).
% max_b_iter = maximum number of iterations in the bracketing phase.
% max_s_iter = maximum number of iterations in the sectioning phase.

if nargin < 8; max_s_iter = 1000; end
if nargin < 7; max_b_iter = 1000; end
if nargin < 6; tol = 1.0e-8; end
if nargin < 5; h = 0.1; end
if nargin < 4; ub = Inf; end
if nargin < 3; lb = -Inf; end

% Golden ratio:
p = [-1 -1 1];
r = roots(p);
c = -r(1);
R = r(2);
C = 1.0 - R;

% Bracketing phase.
% Find the direction downhill from x1.
f1 = f(x1);
x2 = x1 + h; 
f2 = f(x2);

if f2 > f1
    h = -h;
    x2 = x1 + h; 
    f2 = f(x2);
end

% Iterate.
b_iter = 0;
while f2 <= f1 && b_iter<=max_b_iter
    b_iter = b_iter + 1;
    h = c*h;
    x1 = x2; f1 = f2;
    x2 = x2 + h; 
    f2 = f(x2);
end
a = x1; b = x2; 

if b_iter >= max_b_iter       
    disp('Error: maximum number of bracketing iterations exceeded')
end
if a <= lb || a >= ub || b <= lb || b >= ub      
    disp('Error: minimum is not within bounds')
end

% Sectioning phase.
% Initialize x1 and x2:
x1 = R*a + C*b;
x2 = C*a + R*b;

% The following shows the equivalence between the previous two lines and
% the graphical derivation of x1 and x2.
% x1 = Ra + Cb
% x1 = Ra + (1 - R)b
% x1 = Ra + b - Rb
% x1 = b - R(b - a) = b - Rh 
% x2 = (1 - R)a + Rb
% x2 = a - Ra + Rb
% x2 = a + R(b - a) = a + Rh

f1 = f(x1);
f2 = f(x2);

% Iterate.
s_iter = 0;
while abs(b-a) > tol*(abs(x1) + abs(x2))
    s_iter = s_iter + 1;
    if f1 > f2
        a = x1; x1 = x2; f1 = f2;
        x2 = C*a + R*b;
        f2 = f(x2);
    else
        b = x2; x2 = x1; f2 = f1;
        x1 = R*a + C*b;
        f1 = f(x1);
    end
end

if s_iter >= max_s_iter       
    disp('Error: maximum number of sectioning iterations exceeded')
end

if f1 < f2 
    f_min = f1; x_min = x1;
else
    f_min = f2; x_min = x2;
end

end

