function [a, b, N, n] = min_fibonacci(fun_expression, alpha, beta, epsilon, lambda)
% Fibonacci method for finding the local minimum of a function.
%
% fun_expr: (symbolic expression over x) The symbolic expression of the
%           objective function
% alpha:    (number) The starting point of the interval in which we seek
%           for minimum
% beta:     (number) The ending point of the interval in which we seek
%           for minimum
% epsilon:  (number) The epsilon value (the interval of the last step)
% lambda:   (number) The lambda value (accuracy)
%
% return:
%  a:       (vector) Starting points of the interval for each iteration
%  b:       (vector) Ending points of the interval for each iteration
%  N:       (number) The number of iterations needed.
%  nn:      (number) The calls of objective function fun_expr
%

% Error checking
if alpha > beta || lambda <= 0 || epsilon <= 0
    error ('Input criteria not met')
end

% Use Binet's formula instead of matlab's recursive fibonacci
% implementation
fibonacci = @(n) ( ((1 + sqrt(5))^n - (1 - sqrt(5))^n) / (2^n * sqrt(5)) );

% Init variables
a = alpha;
b = beta;
n = 0;
fun = matlabFunction(fun_expression);

% wrapper call count function
function r = count_fun(x)
    n = n + 1;
    r = fun(x);
end

% calculate number of iterations
N = 0;
while fibonacci(N) < (b(1) - a(1)) / lambda
    N = N + 1;
end


% calculate x1, x2 of the first iteration, since the following iteration
% will not require to calculate both
x_1 = a(1) + (fibonacci(N-2) / fibonacci(N)) * (b(1) - a(1));
x_2 = a(1) + (fibonacci(N-1) / fibonacci(N)) * (b(1) - a(1));
f1 = count_fun(x_1);
f2 = count_fun(x_2);

% All but the last calculation
for k = 1:N-2
    % set new search interval
    if f1 <= f2
        a(k+1) = a(k);
        b(k+1) = x_2;
        x_2 = x_1;
        f2 = f1;
        x_1 = a(k+1) + (fibonacci(N-k-2) / fibonacci(N-k)) * (b(k+1) - a(k+1));
        f1 = count_fun(x_1);
    else
        a(k+1) = x_1;
        b(k+1) = b(k);
        x_1 = x_2;
        f1 = f2;
        x_2 = a(k+1) + (fibonacci(N-k-1) / fibonacci(N-k)) * (b(k+1) - a(k+1));
        f2 = count_fun(x_2);
    end
end

% Last calculation
x_2 = x_1 + epsilon;
f2 = count_fun(x_2);
if f1 <= f2
    a(N) = a(N-1);
    b(N) = x_1;
else
    a(N) = x_1;
    b(N) = b(N-1);
end

end