clear all, close all, clc

load /Users/stefanliu/Downloads/allFaces.mat

% Define the number of individuals and images per individual for training
numIndividuals = 38;
imagesPerIndividual = 50;

% Preparing training and testing sets
trainingFaces = [];
testingFaces = [];
labelTest = [];
for i = 1:numIndividuals
    startIndex = sum(nfaces(1:i-1)) + 1;
    endIndex = startIndex + nfaces(i) - 1;
    trainingFaces = [trainingFaces, faces(:, startIndex:startIndex + imagesPerIndividual - 1)];
    testingFaces = [testingFaces, faces(:, startIndex + imagesPerIndividual:endIndex)];
    labelTest(1+sum(nfaces(1:i-1))-(i-1)*50:sum(nfaces(1:i))-i*50)= i;
end

% Compute the average face from training data
avgFace = mean(trainingFaces, 2);

% Center the training data
X = trainingFaces - avgFace * ones(1, size(trainingFaces, 2));

% Perform SVD
[U, S, V] = svd(X, 'econ');

% Define dimensions to test
r = 25;

% Project training data into the top r dimensions
U_r = U(:, 1:r);
projectedTraining = U_r * U_r' * X;

% Compute mean of projections for each individual
means = zeros(size(projectedTraining, 1), numIndividuals);
for j = 1:numIndividuals
    means(:, j) = mean(projectedTraining(:, (j-1)*imagesPerIndividual+1:j*imagesPerIndividual), 2);
end

% Subtract the average face from testing data
testX = testingFaces - avgFace * ones(1, size(testingFaces, 2));
projectedTesting = U_r * U_r' * testX;

% Classify test images
Predition = [];
for j = 1:size(projectedTesting,2)
    distances = vecnorm(means-projectedTesting(:,j));
    [~,Predition(j)] = min(distances);
end
accuracy = sum(Predition == labelTest) / length(Predition) * 100;
fprintf('Classification accuracy: %.2f%%\n', accuracy);





