% Copyright (c) 2007 Antti Ukkonen <Antti.Ukkonen@tkk.fi>
%
% Permission to use, copy, modify, and distribute this software for any
% purpose with or without fee is hereby granted, provided that the above
% copyright notice and this permission notice appear in all copies.
%
% THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
% WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
% MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
% ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
% WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
% ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
% OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
%
% This code was provided as supplementary material for the course
% Machine Learning: Basic Principles given at Helsinki University of
% Technology. It probably shouldn't be used AS IS for any other
% purpose than studying LDA and Matlab.

function [Z, W] = simple_lda( X, c );

d  = size(X, 2);
k  = max(c);
Sw = zeros(d);
Sb = zeros(d);

for i=1:k
  Xi = X( c==i, : );
  mi(:,i) = mean(Xi)';
  N(i) = size(Xi, 1);
end
m = (sum(mi')./k)';
for i=1:k
  Xi = X(c == i, :);
  Xi = Xi - ones(N(i), 1) * mi(:,i)';
  Sw = Sw + Xi'*Xi;
  Sb = Sb + N(i)*((mi(:,i)-m)*(mi(:,i)-m)');
end
[W, e] = eig( Sw^(-1) * Sb );
e = diag(e);
[e, idx] = sort(e, 'descend');
W = W(:,idx);
Z = X*W;