function y = f_bar(x1,x2)

b = [1 1 1 1 1 1]';
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];

if sum(b-A*[x1;x2]<0) == 0
    y = -sum(log(b-A*[x1;x2]));
else
    y =NaN;
end

end

%%%% script

clear all; close all;

%consts
V =[ 0.1562 0.9127 1.0338 0.8086 -1.3895 -0.8782
    -1.0580 -0.6358 0.1386 0.6406 2.3203 -0.8311];
c = [ 1 -1]';
b = [1 1 1 1 1 1]';
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];

%params
x0 = [0 0]';
%x0 = [-0.6 1.6]';
Tinit = 0.2;
m=6;
bb=0.8;
aa=0.5;
gamma=2;
eps = 10^-10;

%macros
psi=@(x,t) t*c'*x - sum(log(b-A*x));
grad=@(x,t) t*c +(A')*(1./(b-A*x));
hessian=@(x) ((diag(1./(b-A*x))*A)')*(diag(1./(b-A*x))*A);


t=Tinit;
x=x0;
X(:,1)= x0;

k=1;
while(m/t>eps)
    s=1;
    while psi( x-s*(hessian(x)^-1)*grad(x,t),t ) >= psi(x, t)-s*aa*grad(x, t)'*(hessian(x)^-1)*grad(x,t)
        s = bb*s;
    end
    x = x-s*(hessian(x)^-1)*grad(x,t);
    t=t*gamma;

    k=k+1;
    X(:,k)=x;
end

Xlin= linprog(c,A,b);

%Display results
disp("*********************************************")
disp("Rozwiazanie:")
disp("Parametry:  gamma: " + gamma + ",  beta: " + bb + ",  alfa: " + aa + ",  eps: " + eps + ", t_init: " + Tinit + ",  t_koncowe: " + t)
disp("Punkt startowy: [" + x0(1)+","+x0(2) + "]")
disp("Liczba iteracji: " + k)
disp("Punkt koncowy: [" + X(1,k)+","+X(2,k) + "]")
disp("Rozw linprog: [" + Xlin(1,1)+","+Xlin(2,1) + "]")

%Plots
x1 = -2:0.01:1.5; x2 = -1.5:0.01:3;
[X1,X2]  = meshgrid(x1,x2);
Z=arrayfun(@f_bar, X1, X2);
values = arrayfun(@f_bar, X(1,:), X(2,:));
curves = values(1,2:2:10);

figure
fill(V(1,:), V(2,:), 'g'); hold on;
%contour(X1, X2, Z, curves, 'ShowText', 'on')
contour(X1, X2, Z, [-0.6 0.12 1.24], 'ShowText', 'on')

for k=1:size(X,2)
   hold on
   plot(X(1,k), X(2,k), 'or')
end
title('Wykres obszaru dopuszczalnego i punktów centralnych.')
xlabel('X1')
ylabel('X2')