function p4(nfaces, sfaces)

OLDCOLORMAP = COLORMAP;
COLORMAP(repmat([0:0.01:1]',[1 3]));

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, testfaces),
%% subtracts the mean face from the set of faces
   meanface = mean(trainfaces,2);
   faces = (testfaces - repmat(meanface,[1 size(testfaces,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,trainfaces),0);   
   eigencumsum = cumsum(diag(s));
   numvecsneeded = find( (eigencumsum > PCApercentage*(eigencumsum(length(eigencumsum)))) );
   minvecsneeded = numvecsneeded(1);

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

function performance = testPCAClassify(nfaces, sfaces, pcapercentage) 
% nfaces is our training set, sfaces our test set

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

  % project nfaces on to our reduced basis: 
  projectNfaces = basis'*minusMean(nfaces,nfaces);
  size(projectNfaces)
  projectSfaces = basis'*minusMean(sfaces,nfaces);
  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;

  NUMPROJECTED_CORRECT = length(PROJECTED_CORRECT);
  performance = NUMPROJECTED_CORRECT/numtestfaces;
end;

% == part a == 

nfaces = standardize(nfaces);
sfaces = standardize(sfaces);

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

% == part b == 
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)

%% part 4d == now try all the various percentages ==
performance = [];
percentage = [0.9:-0.1:0.1];
for p=percentage,
  performance = [performance testPCAClassify(nfaces,sfaces,p)];
end;
figure(3);
performance
plot(percentage,performance);
title('part 4d) performance versus percent of PCA data preserved');

%% == 5a ==

figure(11);
faces = [nfaces sfaces];
meanface = mean(faces,2);
subplot(1,4,1,'replace');
showFace(meanface);
title('5a) Mean face of entire set (smiling + neutral)');
[basis evals] = getPCABasis(faces,.9);
for eigenface_i=1:3,
  subplot(1,4,1+eigenface_i,'replace');
  showFace(basis(:,eigenface_i));
  title(sprintf('eigenface %d',eigenface_i));
end;

%% == 5b ==

function reconstructionSearch(image, midpoint, trainset),
% pass midpoint in as row, column !

   SEARCH_WINDOWSIZE = [80 80];
   WINDOWSIZE = [51 43];
   NUMBASES = 10;

   [basis evals] = getPCABasis(trainset,.99);
   basis = basis(:,[1:NUMBASES]);
   evals = basis(:,[1:NUMBASES]);

   ERRORS_BY_POSITION = [];
  
   for ii=midpoint(1)-SEARCH_WINDOWSIZE(1)/2:midpoint(1)+SEARCH_WINDOWSIZE(1)/2,
     for jj=midpoint(2)-SEARCH_WINDOWSIZE(2)/2:midpoint(2)+SEARCH_WINDOWSIZE(2)/2,
       % cut out a small snippet from image
       window = image([ii-floor(WINDOWSIZE(1)/2):ii+floor(WINDOWSIZE(1)/2)], ...
	              [jj-floor(WINDOWSIZE(2)/2):jj+floor(WINDOWSIZE(2)/2)] );

       window = standardize(window(:));
       %% project it on our bases, and reconstruct
       %% do we want to first subtract the mean?
       %% subtract mean
       window = window - repmat(mean(window),[length(window) 1]);

       coeffs = basis'*window;
       % how do we reconstruct?
       reconstruction = sum(diag(coeffs)*basis',1)';
%        figure(15);
%        subplot(1,2,1);
%        showFace(window);
%        subplot(1,2,2);
%        showFace(reconstruction);
       DIFFERENCE = sum(abs(reconstruction-window));
       ERRORS_BY_POSITION = [ERRORS_BY_POSITION DIFFERENCE];
%       pause;
       
     end
   end
   figure(14);
   [x y] = meshgrid([1:81],[1:81]);
   sizeerrors = size(ERRORS_BY_POSITION)
   plot3(x,y,reshape(ERRORS_BY_POSITION, [81 81]));
   figure(15);
   imagesc(reshape(ERRORS_BY_POSITION, [81 81]));
   [minv minidx] = min(ERRORS_BY_POSITION);
   [MINX MINY] = sub2ind(minidx,81,81)
   
end;

figure(12);
mad = im2double(imread('data/mad.png'));
imshow(mad);

% empirically obtained by using ginput on mad
face1 = [148 127; 207 172]'; 
face2 = [392 131; 443 179]';
midpointf1 = round(face1(:,1)+(face1(:,2)-face1(:,1))/2.0);
midpointf2 = round(face2(:,1)+(face2(:,2)-face2(:,1))/2.0);

% now find 80x80 windows around the midpoint
figure(13);
windowf1 = mad( [midpointf1(2)-40:midpointf1(2)+40], ...
    [midpointf1(1)-40:midpointf1(1)+40] );
subplot(1,2,1); imshow(windowf1);

windowf2 = mad( [midpointf2(2)-40:midpointf2(2)+40], ...
    [midpointf2(1)-40:midpointf2(1)+40] );
subplot(1,2,2); imshow(windowf2);

reconstructionSearch(mad, flipdim(midpointf1,1), faces);


end