function [ LDs ] = fisherLD(X, y)
%% FISHERLD   Calculates the fisher linear discriminants for given labelled data
%%            The data is in D dimensions. There are G groups.
%%
%% [ LDs ] = fisherLD( X, y )
%%
%% LDs - DxG Matrix each column is a linear discriminant, ordered by importance
%% X   - DxN Matrix. Each column is a d-dimensional datum.
%% y   - 1xN Vector. Each entry is the group label of the corresponding datum.


[D, N] = size(X);  %% Get n. of dimensions and n. of sample points from X


%% Count the number of groups from the label vector
%% The labels may be arbitrary numbers and not in order
yS = sort(y); dyS = diff(yS);
g = sum(dyS > 0) + 1;


%% Construct sorted list of group labels
groupLabels = zeros(g, 1);
groupLabels(1) = yS(1);
groupLabels(2 : g) = yS(find(dyS) + 1);

Ng = zeros(g, 1);   %% vector to contain number of sample points per group
yRelabeled = zeros(size(y)); %% vector to substitute y containing simpler labels

for i = 1 : g,
	idx = find(y == groupLabels(i));
	Ng(i) = length(idx);  %% Calculate the n. of sample points for each group
	yRelabeled(idx) = i * ones(Ng(i), 1); %% Generate label vector with 1:g labels
end

%% Build the G matrix - See Ripley's book page 93
G = zeros(N, g);
for i = 1:N,
  G(i, yRelabeled(i)) = 1;
end
	

T = diag(sqrt(N ./ Ng));

% Make X zeromean
mn = mean(X')';
X = X - mn * ones(1, N);

% Collect Group means
[U L V] = svd(X');

U = U(:, 1 : D);
L = L(1 : D, :);

S = sqrt(N) * V * inv(L);
Xrs = (X' * S)';
M = inv(G' * G) * G' * Xrs';

[UU LL VV] = svd(inv(T) * M);
LDs = S * VV;

return

%keyboard;

% tests

Xrs * Xrs'

(G * M)' * (G * M)
N * VV * LL .^ 2 * VV'

T * G' * G * T

M
inv(G' * G) * G' * Xrs'
1/N * T.^2 * G' * Xrs'




XRS = X' * S * VV';

(XRS - G * M)' * (XRS - G * M)


%A = S' * X * X' * S;

%(V' * invL)' * (U * L * V