clear all, close all, clc

load allFaces.mat

% We use the first 36 people for training data
trainingFaces = faces(:,1:sum(nfaces(1:36)));
avgFace = mean(trainingFaces,2);  % size n*m by 1;

% Compute eigenfaces on mean-subtracted training data
X = trainingFaces-avgFace*ones(1,size(trainingFaces,2));
[U,S,V] = svd(X,'econ');

figure(1)
imagesc(reshape(avgFace,n,m)), colormap(gray) % Plot avg face
figure(2)
imagesc(reshape(U(:,1),n,m)), colormap(gray)  % Plot first eigenface

%% Now show eigenface reconstruction of image that was omitted from test set

testFace = faces(:,1+sum(nfaces(1:36))); % First face of person 37
ff = 3;
figure(ff)
imagesc(reshape(testFace,n,m)), colormap(gray)

testFaceMS = testFace - avgFace;
ff = 4;
for r=[25 50 100 200 400 800 1600] 
    reconFace = avgFace + (U(:,1:r)*(U(:,1:r)'*testFaceMS));
    figure(ff)
    imagesc(reshape(reconFace,n,m)), colormap(gray)
    ff = ff + 1;
end


%% Project person 2 and 7 onto PC5 and PC6

P1num = 2;  % Person number 2
P2num = 7;  % Person number 7

P1 = faces(:,1+sum(nfaces(1:P1num-1)):sum(nfaces(1:P1num)));
P2 = faces(:,1+sum(nfaces(1:P2num-1)):sum(nfaces(1:P2num)));

P1 = P1 - avgFace*ones(1,size(P1,2));
P2 = P2 - avgFace*ones(1,size(P2,2));

PCAmodes = [5 6];    % Project onto PCA modes 5 and 6
PCACoordsP1 = U(:,PCAmodes)'*P1;
PCACoordsP2 = U(:,PCAmodes)'*P2;

figure(100)
plot(PCACoordsP1(1,:),PCACoordsP1(2,:),'kd')
hold on
plot(PCACoordsP2(1,:),PCACoordsP2(2,:),'r^')