function [pc, eigenVord, eigenDord, loadings] = pca(X, c)
%PCA Conducts a principal components analysis on a data set by using the covariance matrix.
%
% USE:
%
% [pc, U, Lambda, loadings] = pca(X, cov_option)
%
% OUTPUTS:
% pc is the matrix of principal components.  
% U is the eigenvectors of the covariance matrix of X
% Lambda is the eigenvalues of the covariance matrix of X
% loadings is the correlation between components and original variables.
%          Rencher (1995) does not recommend the use of this quantity in interpreting results.
%          loadings = e*sqrt(lambda)/sqrt(var(Xk))
%          Each column indicates a principal component and each row indicates an original variable.
%          So, the i-th row and j-th column of a loadings matrix is the correlation between the
%          i-th variable and j-th principal component.
%
% INPUTS:
% X should be n x p matrix, where p is the number of variables.
% set cov_option = 0 if the correlation matrix of X is used instead of covariance matrix of X.
%           default is cov_option = 1.

% Tomo Eguchi
% 5 February 2002
% 25 February 2002 (added correlation matrix option)
% 27 February 2002 (added loadings output)

% test data set:
%means = repmat([3,5,2,9], 25, 1); var = repmat([4,2,4,10], 25, 1);
%X = normrnd(means, var, 25, 4);

if nargin == 1,
   c = 1;
end

[nrX, ncX] = size(X);

% Make sure X has full column rank
rankX = rank(X);
if (rankX ~= ncX),
   fprintf(1, 'The input matrix does not have full column rank.\n');
   fprintf(1, '# columns: %d, rank(X): %d --- pca.m\n', ncX, rankX);
   
end

% compute the covariance matrix of X:
%covX = (X - repmat(mean(X), nrX, 1))' * (X - repmat(mean(X), nrX, 1))./(nrX - 1);
if c == 1,
   covX = cov(X);
else
   covX = corrcoef(X);
end

[nrCX, ncCX] = size(covX);

% find eigenvalues and eigenvectors of covX:
[eigenVord, eigenDord] = eigOrder(covX);

% Compute principal components:
if c == 1,
   pc = demean(X) * eigenVord;
else
   pc = (demean(X)./(repmat(std(X), nrX, 1))) * eigenVord;
end

% compute factor loadings
eValuesMat = repmat(sqrt(sum(eigenDord)), [nrCX, 1]);    % eigen values repeated nrCX times
varMatrix = repmat(sqrt(sum(diag(diag(covX)))), [nrCX, 1]);  % variances repaeted nrCX times
loadings = (eigenVord .* eValuesMat)./varMatrix;