function [whitesig, V, Lambda, whiteM,dewhiteM]=nic_rsdd_pca(data2,numpc)
% Do PCA 
% Inputs: 
%	data: a 2-dim matrix
%	numpc: the number of PC after PCA
%-----------------------------------------------------------
%   Copyright (C) 2009
%   Neuroimage Computing Group, State Key Laboratory of
%   Cognitive Neuroscience and Learning, Beijing normal University
%	Written by Guiwen Chen
%	Version=1.0;
%	Release=20090720;
%-----------------------------------------------------------
%

% do PCA
if isreal(data2)
    data2 = nic_icatb_remove_mean(data2, 1);
end
% end for removing mean

saveEig = 0;

% compute covariance matrix and the eigen values
[V, Lambda] = nic_icatb_v_pca(data2, 1, numpc, saveEig, 'untranspose', 'no');

% whiten matrix
[whitesig, whiteM, dewhiteM] = nic_icatb_v_whiten(data2, V, Lambda, 'untranspose');

clear data2;

whitesig = whitesig';

return


function [V, Lambda, V_orig, Lambda_orig] = nic_icatb_v_pca(data, firstEig, lastEig, saveEig, doTranspose, removeMean)
%% compute covariance matrix by removing the mean and then calculate the
% eigen values and eigen vectors

% Input: 1. data - 2D array, 1 dimension - voxels, 2 dimension - time points
% 2., 3. and 4.,  firstEig, lastEig and saveEig - scalars
% 5. doTranspose - 'transpose' - computes transpose

firstEig_t2 = firstEig; lastEig_t2 = lastEig; saveEig_t2=saveEig;

if (~exist('saveEig','var'))
    saveEig = 0;
end;

if ~exist('removeMean', 'var')
    removeMean = 'yes';
end

if (saveEig == 3),
    disp('Loading Previous EigenDecomposition 1');edit
    load eig_stuff1;
    firstEig = firstEig_t2;lastEig=lastEig_t2;saveEig=saveEig_t2;
    originalDimension = size(data, 1);
else
    originalDimension = size(data, 1);

    if exist('doTranspose', 'var')
        if strcmp(doTranspose, 'transpose')
            data = data';
        end
        originalDimension = size(data, 2);
    end

    disp(['Calculating Covariance: ' num2str(originalDimension) '^2']);

    %% Skip detrending if the removing the mean is already done
    rmMean = 0;
    if strcmpi(removeMean, 'yes')
        disp('Removing mean from the data ...');
        rmMean = 1;
    end

    %% Calculate covariance matrix
    covarianceMatrix = nic_icatb_cov(data, rmMean);

    clear data;

    %% Calculate the eigenvalues and eigenvectors of covariance matrix.
    disp('Calculating eigendecomposition');

    condCovMat = cond(covarianceMatrix);

    disp(['Condition number of covariance matrix is ', num2str(condCovMat)]);

    [V, Lambda] = eig(covarianceMatrix, 'nobalance');

    %% Sort the eigenvalues - decending.
    disp('Sorting eigenvalues');
    [eigenvalues ind] = sort(diag(Lambda));
    eigenvalues = flipud(eigenvalues);
    ind = flipud(ind);

    if (saveEig == 4),
        disp('Saving Eigen Decomposition 1');
        save eig_stuff1 Lambda V covarianceMatrix eigenvalues;
        clear covarianceMatrix;
    end

end

%estimate the rank (Too long for tica!)...
%maxLastEig = rank(covarianceMatrix);
maxLastEig = lastEig;

disp('Selecting Desired Eigenvalues');
% Remove the smaller eigenvalues
%disp('Removing smaller eigenvalues');
if lastEig < originalDimension
    lowerLimitValue = eigenvalues(lastEig);
else
    lowerLimitValue = eigenvalues(originalDimension);
end
lowerColumns = (diag(Lambda) >= lowerLimitValue);

% Remove the larger eigenvalues
higherLimitValue = eigenvalues(firstEig);
higherColumns = diag(Lambda) <= higherLimitValue;

% Combine the results from above
selectedColumns = lowerColumns & higherColumns;

V_orig = V(:, ind);

V = v_selcol(V,selectedColumns);

Lambda_orig = eigenvalues;
Lambda = v_selcol(v_selcol(Lambda, selectedColumns)',selectedColumns);

sumAll=sum(eigenvalues);
sumUsed=sum(diag(Lambda));
retained = (sumUsed/sumAll)*100;
fprintf('%g%% of (non-zero) eigenvalues retained.\n', retained);


function [X, orig_mean] = nic_icatb_remove_mean(X, verbose)
%% Remove mean of each column and return original mean values in a vector.
%
% Inputs:
% 1. X - X is 2D double array.
% 2. verbose - Options are 0 and 1.
%
% Outputs:
% 1. X - Data with zero mean
% 2. orig_mean - Original mean values
%

if (~exist('verbose', 'var'))
    verbose = 1;
end

orig_mean = mean(X);

try
    if (~isempty(which('bsxfun.m')))
        %% Use bsxfun if possible
        X = bsxfun(@minus, X, orig_mean);
    else
        tempVar = repmat(orig_mean, size(X, 1), 1);
        X = X - tempVar;
        clear tempVar;
    end
catch
    clear tempVar;
    if (verbose)
        disp('Using slow but less memory for removing mean in the data');
    end
    %% Remove mean of sample one at a time
    for nM = 1:size(X, 2)
        X(:, nM) = X(:, nM) - orig_mean(nM);
    end
end

function [newVectors, whiteningMatrix, dewhiteningMatrix] = nic_icatb_v_whiten(data, V, Lambda, doTranspose);
% [newVectors, whiteningMatrix, dewhiteningMatrix] = icatb_v_whiten(data, V, Lambda);


if exist('doTranspose', 'var')
    if strcmp(doTranspose, 'transpose')
        data = data';
    end
end

% Calculate the whitening and dewhitening matrices (these handle
% dimensionality simultaneously).
%whiteningMatrix = inv(sqrtm(Lambda))* V';
whiteningMatrix = sqrtm(Lambda) \ V'; % Use gaussian elimination approach to solve the equations
dewhiteningMatrix = V * sqrtm(Lambda);

% Project to the eigenvectors of the covariance matrix.
% Whiten the samples and reduce dimension simultaneously.
fprintf('Whitening...\n');

%newVectors =  whiteningMatrix*data;
newVectors = data * whiteningMatrix';
newVectors = newVectors';


function A = nic_icatb_cov(X, removeMean)
%% Compute covariance matrix. If out of memory error occurs, covariance matrix is
% computed with least memory requirements.
%
% Inputs:
% 1. X - Data of size m by n
% 2. removeMean - Remove mean. Options are 0 and 1.
%
% Outputs:
%
% 1. A - Covariance matrix
%

if (~exist('X', 'var'))
    error('Data variable is not passed');
end

%% Use defaults
if (~exist('removeMean', 'var'))
    removeMean = 1;
end

[m, n] = size(X);

%% Remove mean
if (removeMean)
    X = icatb_remove_mean(X, 1);
end

try
    A = X'*X;
catch
    %% Use less memory for computing covariance matrix
    A = zeros(n, n);

    %% Loop over cols
    for i = 1:n
        currentVector = X(:, i);
        inds = (i:n);
        % Loop over nInd
        for nInd = inds
            temp = currentVector(:)'*X(:, nInd);
            %% Utilise the symmetric nature of the covariance matrix
            A(i, nInd) = temp;
            A(nInd, i) = temp;
        end
        % End loop over nInd
    end
    %% End loop over cols

end

%% Return the final result
if (m > 1)
    A = A ./ (m - 1);
else
    A = A ./ m;
end



function newMatrix = v_selcol (oldMatrix, maskVector);
% newMatrix = selcol(oldMatrix, maskVector);
%
% Selects the columns of the matrix that marked by one in the given vector.
% The maskVector is a column vector.
%
% This function is needed by PCAMAT

% 15.3.1998

if size (maskVector, 1) ~= size (oldMatrix, 2),
  error ('The mask vector and matrix are of uncompatible size.');
end

numTaken = 0;

% take only the indices in the mask
takingMask = find(maskVector == 1);

numTaken = length(takingMask);

newMatrix = oldMatrix(:, takingMask);

