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 = 200; %% Number of points per cluster
%% VScale = 1-3; %% Scale of the clusters

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


%%%
%%% GENERATE THE SYNTHETIC DATA ON WHICH FISHER DISCRIMINANTS
%%% WILL BE CALCULATED
%%%


%% Generate the means for the clusters - each mean is a random point
mu = 10 * (rand(C, D) - 0.5);  %% Each row of mu is the mean of a cluster
sig = [];

%% Generate the covariance matrices
for i = 1 : C,
	s = rand(D, D); %% Pick random matrix
	s = s * s'; %% Make it positive definite
	sig = [ sig; s ]; %% store it in sig matrix
end  
sig = VScale^2 * sig;

%% Generate the data
%% X has r=C*N rows (number of data points) and c=D (n. of dimensions) columns
X = genData(mu, sig, C * N);

%% 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 


%%%
%%% CALCULATE THE FISHER LINEAR DISCRIMINANTS 
%%% AND DISPLAY THE DATA ACCORDINGLY
%%%
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, %% Iterate over all clusters and plot each in a different color
  plot(X0((i-1)*N+[1:N], 1), X0((i-1)*N+[1:N], 2), mrks{i});
end    
    
line(10*[0 LDs(1, 1)], 10*[0 LDs(2, 1)]);
hold off;
axis equal;

%% Plot points in 1st and 2nd Fisher directions
figure(2); clf; title('Projection on first two FLDs'); hold on;
for i=1:C,
  clusterData = X0((i-1)*N+[1:N],:)';
  plot(LDs(:,1)' * clusterData, LDs(:,2)' * clusterData, mrks{i});
end    
hold off;

%%%
%%% Plot data in directions picked using various PCA techniques
%%%
figure(3); clf; hold on;
title('Projection on first 2 principal components of all data in all clusters');

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
hold off;

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;