% function [S,Eval,Evec,L]=phyl_pca(C,X) performs PCA using tree (C), data 
% (X), and mode ('cov' or 'corr')
function [S,Eval,Evec,L]=phyl_pca(C,X,mode)

    % check mode (covariance or correlation matrix)
    if exist('mode','var')~=1
        mode='cov';
    end
    
    % get m, n
    [n,m]=size(X);

    % first compute the vector of ancestral states
    one=ones(n,1); a=(one'*C^-1*one)^-1*(one'*C^-1*X)';
    
    % now compute the evolutionary VCV matrix, V
    V=(X-one*a')'*C^-1*(X-one*a')/(n-1);

    % if mode==correlation matrix
    if strcmp(mode,'corr')
        
        % standardize X
        X=X./(one*sqrt(diag(V))');
        
        % change V to correlation matrix
        V=V./(sqrt(diag(V))*sqrt(diag(V))');
        
        % recalculate a
        a=(one'*C^-1*one)^-1*(one'*C^-1*X)';
    
    end    
        
    % do eigenanalysis & sort
    [Evec,Eval]=eig(V); [Evec,Eval]=sorteig(Evec,Eval);

    % now compute scores in the original space
    S=(X-one*a')*Evec;
    
    % compute cross covariance matrix and loadings
    Ccv=(X-one*a')'*C^-1*S/(n-1); L=Ccv.*(diag(V)*diag(Eval)').^(-1/2);

% done