# function to perform phylogenetic principal components analysis
# multiple morphological traits in Y
# also can use lambda transformation in which lambda is optimized by ML
# written by Liam Revell 2010/2011, ref. Revell (2009; Evolution)

phyl.pca<-function(tree,Y,method="BM",mode="cov"){
	# check 'ape'
	if(!require(ape)) stop("function requires 'ape' package. please install.")
	# check tree
	if(class(tree)!="phylo") stop("tree must be an object of class 'phylo.'")

	# check and sort data
	Y<-as.matrix(Y)
	if(is.null(rownames(Y))){
		print("Y has no names. function will assume that the row order of Y matches tree$tip.label")
		rownames(Y)<-tree$tip.label
	} else Y<-as.matrix(Y[tree$tip.label,]) # sort

	# preliminaries
	n<-length(tree$tip); m<-ncol(X)

	# analyze
	if(method=="BM"){ 
		temp<-phyl.vcv(X,tree,1.0)
		V<-temp$R; a<-t(temp$alpha); C<-temp$C
		rm(temp)
	} else if(method=="lambda"){
		temp<-optimize(f=likelihood,interval=c(0,1),X=X,phy=tree,maximum=TRUE)
		lambda<-temp$maximum; logL<-as.numeric(temp$objective)
		rm(temp)
		temp<-phyl.vcv(X,tree,lambda)
		V<-temp$R; a<-t(temp$alpha); C<-temp$C
		rm(temp)
	}
	invC<-solve(C)
	# if correlation matrix
	if(mode=="corr"){
		X=X/matrix(rep(sqrt(diag(V)),n),n,m,byrow=T) # standardize X
		V=V/(sqrt(diag(V))%*%t(sqrt(diag(V)))) # change V to correlation matrix
		a<-matrix(colSums(invC%*%X)/sum(invC),m,1) # recalculate a
	}
	es=eigen(V) # eigenanalyze
	result<-list(); result$Eval<-diag(es$values); result$Evec<-es$vectors
	A<-matrix(rep(a,n),n,m,byrow=T)
	result$S<-(X-A)%*%result$Evec # compute scores in the species space
	Ccv<-t(X-A)%*%invC%*%result$S/(n-1) # compute cross covariance matrix and loadings
	result$L<-matrix(,m,m)
	for(i in 1:m) for(j in 1:m) result$L[i,j]<-Ccv[i,j]/sqrt(V[i,i]*result$Eval[j,j])
	if(method=="lambda"){ 
		result$lambda<-lambda
		result$logL.lambda<-logL
	}
	
	# return result
	return(result)
}

# function to compute phylogenetic VCV using joint Pagel's lambda
# written by Liam Revell 2011
phyl.vcv<-function(X,phy,lambda){
	C<-vcv.phylo(phy)
	C<-lambda.transform(lambda,C)
	invC<-solve(C)
	a<-matrix(colSums(invC%*%X)/sum(invC),ncol(X),1)
	A<-matrix(rep(a,nrow(X)),nrow(X),ncol(X),byrow=T)
	V<-t(X-A)%*%invC%*%(X-A)/(nrow(C)-1)
	return(list(C=C,R=V,alpha=a))
}

# lambda transformation of C
# written by Liam Revell 2011
lambda.transform<-function(lambda,C){
	if(lambda==1) return(C)
	else {
		V<-diag(diag(C))
		C<-C-V
		C.lambda<-(V+lambda*C)
		return(C.lambda)
	}
}

# likelihood function
# written by Liam Revell 2011
likelihood<-function(theta,X,phy){
	C<-vcv(phy)
	C<-lambda.transform(theta,C)
	# compute R, conditioned on lambda
	temp<-phyl.vcv(X,phy,theta);
	R<-temp$R; a<-temp$alpha; rm(temp)
	# prep
	n<-nrow(X); m<-ncol(X); D<-matrix(0,n*m,m)
	for(i in 1:(n*m)) for(j in 1:m) if((j-1)*n<i&&i<=j*n) D[i,j]=1.0
	y<-as.matrix(as.vector(X))
	# compute the log-likelihood	
	logL<--t(y-D%*%a)%*%solve(kronecker(R,C))%*%(y-D%*%a)/2-n*m*log(2*pi)/2-log(det(kronecker(R,C)))/2
	return(logL)
}