close all;
clear all;

mu_1 = [0 0]';
mu_2 = [5 0]';
% Sigma_1 = eye(2,2);
Sigma_1 = [1 0; 0 1];
Sigma_2 = [4 0; 0 4];
rad_1 = sqrt(Sigma_1(1,1));
rad_2 = sqrt(Sigma_2(1,1));

x_points_1 = [0:100]*2/100 - 1;
x_points_1 = rad_1*x_points_1 + mu_1(1);
y_points_1 = sqrt(rad_1^2 - (x_points_1 - mu_1(1)).^2);


x_points_2 = [0:100]*2/100 - 1;
x_points_2 = rad_2*x_points_2 + mu_2(1);
y_points_2 = sqrt(rad_2^2 - (x_points_2 - mu_2(1)).^2);
%%
figure;
plot(x_points_1,y_points_1+ mu_1(2),'r-');
hold on;
plot(x_points_1,-y_points_1+ mu_1(2),'r-','LineWidth',2);
plot(x_points_2,y_points_2+ mu_2(2),'g-','LineWidth',2);
plot(x_points_2,-y_points_2+ mu_2(2),'g-','LineWidth',2);
axis equal;
xlabel('x_1'); ylabel('x_2');

%%
%%

g = inline('X'' *( -1/2 * iE) * X + (iE * m)'' * X - 1/2 * m'' * iE * m - 1/2 * log(det(E)) + log(prior) ','X','m','E','iE','prior');
prior_1 = 20/3%1/4;%10/3;
prior_2 = 1%1 - prior_1;
x_axis = get(gca,'xLim'); y_axis = get(gca,'yLim');
x_vect = linspace(x_axis(1),x_axis(2),60);
y_vect = linspace(y_axis(1),y_axis(2),60);
Z_1 = zeros(length(x_vect),length(y_vect));
Z_2 = Z_1;
iSigma_1 = inv(Sigma_1);
iSigma_2 = inv(Sigma_2);
for i = 1:length(x_vect)
    for j = 1:length(y_vect)
        [i j ]
        Z_1(i,j) = g([x_vect(i) y_vect(j)]',mu_1,Sigma_1,iSigma_1,prior_1);
        Z_2(i,j) = g([x_vect(i) y_vect(j)]',mu_2,Sigma_2,iSigma_2,prior_2);
    end
end

%%
[X Y] = meshgrid(x_vect, y_vect);

figure;

imagesc(X(:),Y(:),abs(Z_2-Z_1)')
set(gca,'YDir','normal')
hold on;

[C,h] = contour(X,Y,(Z_2-Z_1)',[0 0]);
% set(h,'ShowText','on','TextStep',get(h,'LevelStep')*2,'Edgecolor',[1 1 1]);
colormap cool

plot(x_points_1,y_points_1+ mu_1(2),'r-','LineWidth',2);
plot(x_points_1,-y_points_1+ mu_1(2),'r-','LineWidth',2);
plot(x_points_2,y_points_2+ mu_2(2),'g-','LineWidth',2);
plot(x_points_2,-y_points_2+ mu_2(2),'g-','LineWidth',2);
axis equal;
xlabel('x_1'); ylabel('x_2');




%%

p1 = prior_1;
p2 = prior_2;



x_points = [-3:0.1:5];
% R = -4 - 4/5 *log(16) + 8/5 *log(p2) - 8/5 *log(p1)
% R = 25/3 + 4/3 *log(16) - 8/3 *log(p2) + 8/3 *log(p1) - 1
R = 8/3 * (25/8 + 1/2 *log(16) - log(p2) + log(p1)) + 25/9
y_points = sqrt((R) - (x_points + 5/3).^2);

% R2 =  - (25/8 + 1/2 * log(16) - log(p2) + log(p1))
% y_points = sqrt(R2 + 5/4 .* x_points + 3/8 .* x_points.^2)*8/3;
% figure;
hold on;
plot(x_points,y_points,'-k')
hold on;
plot(x_points,-y_points,'-k')
axis equal
