function P = principalComponents(X,q)
%principalComponents Principal components of a vector population.
%   P = principalComponents(X,Q) computes the principal-component
%   vectors of the vector population contained in the rows of X, a
%   matrix of size K-by-n where K (assumed to be > 1) is the number of
%   vectors and n is their dimensionality. Q, with values in the range
%   [0,n], is the number of eigenvectors used in constructing the
%   principal-components transformation matrix. P is a structure with
%   the following fields:
%
%     P.Y      K-by-Q matrix whose columns are the principal-
%              component vectors.
%     P.A      Q-by-n principal components transformation matrix
%              whose rows are the Q eigenvectors of Cx corresponding
%              to the Q largest eigenvalues. 
%     P.X      K-by-n matrix whose rows are the vectors 
%              reconstructed from the principal-component vectors. 
%              P.X and P.Y are identical if Q = n. 
%     P.ems    The mean square error incurred in using only the Q
%              eigenvectors corresponding to the largest
%              eigenvalues. P.ems is 0 if Q = n.  
%     P.Cx     The n-by-n covariance matrix of the population in X.
%     P.mx     The n-by-1 mean vector of the population in X.
%     P.Cy     The Q-by-Q covariance matrix of the population in
%              Y. The main diagonal contains the eigenvalues (in
%              descending order) corresponding to the Q
%              eigenvectors.
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based 
%   on the theoretical and practical foundations established in the 
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark 
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt

K = size(X,1);
X = double(X);

% Obtain the mean vector and covariance matrix of the vectors in X.
[P.Cx,P.mx] = covmatrix(X);
% Convert mean vector to a row vector.
P.mx = P.mx'; 

% Obtain the eigenvectors and corresponding eigenvalues of Cx.  The
% eigenvectors are the columns of n-by-n matrix V.  D is an n-by-n
% diagonal matrix whose elements along the main diagonal are the
% eigenvalues corresponding to the eigenvectors in V, so that X*V =
% D*V.
[V,D] = eig(P.Cx);

% Sort the eigenvalues in decreasing order. Rearrange the eigenvectors
% to match.
d = diag(D);
[d,idx] = sort(d);
d = flipud(d);
idx = flipud(idx);
D = diag(d);
V = V(:,idx);

% Now form the q rows of A from the first q columns of V.
P.A = V(:,1:q)';

% Compute the principal component vectors.
P.Y = P.A*(X - P.mx)'; % q-by-K matrix.

% Obtain the reconstructed vectors.
P.X = (P.A'*P.Y)' + P.mx;

% Convert P.Y to a K-by-q array and P.mx to an n-by-1 vector.
P.Y = P.Y';
P.mx = P.mx';

% The mean square error is the sum of all the eigenvalues minus the sum
% of the q largest eigenvalues.
d = diag(D);
P.ems = sum(d(q + 1:end));

% Covariance matrix of the Y's.
P.Cy = P.A*P.Cx*P.A';

