function [pcacomp,pcaraw] = dismatrix_pca(data,comp,comptype)
%% pca-based distance matrix construction
% Input
%   data      - subjid x voxels
%   comp      - (0-1) / 10. example 90% components / the first 10 components
%   comptype  - accumulative contribution rate or  number of components
% Output
%   pcacomp   - components
%   pcaraw    - components back to raw data space
disp('Do PCA. Waiting......');
data         = data - repmat(mean(data,1),size(data,1),1); % centering data
covmatrix    = data'*data;
[u,d,v]      = svd(covmatrix);
d            = diag(d);
rato         = cumsum(d)./sum(d);

if isequal(comptype,'number')
    pcacomp = data*u(:,1:ceil(comp));
    fprintf('Extract the first %d components\n',ceil(comp));
elseif isequal(comptype,'percent')
    for i=1:length(rato)
        if rato(i) >= comp
            comp = i;
            break;
        end
    end
    pcacomp = data*u(:,1:ceil(comp));
    fprintf('Extract the first %d components\n',ceil(comp));
    clear i
end

pcaraw = data*u(:,1:ceil(comp))*u(:,1:ceil(comp))';
disp('PCA finished');
end