# This function computes the data for a lineage-through-time plot and (optionally) creates this plot
# The function does not require a tree that is ultrametric.  Optionally, the function can remove extinct
# species from the phylogeny
# Written by Liam J. Revell 2010

ltt<-function(tree,plot=TRUE,drop.extinct=FALSE,log.lineages=TRUE){

	# require
	require(ape)

	# global
	tol<-1e-6

	# internal functions
	# drop extinct tips
	drop.extinct.tips<-function(phy){
		tol<-1e-6
		temp<-diag(vcv(phy))
		pruned.phy<-drop.tip(phy,names(temp[temp<(max(temp)-tol)]))
		return(pruned.phy)
	}

	# first, if drop.extinct==TRUE
	if(drop.extinct==TRUE)
		tree<-drop.extinct.tips(tree)

	# compute node heights
	root<-length(tree$tip)+1
	node.height<-matrix(NA,nrow(tree$edge),2)
	for(i in 1:nrow(tree$edge)){
		if(tree$edge[i,1]==root){
			node.height[i,1]<-0.0
			node.height[i,2]<-tree$edge.length[i]
		} else {
			node.height[i,1]<-node.height[match(tree$edge[i,1],tree$edge[,2]),2]
			node.height[i,2]<-node.height[i,1]+tree$edge.length[i]
		}
	}

	tree.length<-max(node.height) # tree length
	ltt<-vector()
	time<-c(0,node.height[,2]); names(time)<-as.character(c(root,tree$edge[,2]))
	temp<-vector()
	time<-time[order(time)]
	i=1; while(time[i]<(tree.length-tol)) i=i+1 
	time<-time[1:i] # times

	# get numbers of lineages
	for(i in 1:(length(time)-1)){
		ltt[i]<-0
		for(j in 1:nrow(node.height))
			ltt[i]<-ltt[i]+(time[i]>=(node.height[j,1]-tol)&&time[i]<=(node.height[j,2]-tol))
	}
	ltt[i+1]<-0
	for(j in 1:nrow(node.height))
		ltt[i+1]<-ltt[i+1]+(time[i+1]<=(node.height[j,2]+tol))

	# set names (these are the node indices)
	names(ltt)<-names(time)

	# append 0,1 for time 0
	ltt<-c(1,ltt); time<-c(0,time)

	# plot ltt
	if(plot==TRUE){
		if(log.lineages==TRUE)
			plot(time,log(ltt),"s",xlab="time",ylab="log(lineages)")
		else
			plot(time,ltt,"s",xlab="time",ylab="lineages")
	}

	return(list(ltt=ltt,times=time))

}