function run_experiment()

%Initialize variables
global w; global h;
w = 400; h = 450;
[X] = load_imgs('training');

%Reduce dimensionality while preserving .95 of the energy
[W mX mu] = getPCA(X,0.95,1);

%Load and recover corrupted image
[I] = load_imgs('corrupted');
recover(I,mu,W,1);

%Load and make prediction image
[I] = load_imgs('testing');
predict(I,mu,W,mX,1);



function [P] = predict(Itest,mu,W,mX, show_info)

%Subtract the mean and project
Itest = Itest(:) - mu;
Y = W' * mX;
ynew = W' * Itest;

%Compute Euclidean distance
Dist =  sqrt(sum((Y - repmat(ynew,1,size(Y,2))).^2));
[v ord] = sort(Dist);

%Closest Match
P = mX(:,ord(1)) + mu;

if show_info   
    global w; global h;
    figure;
    subplot(121); imagesc(reshape(Itest+mu,[h w]));
    colormap gray; set(gca,'XTick',[],'YTick',[]);
    axis equal; title('Testing Image');

    subplot(122);imagesc(reshape(P,[h w]));
    colormap gray; set(gca,'XTick',[],'YTick',[]);
    title('Closest Image'); axis equal;
end



function [rI] = recover(I,mu,W,show_info)
%In order to fill the missing information, we project I onto a lower dimensional 
%space and then get it back
ynew = W' * (I - mu);
rI = W * ynew;
%We add the mean back
rI = rI + mu;

if show_info
    global w; global h;
    figure;
    subplot(121); imagesc(reshape(I,[h w]));
    colormap gray; set(gca,'XTick',[],'YTick',[]);
    axis equal; title('Before');

    subplot(122);imagesc(reshape(rI,[h w]));
    colormap gray; set(gca,'XTick',[],'YTick',[]);
    title('After'); axis equal;
end




%PCA
function [W mX mu] = getPCA(X, energy, show_info)
n = size(X,2);

%-Subtract the mean
mu = mean(X,2);
mX = X - repmat(mu,1,n);
%-There is no need to make unit variance
%because all of the features have the same scale
%-Compute projection matrix
%(eigenvalues are already sorted)
[U,S,V] = svd(mX,0);

%As we saw in class, we can the eigenvectors based on their energy
%(i.e. normalized value of its eigenvalue)
D = diag(S);
tmp = D./sum(D);
good_idx = find(cumsum(tmp)<=energy);

W = U(:,good_idx);

if show_info
    global w; global h;
    %Show first 3 eigenfaces
    figure;
    for i = 1:3
        subplot(1,3,i);
        imagesc(reshape(W(:,i),[h w]));
        title(sprintf('PC %d',i)); colormap gray;
        set(gca,'XTick',[],'YTick',[]); axis equal
    end

    %Show energy for all eigenvalues
    figure;
    plot(tmp);
    hold on;
    plot([length(good_idx) length(good_idx)],[0 max(tmp)],'r');
    xlim([1 n]);
    ylabel('Energy');
    xlabel('Eigenvalues');
end



