%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Demo code for ParzenWindows - 2D distribution
% S.P. Adam
% ParzenWindowsDemo.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear
clc

% The following is used to define the demo input data
randn('seed',0);
l=2;

m1 = [-1.0 0; 0 -1.0; 1.0 0; 0 1.0]';
m2 = [-1.0 -1.0; 0 0; 1.0 -1.0; -1.0 1.0; 1.0 1.0]';

[l,c1] = size(m1);
[l,c2] = size(m2);
P1 = ones(1,c1)/c1;
P2 = ones(1,c2)/c2;
s = 0.04;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N1 = 80;
for i=1:c1
    S1(:,:,i) = s*eye(l);
end
sed = 0.5;
[X1,y1] = mixt_model(m1,S1,P1,N1,sed);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N2 = 100;
for i=1:c2
    S2(:,:,i) = s*eye(l);
end
sed = 0.5;
[X2,y2] = mixt_model(m2,S2,P2,N2,sed);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
XX = [X1 X2];
yy = [ones(1,N1) -ones(1,N2)];

figure(1);
%plot(XX(1,yy==1),XX(2,yy==1),'ro',XX(1,yy==-1),XX(2,yy==-1),'bx');

% Use these limits to compute restricted distribution
% x1min = min(XX(1,:));
% x1max = max(XX(1,:));
% x2min = min(XX(2,:));
% x2max = max(XX(2,:))

axis([-2 2 -2 2]);
fig_num = 1;

% Generate the points
%axis equal
 
n1 = N1;
%n2 = N2;
n = n1; %+n2;
x1 = X1;
plot(x1(1,:), x1(2,:), 'ro');
%x2 = X2;
%hold on
%plot(x2(1,:), x2(2,:), 'b+');
points = x1;
xlim([-2 2]);
ylim([-2 2]);
 
M = 100;
xs = linspace(-2, 2, M);
ys = linspace(-2, 2, M);
 
pXY = zeros(M);
[X Y] = meshgrid(xs,ys);

% Parzen window algorithm iteration
Vn = 1/sqrt(n);
h = sqrt(Vn);
h2 = 1;
hold on
for i = 1:M
    for j= 1:M
        pXY(i,j) = 0;
        skn = 0;
        for k=1:n
            inside = (abs((xs(i)-points(1,k))/h) < h2) &&...
                     (abs((ys(j)-points(2,k))/h) < h2);
            kn = inside;
            skn = skn+kn;
        end
        pXY(i,j) = (1/(n*Vn))*skn;
    end
end

% Compute distribution at specific dataset points
pP = zeros(1,n);
for i=1:n
    skn = 0;
    for k=1:n
        inside = (abs((points(1,i)-points(1,k))/h) < h2) &&...
                 (abs((points(2,i)-points(2,k))/h) < h2);
        kn = inside;
        skn = skn+kn;
     end
     pP(1,i) = (1/(n*Vn))*skn;
end

% plot'em all
figure(2);
surf(X,Y,pXY');
colormap(copper);
title(['h = ', num2str(h)]);
xlim([-2 2]);
ylim([-2 2]);

hold on;
scatter3(points(1,:),points(2,:),pP(1,:),'r','filled');
hold off;