function [Xm, P, b, pcaDat] = lda(X, N)
% function [Xm, P, b, pcaDat] = lda(X, N)
%
% Note that in order for LDA to perfrom we need to have at least t + c
% samples so that the inter-variance matrix does not become singular. Here
% t denotes the dimension of our original data, and c denotes the number of
% classes. So if we have 3 groups of faces where each faces is given by 100
% points we need at least 103 faces. This in practice is unfeasible so we
% generally project down onto a lower dimensional space using PCA.
%
% Inputs:
%  X - a cell array of groups, each element of the cell array is a matrix
%      where every row is an observation and every column is a variable
%  N - the number of dimensions to project down to (default = 2)
%
% Outputs:
%  Xm - the mean of the projected data
%  P - the projection matrix
%  b - the variance covered along those axes
%  pcaDat - the full P and b info.
%
% Example:
%
% figure(1); cla; hold on;
% A = randn(100,3);
% B = randn(100,3);
% B(:,1:2) = B(:,1:2)+3;
% plot3(A(:,1), A(:,2), A(:,3), 'o');
% plot3(B(:,1), B(:,2), B(:,3), 'ro');
% X{1} = A; X{2} = B;
% [Xm, P, b, pcaDat] = lda(X, 2);
% a = (P*A')';
% b = (P*B')';
% figure(2); clf; hold on;
% plot(a(:,1), a(:,2), 'o');
% plot(b(:,1), b(:,2), 'ro');
%
% Dr. A. I. Hanna (2006).
if nargin<2
    N = 2;
end
for i=1:length(X)
    if i==1
        Sw = cov(X{i});
    else
        Sw = Sw + cov(X{i});
    end
    mean_g(i,:) = mean(X{i});
end
Sb = cov(mean_g);
[V, L] = eig(pinv(Sw)*Sb);
L = diag(L);
V = real(V);
L = real(L);
[L idx] = sort(L, 'descend'); 
V = V(:,idx);
if det(Sw)<.00001
    fprintf('Matrix Sw is near singular, watch out for complex values\n');
end
U = V(:,1:min(N, size(V, 2)));
Lu = L(1:min(N, size(V, 2)));
vr = cumsum(L) / sum(L);
Xm = (U'*mean(mean_g)')';
P = U';
b = Lu';
pcaDat.P = V';
pcaDat.b = L';
pcaDat.v = vr;