%Does PCA and Fisher's Linear Discriminant (FLD)
clear all; hold off; close all;
N = 100;
Dog = zeros(N, 2); %Make fake two class data
Dog(:, 1) = 2*rand(N, 1) - 1;
Dog(:, 2) = (Dog(:, 1) + 0.25*randn(N, 1)) + 0.5;
Cat = zeros(N, 2);
Cat(:, 1) = 1 - 2*rand(N, 1) + 0.5*randn(N, 1);
Cat(:, 2) = (Cat(:, 1) - 0.4*rand(N, 1));
figure(1); plot(Dog(:, 1), Dog(:, 2), 'ro', Cat(:, 1), Cat(:, 2), 'bo');
hold on;
axis equal;
MD = mean(Dog);
MC = mean(Cat);
%Now do PCA on entire set of data
Dta(1:N, :, :) = Dog;
Dta((N+1):2*N, :, :) = Cat;
CovBetween = cov(Dta);
[PCA, EV] = eig(CovBetween);
line(2*[-PCA(1, 2), PCA(1, 2)], 2*[-PCA(2, 2), PCA(2, 2)], 'Color', 'k', 'LineWidth', 2, 'LineStyle', '-');

%Now do Fisher Linear Discriminant
CovDog = cov(Dog);
CovCat = cov(Cat);
CovWithin = (CovDog + CovCat)/2;
[FLD, FV] = eig(inv(CovWithin)*CovBetween);
NmFLD = FLD(:, 2)/norm(FLD(:, 2));
line([-NmFLD(1), NmFLD(1)], [-NmFLD(2), NmFLD(2)], 'Color', 'r', 'LineWidth', 2);
FDog = Dog*NmFLD; %Project onto decision axis
FCat = Cat*NmFLD;
MnD = mean(FDog);
MnC = mean(FCat);
if MnD < MnC; %Make sure category #1 has biggest projection
    FDog = -FDog;
    FCat = - FCat;
    MnD = -MnD;
    MnC = -MnC;
end;
Correct = zeros(1, 100);
for KK = 1:100;
    Decide(KK) = MnD + (MnC - MnD)*KK/100; % Candidate thresholds between MnD & MnC
    Correct(KK) = sum(FDog > Decide(KK)) + sum(FCat < Decide(KK));
end;
[MX, II] = max(Correct);
[MX, FF] = max(flip(Correct)); %get other end of plateau of maxima
Theta = (Decide(II) + Decide(length(Correct) - FF))/2
Vthresh = [1; -NmFLD(1)/NmFLD(2)]; %Vector for threshold separatrix
Tshift = Theta/NmFLD(2);
Vthresh = Vthresh/norm(Vthresh);
line(2*[-Vthresh(1), Vthresh(1)], 2*[-Vthresh(2) + Tshift, Vthresh(2) + Tshift], 'Color', 'r', 'LineWidth', 2, 'LineStyle', '--');
% XX = -2:0.1:2;
% YY = -(NmFLD(1)*XX + NmFLD(3)*XX.^2)/NmFLD(2) - Theta/NmFLD(2);
% plot(XX, YY, 'k', 'LineWidth', 2);
NumCorrect = sum((FDog > Theta) + (FCat < Theta))
Errors = sum((FDog < Theta) + (FCat > Theta))
% figure(2); plot(Dog(:, 1), Dog(:, 2), 'ko', Cat(:, 1), Cat(:, 2), 'ko');