%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function fda()
%   z.li, 04-12-2004
%   fisher discriminant analysis
% function dependency:
%   - n/a
% input:
%   K - number of  classes in data
%   x - data: d x n
%   y - label: 1 x n
% output:
%   A  - the transform: K x d
%   ev  - eigen values
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%function [A,ev]=fda(K, x, y)
function [A,ev]=fda(K, x, y)
 dbg = 'n';
 if dbg == 'y'
   K = 2; m=16;
   styl = ['+', '.', 'x'];
   axis([0 100 0 100]); hold on;
   x = zeros(2, K*m);
   y = zeros(1, K*m);
   for j=1:K
       fprintf('\n new class: ');
       for t=1:m
          [px py]=ginput(1);
          plot(px, py, styl(j));
          fprintf('%c', styl(j));
          x(:,(j-1)*m+t) = [px py]';
          y((j-1)*m+t) = j;
       end
   end
 end

 % const
 dbgPlot = 'n';

 % collect labeling
 [d, N] = size(x);          % sample dimension and total number
 M = mean(x')';             % sample mean for all
 mj = zeros(d, K);          % class means
 nj = zeros(1, K);          % number of samples per class

 % compute Sw - intra class scatter: dxd
 Sw = zeros(d,d);
 % process all classes
 for j=1:K
   indx =[];indx = find(y==j);
   nj(j) =length(indx);
   xj=zeros(d, nj(j)); xj = x(:,indx);
   % compute sample mean
   mj(:,j) = mean(xj')';
   for t=1:nj(j)
      Sw = Sw + (xj(:,t)-mj(:,j))*(xj(:,t)-mj(:,j))';
   end
 end

 % compute Sb - inter class scatter: dxd
 Sb = zeros(d, d);
 for j=1:K
   Sb = Sb + nj(j)*(mj(j)-M)*(mj(j)-M)';
 end

 % find fisher discriminant subspace
 [v,d]=eig(Sb, Sw);
  
 A = v; ev = diag(d);

 % dbg plot
 if dbgPlot == 'y'
   hold on;
   for j=1:K
      plot(mj(1,j), mj(2,j), 'or');
   end
   plot(M(1), M(2), 'or'); plot(M(1), M(2), '+r');
   dplot=800;
   px = [M(1)-dplot*A(1,1), M(1)+dplot*A(1,1)];
   py = [M(2)-dplot*A(2,1), M(2)+dplot*A(2,1)];
   plot(px, py, '.', px,py, '-k');
   
   px = [M(1)-dplot*A(1,2), M(1)+dplot*A(1,2)];
   py = [M(2)-dplot*A(2,2), M(2)+dplot*A(2,2)];
   plot(px, py, '.', px,py, '-r');

   figure; hold on;
   styl = ['+', '.', 'x'];
   for j=1:K
     indx =[];indx = find(y==j);
     xj=zeros(d, nj(j)); xj = x(:,indx);
     zj=A(:,1:2)*xj;
     plot(zj(1,:), zj(2,:), styl(j));
   end
 end