%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% TDI_gaussianclas  -  Demo de clasificaci?n para el caso binario gaussiano
%                      en dos dimensiones.  La demo muestra gr?ficas 3D
%                      de las verosimilitudes de cada hip?tesis, as? como la
%                      frontera de clasificaci?n obtenida por el clasificador
%                      de m?nima probabilidad de error (MAP)
%
% Version:             1.0 (Jerónimo Arenas García, 26/02/2009)
% Asignatura:          Teoría Moderna de la Detección y la Estimación 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear
close all

% Ajustes - Modificar estos par?metros para cambiar las verosimilitudes, 
% probabilidades a priori y opciones de visualizaci?n

example = 4; % Ajuste de las verosimilitudes y probabilidades a priori
visual = 1; % (0/1) Cambia la apariencia de color del gr?fico 3D

% Study the covariance matrices, means, and a priori probabilities for the
% different available examples.  You can also modify these settings to
% study other cases.

switch example

    case 1   %%%% Circunferencia ---  CASO s1 = s0 = s y Vi = vi I
        
        V0 = [1 0;
              0 1];
        V1 = [2 0;
              0 2];
        s0 = [0 0]; s1 = -s0;
        P0 = 0.5;
        P1 = 1 - P0;
          
    case 2  %%%  Aumentamos la probabilidad a priori P1
        
        V0 = [1 0;
              0 1];
        V1 = [2 0;
              0 2];
        s0 = [0 0]; s1 = -s0;
        P1 = 0.65;
        P0 = 1 - P1;
        
    case 3  % Caso general con s1 neq s0
        
        V0 = [5 0;%%%%%  Incrementar este n?mero para observar otros casos (=2 -> parabola; >2 -> hiperbola)
            0 1];
        V1 = [2 0; 
            0 2];
        s0 = [-.8 0]; s1 = -s0;
        P0 = 0.5;
        P1 = 1 - P0;

    case 4 % Caso lineal
        V0 = [2 .7;
            .7 2];
        V1 = V0;
        s0 = [1 -1]; s1 = -s0;
        P0 = 0.99;  %%%  Aumentar hasta 0.999 para ilustrar que se puede decidir D0 incluso en s1
        P1 = 1-P0;
        
end

eta = P0/P1; %(Minimum error probability (MAP) -> C00=C11=0; C01=C10=C)
inv_V0 = inv(V0);
inv_V1 = inv(V1);
det0_sqr = det(V0)^0.5;
det1_sqr = det(V1)^0.5;

%Number of points for the representation
n_points = 200;

x1 = linspace(-5,5,n_points); x1 = x1(:);
x2 = linspace(-5,5,n_points); x2 = x2(:);
Z = zeros(n_points,n_points); % Will save the LRT output for each point
p0x = zeros(n_points,n_points); % Will save p(x|H0)
p1x = zeros(n_points,n_points); % Will save p(x|H1)

% Now, p(x|H0), p(x|H1) and the value of the LRT are calculated for points
% in grid [-5:5]x[-5:5] with n_points per dimension
for k = 1:n_points
    X0 = [x1(k)*ones(n_points,1) x2] - ones(n_points,1)*s0;
    X1 = [x1(k)*ones(n_points,1) x2] - ones(n_points,1)*s1;
    Z(:,k) = sum((X0 * inv_V0) .* X0,2) - sum((X1 * inv_V1) .* X1,2) - 2*log(eta) - log(det(V1)/det(V0));
    p0x(:,k) = (1/(2*pi*det0_sqr)) * exp(-0.5*sum((X0 * inv_V0) .* X0,2));
    p1x(:,k) = (1/(2*pi*det1_sqr)) * exp(-0.5*sum((X1 * inv_V1) .* X1,2));
end

% % Plots the output of the LRT, as well as the classification border
% figure
% H2 = imagesc(x1,x2,Z); colorbar
% hold on
% [kk, H3] = contour(x1,x2,Z,[0 0],'w');  set(H3,'LineWidth',2)
% axis square
% set(get(H2,'Parent'),'Ydir','normal')
% H1 = plot(s0(1),s0(2),'w.'); set(H1,'MarkerSize',18)
% H1 = plot(s1(1),s1(2),'w.'); set(H1,'MarkerSize',18)
% H1 = xlabel('x_1'); set(H1,'FontSize',16);
% H1 = ylabel('x_2'); set(H1,'FontSize',16);
% H1 = title('Output of the LRT test - Classification border'); set(H1,'FontSize',12);

% Plots the level curves of both likelihoods, together with the
% classification border
figure
H2 = imagesc(x1,x2,max(p0x,p1x)); colorbar
hold on
[kk , H3] = contour(x1,x2,Z,[0 0],'w'); set(H3,'LineWidth',2)
axis square
set(get(H2,'Parent'),'Ydir','normal')
H1 = plot(s0(1),s0(2),'w.'); set(H1,'MarkerSize',18)
H1 = plot(s1(1),s1(2),'w.'); set(H1,'MarkerSize',18) 
H1 = xlabel('x_1'); set(H1,'FontSize',16);
H1 = ylabel('x_2'); set(H1,'FontSize',16);
H1 = title('p(x|H_0) & p(x|H1) level curves - Classification border'); set(H1,'FontSize',12);

figure
H2 = surfc(x1,x2,p1x);
shading flat
hold on
if (visual)
    H2 = surfc(x1,x2,p0x,-p0x);
else
    H2 = surfc(x1,x2,p0x);
end
shading flat
axis square
H1 = xlabel('x_1'); set(H1,'FontSize',16);
H1 = ylabel('x_2'); set(H1,'FontSize',16);
H1 = title('p(x|H_0) & p(x|H1) 3D plot'); set(H1,'FontSize',12);

figure(1)

% figure
% H2 = imagesc(x1,x2,p0x); colorbar
% hold on
% contour(x1,x2,Z,[0 0],'k')
% axis square
% set(get(H2,'Parent'),'Ydir','normal')
% H1 = plot(s0(1),s0(2),'k.'); set(H1,'MarkerSize',18)
% H1 = plot(s1(1),s1(2),'k.'); set(H1,'MarkerSize',18) 
% 
% figure
% H2 = imagesc(x1,x2,p1x); colorbar
% hold on
% contour(x1,x2,Z,[0 0],'k')
% axis square
% set(get(H2,'Parent'),'Ydir','normal')
% H1 = plot(s0(1),s0(2),'k.'); set(H1,'MarkerSize',18)
% H1 = plot(s1(1),s1(2),'k.'); set(H1,'MarkerSize',18)