clear;
close all;

A = [0.4873 -0.8732;
    0.6072 0.7946;
    0.9880 -0.1546;
    -0.2142 -0.9768;
    -0.9871 -0.1601;
    0.9124 0.4093];

b = ones(1,6)';
c = [1 ;-1];
x0 = [0 0]';

[m,n] = size(A);

t_init = 0.2;
t = t_init;
gamma = 2;

eps = 10^-15;

alfa = 0.5;
beta = 0.8;

xk = x0;
sk = 0;

psi = @(x,t) t*c'*x - sum(log(b-A*x));
grad_psi = @(x,t) t*c +(A')*(1./(b-A*x));
hessian_psi = @(x)((diag(1./(b-A*x))*A)')*(diag(1./(b-A*x))*A);

results = [];

while(m/t>eps)

    s = 1;
    while(psi(xk-s*(hessian_psi(xk)^-1)*grad_psi(xk,t),t) >= psi(xk,t)-s*alfa*grad_psi(xk,t)'*(hessian_psi(xk)^-1)*grad_psi(xk,t))
        s = s*beta;
    end

    sk = s;
    xk = xk-sk*(hessian_psi(xk)^-1)*grad_psi(xk,t);
    results = [results,xk];
    t = gamma * t;

end


xk

show_figure(A,b,results,x0);

function y = fbar (x1, x2)
A = [0.4873 -0.8732;
    0.6072 0.7946;
    0.9880 -0.1546;
    -0.2142 -0.9768;
    -0.9871 -0.1601;
    0.9124 0.4093];

b = ones(1,6)';
x = [x1 ; x2];
    if sum(b - A * x < 0) == 0
        y = -sum(log(b - A * x));
    else
        y = NaN;
    end
end

function show_figure(A,b,results,x0)

barrier_f = @(x) -sum(log(b-A*x));

P = polytope(A,b);
plot(P);
xlim([-2.5 2])
ylim([-1.5 3])
xlabel('$x_{1}$','Interpreter','latex')
ylabel('$x_{2}$','Interpreter','latex')
hold on;
scatter(results(1,:),results(2,:),20,'MarkerEdgeColor','b','MarkerFaceColor','b');
scatter(x0(1),x0(2),20,'MarkerEdgeColor','black','MarkerFaceColor','black');

x1_lim = xlim;
x2_lim = ylim;

x1_lim = [x1_lim(1)-0.5,x1_lim(2)+0.5];
x2_lim = [x2_lim(1)-0.5,x2_lim(2)+0.5];

x1 = x1_lim(1):0.01:x1_lim(2);
x2 = x2_lim(1):0.01:x2_lim(2);

[X1, X2] = meshgrid(x1, x2);
Z = arrayfun(@fbar, X1, X2);

v = real(barrier_f(results));
[C,h] = contour(X1, X2, Z, v,'color','black');
clabel(C,h)

end