function [Y,Uk,Vk,Pk,reconErr,gain] = PCA(X,k,recenter,rescale,useGram)
%% Input:
% - X (d by N doubles matrix), d dimensions and N examples.
% - k (integer) the number of principal components to return
% - recenter (logical) whether the data should be centered
% - rescale (logical) whether the data should be rescaled
% - useGram (logical) whether to use Gram Matrix or not
%
% Output:
% - Y  (d by N doubles matrix) d dimensions and N examples.
% - Uk (d by k doubles matrix) k first principal components.
% - Vk (k by k doubles matrix) k largest eigenvalues of (1/N)XX^T
% - Pk (k by N doubles matrix) projection of data points unto the PCs.
% - reconErr (double) the reconstruction error, obtained from e.values.
% - gain (double) the compression gain.


[d,N]  = size(X);
gain = d*N/(k*(d+N));
%Normalize data as needed
if(recenter)
    Y = X - repmat(mean(X,2),1,N);
end

if(rescale)
    Y = diag(1./sqrt(var(X,[],2)))*Y;
end

if(useGram)
    [U,V] = eig(Y'*Y);
    U = Y*U*(diag(1./sqrt(diag(V))));
    V = V/N;
    [sortedEvalues,index] = sort(diag(V));
    Uk = U(:,index((end-k+1):end));
    Vk = V(index((end-k+1):end),index((end-k+1):end));
else
    [U,V] = eig((Y*Y')/N);
    [sortedEvalues,index] = sort(diag(V));
    Uk = U(:,index((d-k+1):d));
    Vk = V(index((d-k+1):d),index((d-k+1):d));
end

Pk = Uk'*Y;
Y = Uk*Pk;
reconErr = sum(sortedEvalues(1:(end-k)));
%Note the reconstruction error is incorrect if rescaling is performed.
%The correct formula was not covered in class.

%Denormalize, in reverse order
if(rescale)
    Y = diag(sqrt(var(X,[],2)))*Y;
end
    
if(recenter)
    Y = Y + repmat(mean(X,2),1,N);
end

end


