function [classification_rate] = faces(k)

fprintf('Loading data...\n');
load('ORL_32x32.mat'); % matrix with face images (fea) and labels (gnd)
load('train_test_orl.mat'); % training and test indices (trainIdx, testIdx)
fea = double(fea / 255);

% partition the data into training and test subset
fprintf('Partitioning data...\n');
n_train = size(trainIdx,1);
n_test = size(testIdx,1);
train_data = fea(trainIdx,:);
train_label = gnd(trainIdx,:);
test_data = fea(testIdx,:);
test_label = gnd(testIdx,:);

%get the mean face image
u = mean(train_data);
uspace = repmat(u,n_train,1); 

fprintf('Running PCA...\n');
components = princomp(train_data); % find principal components (use princomp function

fprintf('Projecting training data...\n');
ucomponents = components(:,1:k); % low-dim coefficients for training data (projection onto components)

train_data_pca = (train_data - uspace)*ucomponents;
fprintf('Reconstructing training data...\n');
train_data_reconstructed = uspace + train_data_pca*ucomponents';

uspace = repmat(u,n_test,1); 
fprintf('Projecting test data...\n');
test_data_pca = (test_data - uspace)*ucomponents;
fprintf('Reconstructing testing data...\n');
test_data_reconstructed = uspace + test_data_pca*ucomponents';

fprintf('Running nearest-neighbor classifier...\n');
mindist = 100000;
idx = 0;
nn_ind = [];
estimated_label = [];
for i=1:n_test
    for j=1:n_train
        temp = abs(test_data_pca(i,:) - train_data_pca(j,:));
        if(sum(temp) < mindist)
           mindist = sum(temp);
           idx = j;
        end
    end
    nn_ind(i,:) = idx;
    estimated_label(i,:) = train_label(idx,:);
    mindist = 100000;
    idx = 0;
end

classification_rate = sum(estimated_label == test_label)/n_test;

fprintf('Classification rate: %f\n', classification_rate);
j = 1;
for i=1:n_test
    if estimated_label(i) ~= test_label(i)
        blah(j,:) =  test_data(i,:);
        j = j + 1;
    end
end

for i=1:n_test
    if estimated_label(i) ~= test_label(i)
        blah(j,:) =  train_data_reconstructed(nn_ind(i),:);
        j = j + 1;
    end
end

display_faces(blah,size(blah,1)/2,2);

% figure;
% for batch = 1:10
%     clf;
%     for i = 1:12
%         test_ind = (batch-1)*12+i;
%         subplot(4,12,i);
%         imshow(reshape(test_data(test_ind,:),[32 32]),[]);
%         if i == 6
%             title('Orig. test img.');
%         end
%         subplot(4,12,i+12);
%         imshow(reshape(test_data_reconstructed(test_ind,:),[32 32]),[]);
%         if i == 6
%             title('Low-dim test img.');
%         end
%         subplot(4,12,i+24);
%         imshow(reshape(train_data_reconstructed(nn_ind(test_ind),:),[32 32]),[]);
%         if i == 6
%             title('Low-dim nearest neighbor');
%         end 
%         subplot(4,12,i+36);
%         imshow(reshape(train_data(nn_ind(test_ind),:),[32 32]),[]);
%         if i == 6
%             title('Orig. nearest neighbor');
%         end
%         if estimated_label(test_ind)~=test_label(test_ind)
%             xlabel('incorrect');
%         end
%     end
%     pause;
% end