function [data, transfm, origmean, origsd, evs] = ...
    performPCA(x, newdim, optional_origmean, optional_origsd)
% Reduces the dimension of a data set using Principal Component Analysis
% (PCA). Before PCA the data is centralized to have zero mean, and then
% divided by standard deviation to force unit variance. The whole information
% necessary for transformation with original mean and standard
% deviation is returned.
%
%  There are some different ways to call this:
%
%     [y,tM,oM,oSD] = performPCA(x, newdim) 
%         Yields a PCA projection of the columns of x as an newdim-by-N
%         matrix; also returns necessary information for repeating the 
%         transformation.
%
%     [y,tM,oM,oSD] = performPCA(x) 
%         Yields a PCA projection of the columns of x with original
%         dimension; also returns necessary information for repeating the 
%         transformation.
%
%     y = performPCA(x, tM, oM, oSD)
%         re-uses a formerly found transformation.
% 
%
% We implement this on our own to avoid dependency of the Statistical or
% Neural Network Toolboxes of Mathworks, which do provide similar
% functionality.
%
% Input:
%   x      - Data matrix; columns are observations, rows are variables
%   newdim - To how many dimensions the input is to be dropped.
%
% Output:
%   data     = transfm' * ((x-REP(origmean))./REP(origsd))
%   transfm  = PCA transformation matrix based on the eigenvectors of the
%              covariance matrix of the normalized input dataset
%   origmean = The original coordinatewise means of the data.
%   origsd   = The original standard deviations of the data.
%
%
% Notes: 
%
%   We represent an observation by a column vector; hence use transposes
%   here because built-in Matlab functions interpret observations as rows.
%   Acknowledge this when you apply the transformation matrix. In essence,
%   you'd do 
%
%       transfm' * ((yourdata - origmean)*origsd) 
%

% By default, create the transform but don't reduce dimension:
if nargin==1
    newdim = size(x,1);
end

% Optionally re-use a transformation:
if nargin==4
    transfm = newdim; %Naming oddity due to DIY-overloading.
    data = bsxfun(@minus,x,optional_origmean);
    data = transfm' * (bsxfun(@rdivide,data,optional_origsd));
    return
end

% Make sure we have zero coordinatewise means:
origmean = mean(x, 2);
x = x - repmat(origmean, 1, size(x,2));

% Make sure we have unit variance (except for constant variables)
origsd = std(x,0,2);
constants = (origsd == 0);
if sum(constants) > 0
    warning('Some variables are constant (%d out of %d)', sum(constants),numel(constants));
    origsd(constants) = 1; % Divide by one -> transformed to 0.
end
x = bsxfun(@rdivide,x,origsd);

% Reduce number of dimensions using Principal Component Analysis:
c = cov(x'); 
[u, lambda] = eig(c);

% Sort eigenvalues (at least Octave won't do it for us...)
[evs, i] = sort(abs(diag(lambda)),1,'descend');
transfm = u(:,i(1:newdim));

data = (x'*transfm)';

