%solutions to project 3
clear all

%1
%load data
load X1; load X2;

%1a
N1 = size(X1,2)
N2 = size(X2,2)
m1 = mean(X1,2)
m2 = mean(X2,2)

%1b
Sw = (X1-repmat(m1,1,N1))*(X1-repmat(m1,1,N1))' + ...
    (X2-repmat(m2,1,N2))*(X2-repmat(m2,1,N2))';
w = Sw\(m2-m1);
w = w/norm(w)

%1c
Z1 = w'*X1;
Z2 = w'*X2;
mu1 = mean(Z1)
mu2 = mean(Z2)
var1 = sum( (Z1-mu1).^2 )/N1
var2 = sum( (Z2-mu2).^2 )/N2

%1d
%find intersection of Gaussians between mu1 and mu2
A = -var1 + var2;
B = 2*(-mu1*var2 + mu2*var1);
C = var1*var2*log(var1/var2) + var2*mu1^2 - var1*mu2^2;
x = (-B + sqrt(B^2 - 4*A*C))/(2*A);
w0 = -x

%1e
correct1 = sum(w'*X1+w0 <= 0);
fraction_correct1 = correct1/N1
correct2 = sum(w'*X2+w0 > 0);
fraction_correct2 = correct2/N2


%2
clear all
%load data
load digits
N3 = size(X3image,3);
N5 = size(X5image,3);
M3 = size(Y3image,3);
M5 = size(Y5image,3);
r = 28; %28 by 28 images

%2a
%try using the entire images as high dimensional features
X3 = reshape(X3image,r^2,N3);
X5 = reshape(X5image,r^2,N5);
Y3 = reshape(Y3image,r^2,M3);
Y5 = reshape(Y5image,r^2,M5);

%2b
%compute Fisher's discriminant to get w
m3 = mean(X3,2);
m5 = mean(X5,2);

Sw = (X3-repmat(m3,1,N3))*(X3-repmat(m3,1,N3))' + (X5-repmat(m5,1,N5))*(X5-repmat(m5,1,N5))';
w = pinv(Sw)*(m5-m3);
w = w/norm(w);

%choose w0 to minimize classification error
Z3 = w'*X3;
Z5 = w'*X5;
mu3 = mean(Z3);
mu5 = mean(Z5);
var3 = sum( (Z3-mu3).^2 )/N3;
var5 = sum( (Z5-mu5).^2 )/N5;

%find intersection of Gaussians between mu1 and mu2
A = -var3 + var5;
B = 2*(-mu3*var5 + mu5*var3);
C = var3*var5*log(var3/var5) + var5*mu3^2 - var3*mu5^2;
x = (-B + sqrt(B^2 - 4*A*C))/(2*A);
w0 = -x;

%2c
%classify Y3 and Y5
Y3fit = w'*Y3+w0>=0;
Y5fit = w'*Y5+w0>=0;

%2d
%make confusion matrix
confusion(1,1) = sum(Y3fit==0);
confusion(1,2) = sum(Y3fit==1);
confusion(2,1) = sum(Y5fit==0);
confusion(2,2) = sum(Y5fit==1);
confusion
fraction_correct3 = confusion(1,1)/M3
fraction_correct5 = confusion(2,2)/M5

%2e
%compare to random classifier
RY3fit = round(rand(1,M3))>=.5;
RY5fit = round(rand(1,M5))>=.5;
Rconfusion(1,1) = sum(RY3fit==0);
Rconfusion(1,2) = sum(RY3fit==1);
Rconfusion(2,1) = sum(RY5fit==0);
Rconfusion(2,2) = sum(RY5fit==1);
Rconfusion


%brute force approach worked well, but gave no insight about good features
%try a 2D feature vector consisting of central moments instead
X3 = zeros(2,N3);
X5 = zeros(2,N5);
Y3 = zeros(2,M3);
Y5 = zeros(2,M5);
[J,I] = meshgrid(1:28,1:28); %coords for 28 by 28 images

%feature 1
px = 2;
py = 1;
for i = 1:N3
    area = sum(sum(X3image(:,:,i)));
    mx = sum(sum(J.*X3image(:,:,i)))/area;
    my = sum(sum(I.*X3image(:,:,i)))/area;
    X3(1,i) = sum(sum(((J-mx).^px).*((I-my).^py).*X3image(:,:,i)))/area;
end
for i = 1:N5
    area = sum(sum(X5image(:,:,i)));
    mx = sum(sum(J.*X5image(:,:,i)))/area;
    my = sum(sum(I.*X5image(:,:,i)))/area;
    X5(1,i) = sum(sum(((J-mx).^px).*((I-my).^py).*X5image(:,:,i)))/area;
end
for i = 1:M3
    area = sum(sum(X3image(:,:,i)));
    mx = sum(sum(J.*X3image(:,:,i)))/area;
    my = sum(sum(I.*X3image(:,:,i)))/area;
    Y3(1,i) = sum(sum(((J-mx).^px).*((I-my).^py).*Y3image(:,:,i)))/area;
end
for i = 1:M5
    area = sum(sum(X5image(:,:,i)));
    mx = sum(sum(J.*X5image(:,:,i)))/area;
    my = sum(sum(I.*X5image(:,:,i)))/area;
    Y5(1,i) = sum(sum(((J-mx).^px).*((I-my).^py).*Y5image(:,:,i)))/area;
end

%feature 2
px = 1;
py = 2;
for i = 1:N3
    area = sum(sum(X3image(:,:,i)));
    mx = sum(sum(J.*X3image(:,:,i)))/area;
    my = sum(sum(I.*X3image(:,:,i)))/area;
    X3(2,i) = sum(sum(((J-mx).^px).*((I-my).^py).*X3image(:,:,i)))/area;
end
for i = 1:N5
    area = sum(sum(X5image(:,:,i)));
    mx = sum(sum(J.*X5image(:,:,i)))/area;
    my = sum(sum(I.*X5image(:,:,i)))/area;
    X5(2,i) = sum(sum(((J-mx).^px).*((I-my).^py).*X5image(:,:,i)))/area;
end
for i = 1:M3
    area = sum(sum(X3image(:,:,i)));
    mx = sum(sum(J.*X3image(:,:,i)))/area;
    my = sum(sum(I.*X3image(:,:,i)))/area;
    Y3(2,i) = sum(sum(((J-mx).^px).*((I-my).^py).*Y3image(:,:,i)))/area;
end
for i = 1:M5
    area = sum(sum(X5image(:,:,i)));
    mx = sum(sum(J.*X5image(:,:,i)))/area;
    my = sum(sum(I.*X5image(:,:,i)))/area;
    Y5(2,i) = sum(sum(((J-mx).^px).*((I-my).^py).*Y5image(:,:,i)))/area;
end

%scatter plot
figure(1)
scatter(X3(1,:),X3(2,:),'b.');
hold on
scatter(X5(1,:),X5(2,:),'r.');

%compute Fisher's discriminant to get w
m3 = mean(X3,2);
m5 = mean(X5,2);

Sw = (X3-repmat(m3,1,N3))*(X3-repmat(m3,1,N3))' + (X5-repmat(m5,1,N5))*(X5-repmat(m5,1,N5))';
w = Sw\(m5-m3);
w = w/norm(w);

%choose w0 to minimize classification error
Z3 = w'*X3;
Z5 = w'*X5;
mu3 = mean(Z3);
mu5 = mean(Z5);
var3 = sum( (Z3-mu3).^2 )/N3;
var5 = sum( (Z5-mu5).^2 )/N5;

%find intersection of Gaussians between mu1 and mu2
A = -var3 + var5;
B = 2*(-mu3*var5 + mu5*var3);
C = var3*var5*log(var3/var5) + var5*mu3^2 - var3*mu5^2;
x = (-B + sqrt(B^2 - 4*A*C))/(2*A);
w0 = -x;

%classify Y3 and Y5
Y3fit = w'*Y3+w0>=0;
Y5fit = w'*Y5+w0>=0;

%make confusion matrix
confusion(1,1) = sum(Y3fit==0);
confusion(1,2) = sum(Y3fit==1);
confusion(2,1) = sum(Y5fit==0);
confusion(2,2) = sum(Y5fit==1);
confusion
fraction_correct3 = confusion(1,1)/M3
fraction_correct5 = confusion(2,2)/M5