function [ Y, Z, E ] = PCA( X, d, centering )
%
% Corded by D. Kitamura (d-kitamura@nii.ac.jp) on 2 Jan. 2016.
%
% This function calculates principal component analysis of an input data.
% This function supports both PCAs with and without data centering.
%
% see also
% http://d-kitamura.sakura.ne.jp/index.html
%
% [syntax]
%   [ Y, Z, E ] = PCA( X )
%   [ Y, Z, E ] = PCA( X, d )
%   [ Y, Z, E ] = PCA( X, d, centering )
%
% [inputs]
%          X: input data ( K (variables) x N (samples) )
%          d: number of dimensions to which X is projected
%  centering: centering X first or not (true (1) or false (0), default: true)
%
% [outputs]
%          Y: output matrix (d x N)
%          Z: transformation matrix (d x K)
%          E: all eigenvectors (K x K)

% Check errors and set default values
if (nargin < 3)
    centering = true; % Default setting
end

[K,N] = size(X);
if centering
    cX = X - repmat(mean(X,2),1,N); % Data centering
    V = cX*cX.'/N; % Covariance matrix (K x K) (equivalent to V = cov(centerX.'))
else
    V = X*(X')/N; % Covariance matrix without centering
end
[P,D] = eig(V); % Eigenvalue decomposition

% Sort eigenvalues in descending order
eig_val = diag(D);
[eig_val,idx] = sort(eig_val, 'descend');
D = D(idx,idx);
P = P(:,idx);

% Pick up top d eigenvalues and their eigenvectors
D2 = D(1:d,1:d);
P2 = P(:,1:d); 

Y = P2'*X; % Output matrix
Z = P2'; % Transformation matrix
E = P'; % All eigenvectors
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% EOF %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%