function p4(nfaces, sfaces)

function f = standardize(faces) 
   mus = mean(faces,1);
   f = faces-repmat(mus,[size(faces,1) 1]);
   variances = diag(cov(f))';
   f = f./sqrt(repmat(variances,[size(faces,1) 1]));
end

function showFace(facevec)
  imagesc(reshape(facevec, [51 43]));
end;

function [score,idx] = recognize(face, candidates)

   fsub = repmat(face,[1 size(candidates,2)]);
   diff = fsub - candidates;
   dist = sqrt(sum(abs(diff).^2,1));

   [score,idx] = min(dist);

end;

function faces = minusMean(trainfaces),
%% subtracts the mean face from the set of faces
   meanface = mean(trainfaces,2);
   faces = (trainfaces - repmat(meanface,[1 size(trainfaces,2)]));
end;

function [basis, evals] = getPCABasis(trainfaces, PCApercentage)
%% performs SVD decomposition on training set of faces
%% given you want to preserve at least PCApercentage of the 
%% amount of original data

   [u s v] = svd(minusMean(trainfaces),0);   
   eigencumsum = cumsum(diag(s));
   figure(10);
   plot(eigencumsum);
   numvecsneeded = find( (eigencumsum > PCApercentage*(eigencumsum(length(eigencumsum)))) );
   minvecsneeded = numvecsneeded(1);

   basis = u(:,[1:minvecsneeded]);
   evals = s(:,[1:minvecsneeded]);
end;

function performance = testPCAClassify(trainfaces, testfaces, pcapercentage) 

  [basis evals] = getPCABasis(nfaces,pcapercentage);
  SIZE_OF_BASIS = size(basis)

  % project nfaces on to our reduced basis: 
  projectNfaces = basis'*minusMean(nfaces);
  size(projectNfaces)
  projectSfaces = basis'*minusMean(sfaces);
  size(projectSfaces)

  PROJECTED_CORRECT = [];

  % now try classifying each face
  for face_i=1:size(projectSfaces,2),
    [score,rec_idx]=recognize(projectSfaces(:,face_i),projectNfaces);
    if (face_i == rec_idx)
      PROJECTED_CORRECT = [PROJECTED_CORRECT face_i];
    end;
  end;

PROJECTED_CORRECT
NUMPROJECTED_CORRECT = size(PROJECTED_CORRECT)

performance = NUMPROJECTED_CORRECT;

end;

% == main == 

nfaces = standardize(nfaces);
%showFace(nfaces(:,1));
%min(nfaces(:,1))
%hist(nfaces(:,2))

% == part b == 
sfaces = standardize(sfaces);
CORRECT = [];
for face_i=1:size(sfaces,2),
  [score,rec_idx]=recognize(sfaces(:,face_i),nfaces);
  if (face_i == rec_idx)
    CORRECT = [CORRECT face_i];
  end;
end;

CORRECT
size(CORRECT)

% == part c ==

figure(1);
meanface = mean(nfaces,2);
subplot(1,6,1);
showFace(meanface);
title('mean face');

[u s v] = svd((nfaces - repmat(meanface,[1 size(nfaces,2)])),0);
for f=1:5,
  eigenface = u(:,f);
  subplot(1,6,f+1);
  showFace(eigenface);
  title(sprintf('eigenface %d',f));
end;

figure(2);
clf;
%plot([diag(s) cumsum(diag(s))]);
%plot([diag(s)]);
plot([cumsum(diag(s))]);

eigencumsum = cumsum(diag(s));
numvecsneeded = find( (eigencumsum > ...
    0.90*(eigencumsum(length(eigencumsum)))) );
MINVECSNEEDED = numvecsneeded(1)

%% 

end