function fisherDemo ( D, C, N, VScale )
%% FISHERDEMO         Demonstrate Fisher linear discriminants and compare with PCA/SVD
%
% 

%%% Good values:
%D = 10; %% Dimensions
%C = 3;  %% Number of clusters
%N = 50; %% Number of points per cluster
%% VScale = 1-10; %% Scale of the clusters

if nargin < 1,
    D = 10;
    C = 3;
    N = 50;
    VScale=10;
end;

mrks = {'r.', 'g.', 'b.', 'c.', 'm.', 'k.'};  %% Colors of the markers

%% Generate the (random) means for the clusters
mu = 10 * (rand(C, D) - 0.5);
sig = [];

%% Generate the (random) covariance matrices
for i = 1 : C,
	s = rand(D, D);
	s = s * s';
	sig = [ sig; s ];
end  
sig = VScale * sig;

%% Generate the data
%% X has r=C*N rows and c=D columns
X = genData(mu, sig, C * N);
size(X)

%% Generate the cluster labels of the data
%% y is a long row vector with a label for each sample point
%% therefore its lenght is C*N
y = [];
for i = 1 : C,
	y = [y; i * ones(N, 1)];
end 

y
size(y)
%keyboard

figure(100); title('The raw data plotted along the first two dimensions');
hold on;
for c=1:C,
    for i =1:N,
        plot(X((c-1)*N+i,1),X((c-1)*N+i,2),mrks{c});
    end;
end;
hold off; xlabel('1st dimension'); ylabel('2nd dimension');


%% Calculate the Fisher linear discriminants
LDs = fisherLD(X', y);

%%%
%%% Display the results
%%%


%% Subtract the mean from the data
mn = mean(X);
X0 = X - ones(C * N, 1) * mn;

%% Plot points in 1st 2 dimensions
figure(1); clf;
title(['Canonical View (', int2str(C), ' clusters, ', int2str(D), ' dimensions, ', int2str(N), ' points)']); hold on;
for i=1 : C,
	plot(X0((i-1)*N+[1:N], 1), X0((i-1)*N+[1:N], 2), mrks{i});
end      
hold on;
line([0 LDs(1, 1)], [0 LDs(2, 1)]);
hold off;
axis equal;

figure(2); clf;  
subplot(1,2,1);
hold on;
for i=1 : C,
	%plot(LDs(:, 1)' * X0(1 + (i-1) * N:i * N, :)', LDs(:, 2)' * X0(1 + (i-1) * N:i * N, :)', mrks{i});
	plot(LDs(:,1)' * X0((i-1)*N+[1:N],:)', LDs(:,2)' * X0((i-1)*N+[1:N],:)', mrks{i});
end    
hold off;
title('Projection of training data on first two FLDs');
axis 'equal'

Xt = genData(mu, sig, C * N); %% Generate test data
Xt0 = Xt - ones(C * N, 1) * mean(Xt);
subplot(1,2,2);hold on;
for i=1 : C,
	plot(LDs(:,1)' * Xt0((i-1)*N+[1:N],:)', LDs(:,2)' * Xt0((i-1)*N+[1:N],:)', mrks{i});
end    
hold off;
title('Projection of test data on first two FLDs');
axis 'equal'




figure(3); clf; hold on;
title('PCA combined');

PCs = PCAspace(X', 2);
S = PCs.U' * (X' - PCs.mu * ones(1, C * N));

for j = 1 : C,
   plot(S(1, 1 + (j-1) * N:j * N), S(2, 1 + (j-1) * N:j * N), mrks{j});
end

figure(4); clf;

for i = 1 : C,
    subplot(ceil(C/2), 2, i); title(['Projection on PCs of clst. ', int2str(i)]);
    hold on;
    
    PCs = PCAspace(X(1 + (i-1) * N:i * N, :)', 2);
    S = PCs.U' * (X' - PCs.mu * ones(1, C * N));
    
    for j = 1 : C,
    	
        plot(S(1, 1 + (j-1) * N:j * N), S(2, 1 + (j-1) * N:j * N), mrks{j});
    end
end    


return;