R/utilities.R

Defines functions untangle matchType paste.tree multiC getDescendants getSisters orderMappedEdge maxLambda matchTreetoData matchDatatoTree likMlambda lambda.transform phyl.vcv applyBranchLengths matchNodes countSimmap getCladesofSize extract.clade.simmap findMRCA collapse.to.star bind.tip sampleFrom rotateNodes mergeMappedStates roundBranches drop.leaves nodeHeights di2multi.simmap lambdaTree vcvPhylo midpoint.root midpoint_root getAncestors fastHeight fastMRCA nodeheight plot.describe.simmap print.describe.simmap drop.tip.singleton rep.multiPhylo rep.phylo repPhylo ladderize.simmap reroot drop.clade splitTree getExtinct getExtant labelSubTree cladelabels markChanges rstate whichorder reorder.simmap reorderSimmap fastDist get.treepos edgeProbs matchLabels mapped.states as.multiPhylo as.multiPhylo.phylo as.multiPhylo.multiSimmap as.phylo.simmap bind.tree.simmap rotate.multi allRotations modified.Grafen node.paths print.aic.w aic.w bd labelnodes getnode arc.cladelabels linklabels cladebox get.asp print.geo.palette geo.palette print.geo.legend geo.legend plot.expand.clade print.expand.clade expand.clade ebTree multi2di.multiSimmap di2multi.multiSimmap multi2di.simmap multi2di.densityMap multi2di.contMap di2multi.densityMap di2multi.contMap getStates rescaleSimmap rescale.multiSimmap rescale.simmap rescale.default rescale Nedge.densityMap Nnode.densityMap Ntip.densityMap Nedge.contMap Nnode.contMap Ntip.contMap add.arrow extract.clade.simmap describe.simmap c.multiSimmap c.simmap force.ultrametric

Documented in add.arrow aic.w allRotations applyBranchLengths arc.cladelabels as.multiPhylo as.multiPhylo.multiSimmap as.multiPhylo.phylo bd bind.tip bind.tree.simmap cladelabels collapse.to.star countSimmap describe.simmap di2multi.contMap di2multi.densityMap di2multi.multiSimmap di2multi.simmap drop.clade drop.leaves drop.tip.singleton edgeProbs expand.clade extract.clade.simmap fastDist fastHeight fastMRCA findMRCA force.ultrametric geo.legend geo.palette getCladesofSize getDescendants getExtant getExtinct getnode getSisters getStates get.treepos labelnodes ladderize.simmap lambda.transform likMlambda linklabels mapped.states markChanges matchLabels matchNodes mergeMappedStates midpoint_root midpoint.root modified.Grafen multi2di.contMap multi2di.densityMap multi2di.multiSimmap multi2di.simmap multiC nodeheight nodeHeights node.paths orderMappedEdge paste.tree phyl.vcv plot.describe.simmap plot.expand.clade reorderSimmap rep.multiPhylo rep.phylo repPhylo reroot rescale rescale.multiSimmap rescale.simmap rescaleSimmap rotate.multi rotateNodes roundBranches rstate sampleFrom splitTree untangle vcvPhylo

## some utility functions
## written by Liam J. Revell 2011, 2012, 2013, 2014, 2015, 2016, 2017, 
## 2018, 2019, 2020, 2021, 2022, 2023, 2024

## function forces a tree to be ultrametric using two different methods
## written by Liam J. Revell 2017, 2021, 2022, 2023

force.ultrametric<-function(tree,method=c("nnls","extend"),...){
	if(hasArg(message)) message<-list(...)$message
	else message<-TRUE
	if(message){
		cat("***************************************************************\n")
		cat("*                          Note:                              *\n")
		cat("*    force.ultrametric does not include a formal method to    *\n")
		cat("*    ultrametricize a tree & should only be used to coerce    *\n")
		cat("*   a phylogeny that fails is.ultrametric due to rounding --  *\n")
		cat("*    not as a substitute for formal rate-smoothing methods.   *\n")
		cat("***************************************************************\n")
	}
	method<-method[1]
	if(method=="nnls") tree<-nnls.tree(cophenetic(tree),tree,
		method="ultrametric",rooted=is.rooted(tree),trace=0)
	else if(method=="extend"){
		h<-diag(vcv(tree))
		d<-max(h)-h
		ii<-sapply(1:Ntip(tree),function(x,y) which(y==x),
			y=tree$edge[,2])
		tree$edge.length[ii]<-tree$edge.length[ii]+d
	} else 
		cat("method not recognized: returning input tree\n\n")
	tree
}

## c "combine" method for "simmap" and "multiSimmap" object classes

## adapted from c.phylo in ape (Paradis & Schliep 2019)
c.simmap<-function(...,recursive=TRUE){
	obj<-list(...)
	classes<-lapply(obj,class)
	isphylo<-sapply(classes,function(x) "simmap" %in% x)
	if(all(isphylo)){
		class(obj)<-c("multiSimmap","multiPhylo")
		return(obj)
	}
	if(!recursive) return(obj)
	ismulti<-sapply(classes, function(x) "multiSimmap" %in% x)
	if(all(isphylo|ismulti)){
		result<-list()
		j<-1
		for(i in 1:length(isphylo)){
			if(isphylo[i]){ 
				result[[j]]<-obj[[i]]
				j<-j+1
			} else {
				n<-length(obj[[i]])
				result[0:(n-1)+j]<-.uncompressTipLabel(obj[[i]])
				j<-j+n
			}
		}
		class(result)<-c("multiSimmap","multiPhylo")
		obj<-result
	} else {
		msg<-paste("some objects not of class \"simmap\" or",
			"\"multiSimmap\": argument recursive=TRUE ignored")
		warning(msg)
	}
	return(obj)
}

## adapted from c.multiPhylo in ape (Paradis & Schliep 2019)
c.multiSimmap<-function(...,recursive=TRUE){
	obj<-list(...)
	if(!recursive) return(obj)
	classes<-lapply(obj,class)
	isphylo<-sapply(classes,function(x) "simmap" %in% x)
	ismulti<-sapply(classes, function(x) "multiSimmap" %in% x)
	if(all(isphylo|ismulti)){
		result<-list()
		j<-1
		for(i in 1:length(isphylo)){
			if(isphylo[i]){ 
				result[[j]]<-obj[[i]]
				j<-j+1
			} else {
				n<-length(obj[[i]])
				result[0:(n-1)+j]<-.uncompressTipLabel(obj[[i]])
				j<-j+n
			}
		}
		class(result)<-c("multiSimmap","multiPhylo")
		obj<-result
	} else {
		msg<-paste("some objects not of class \"simmap\" or",
			"\"multiSimmap\": argument recursive=TRUE ignored")
		warning(msg)
	}
	return(obj)
}

## function to summarize the results of stochastic mapping
## written by Liam J. Revell 2013, 2014, 2015, 2021, 2022, 2023

describe.simmap<-function(tree,...){
	if(hasArg(plot)) plot<-list(...)$plot
	else plot<-FALSE
	if(hasArg(check.equal)) check.equal<-list(...)$check.equal
	else check.equal<-FALSE
	if(hasArg(message)) message<-list(...)$message
	else message<-FALSE
	if(hasArg(ref.tree)) ref.tree<-list(...)$ref.tree
	else ref.tree<-NULL
	if(hasArg(states)) states<-list(...)$states
	else states<-NULL
	if(inherits(tree,"multiPhylo")){
		if(check.equal){
			TT<-sapply(tree,function(x,y) sapply(y,all.equal.phylo,x),y=tree)
			check<-all(TT)
			if(!check) cat("Note: Some trees are not equal.\nA \"reference\" tree will be computed if none was provided.\n\n")
		} else check<-TRUE
		if(is.null(ref.tree)&&check){ 
			YY<-getStates(tree,"both")
			if(is.null(states)) states<-sort(unique(as.vector(YY)))
			ZZ<-t(apply(YY,1,function(x,levels,Nsim) summary(factor(x,levels))/Nsim,
				levels=states,Nsim=length(tree)))
		} else {
			YY<-getStates(tree)
			if(is.null(states)) states<-sort(unique(as.vector(YY)))
			if(is.null(ref.tree)){
				cat("No reference tree provided & some trees are unequal.\nComputing majority-rule consensus tree.\n")
				ref.tree<-consensus(tree,p=0.5)
			}
			YYp<-matrix(NA,ref.tree$Nnode,length(tree),dimnames=list(1:ref.tree$Nnode+Ntip(ref.tree),
				NULL))
			for(i in 1:length(tree)){
				M<-matchNodes(ref.tree,tree[[i]])
				jj<-sapply(M[,2],function(x,y) if(x%in%y) which(as.numeric(y)==x) else NA,
					y=as.numeric(rownames(YY)))
				YYp[,i]<-YY[jj,i]
			}
			ZZ<-t(apply(YYp,1,function(x,levels) summary(factor(x[!is.na(x)],
				levels))/sum(!is.na(x)),levels=states))
		}
		XX<-countSimmap(tree,states,FALSE)
		XX<-XX[,-(which(as.vector(diag(-1,length(states)))==-1)+1)]
		AA<-lapply(unclass(tree),function(x) setNames(
			c(colSums(x$mapped.edge),sum(x$edge.length)),
			c(colnames(x$mapped.edge),"total")))
		foo<-function(x,n){
			y<-setNames(rep(0,length(n)),n)
			y[names(x)]<-x
			y
		}
		AA<-t(sapply(AA,foo,n=c(states,"total")))
		BB<-getStates(tree,type="tips")
		CC<-t(apply(BB,1,function(x,levels,Nsim) summary(factor(x,levels))/Nsim,levels=states,
			Nsim=length(tree)))
		x<-list(count=XX,times=AA,ace=ZZ,tips=CC,tree=tree,ref.tree=if(!is.null(ref.tree)) 
			ref.tree else NULL)
		class(x)<-"describe.simmap"
	} else if(inherits(tree,"phylo")){
		XX<-countSimmap(tree,message=FALSE)
		YY<-getStates(tree)
		if(is.null(states)) states<-sort(unique(YY))
		AA<-setNames(c(colSums(tree$mapped.edge),sum(tree$edge.length)),
			c(colnames(tree$mapped.edge),"total"))
		AA<-rbind(AA,AA/AA[length(AA)]); rownames(AA)<-c("raw","prop")
		x<-list(N=XX$N,Tr=XX$Tr,times=AA,states=YY,tree=tree)
		class(x)<-"describe.simmap"
	}
	if(message) print(x)
	if(plot) plot(x)
	x
}

# function works like extract.clade in ape but will preserve a discrete character mapping
# written by Liam J. Revell 2013, 2023
extract.clade.simmap<-function(tree,node){
	if(!inherits(tree,"simmap")) stop("tree should be an object of class \"simmap\".")
	x<-getDescendants(tree,node)
	x<-x[x<=Ntip(tree)]
	drop.tip.simmap(tree,tree$tip.label[-x])
}

## function to add an arrow pointing to a tip or node in the tree
## written by Liam J. Revell 2014, 2017, 2020, 2023

add.arrow<-function(tree=NULL,tip,...){
	if(length(tip)>1){ 
		object<-sapply(tip,add.arrow,tree=tree,...)
		invisible(object)
	} else {
		lastPP<-get("last_plot.phylo",envir=.PlotPhyloEnv)
		asp<-if(lastPP$type=="fan") 1 else (par()$usr[4]-par()$usr[3])/(par()$usr[2]-
			par()$usr[1])
		if(!is.null(tree)){
			if(inherits(tree,"contMap")) tree<-tree$tree
			else if(inherits(tree,"densityMap")) tree<-tree$tree
		}
		if(is.numeric(tip)){
			ii<-tip
			if(!is.null(tree)&&ii<=Ntip(tree)) tip<-tree$tip.label[ii]
			else tip<-""
		} else if(is.character(tip)&&!is.null(tree)) ii<-which(tree$tip.label==tip)
		if(hasArg(offset)) offset<-list(...)$offset
		else offset<-lastPP$label.offset
		strw<-lastPP$cex*(strwidth(tip)+offset*mean(strwidth(c(LETTERS,letters))))
		if(lastPP$direction%in%c("upwards","downwards")) 
			strw<-strw*asp*par()$pin[1]/par()$pin[2]
		if(hasArg(arrl)) arrl<-list(...)$arrl
		else { 
			if(lastPP$type=="fan") arrl<-0.3*max(lastPP$xx)
			else if(lastPP$type=="phylogram") arrl<-0.15*max(lastPP$xx)
		}
		if(hasArg(hedl)) hedl<-list(...)$hedl
		else hedl<-arrl/3
		if(hasArg(angle)) angle<-list(...)$angle
		else angle<-45
		arra<-angle*pi/180
		if(hasArg(col)) col<-list(...)$col
		else col<-"black"
		if(hasArg(lwd)) lwd<-list(...)$lwd	
		else lwd<-2
		if(lastPP$type=="fan") theta<-atan2(lastPP$yy[ii],lastPP$xx[ii])
		else if(lastPP$type=="phylogram"){
			if(lastPP$direction=="rightwards") theta<-0
			else if(lastPP$direction=="upwards") theta<-pi/2 
			else if(lastPP$direction=="leftwards") theta<-pi
			else if(lastPP$direction=="downwards") theta<-3*pi/2
		}
		segments(x0=lastPP$xx[ii]+cos(theta)*(strw+arrl),
			y0=lastPP$yy[ii]+sin(theta)*(strw+arrl),
			x1=lastPP$xx[ii]+cos(theta)*strw,
			y1=lastPP$yy[ii]+sin(theta)*strw,
			col=col,lwd=lwd,lend="round")
		segments(x0=lastPP$xx[ii]+cos(theta)*strw+cos(theta+arra/2)*hedl,
			y0=lastPP$yy[ii]+sin(theta)*strw+sin(theta+arra/2)*hedl*asp,
			x1=lastPP$xx[ii]+cos(theta)*strw,
			y1=lastPP$yy[ii]+sin(theta)*strw,
			col=col,lwd=lwd,lend="round")
		segments(x0=lastPP$xx[ii]+cos(theta)*strw+cos(theta-arra/2)*hedl,
			y0=lastPP$yy[ii]+sin(theta)*strw+sin(theta-arra/2)*hedl*asp,
			x1=lastPP$xx[ii]+cos(theta)*strw,
			y1=lastPP$yy[ii]+sin(theta)*strw,
			col=col,lwd=lwd,lend="round")
		invisible(list(x0=lastPP$xx[ii]+cos(theta)*(strw+arrl),
			y0=lastPP$yy[ii]+sin(theta)*(strw+arrl),
			x1=lastPP$xx[ii]+cos(theta)*strw,
			y1=lastPP$yy[ii]+sin(theta)*strw))
	}
}

## function to rescale simmap style trees
## written by Liam J. Revell 2012, 2013, 2014, 2015, 2017, 2023

## S3 Ntip etc. methods for "contMap" object class
Ntip.contMap<-function(phy) Ntip(phy$tree)
Nnode.contMap<-function(phy,...) Nnode(phy$tree)
Nedge.contMap<-function(phy) Nedge(phy$tree)

## S3 Ntip etc. methods for "densityMap" object class
Ntip.densityMap<-function(phy) Ntip(phy$tree)
Nnode.densityMap<-function(phy,...) Nnode(phy$tree)
Nedge.densityMap<-function(phy) Nedge(phy$tree)

rescale<-function(x,...) UseMethod("rescale")

rescale.default<-function(x,...){
	warning(paste(
		"rescale does not know how to handle objects of class ",class(x),".\n",
		"if",class(x),"= \"phylo\" load geiger package to rescale.\n"))
}

rescale.simmap<-function(x, model="depth", ...) rescaleSimmap(x,...)

rescale.multiSimmap<-function(x, model="depth", ...) rescaleSimmap(x,...)

rescaleSimmap<-function(tree,...){
	if(inherits(tree,"multiSimmap")){
		cls<-class(tree)
		trees<-unclass(tree)
		trees<-lapply(trees,rescaleSimmap,...)
		class(trees)<-cls
		return(trees)
	} else if(inherits(tree,"simmap")){
		if(hasArg(lambda)) lambda<-list(...)$lambda
		else lambda<-1
		if(hasArg(totalDepth)) depth<-totalDepth<-list(...)$totalDepth
		else if(hasArg(depth)) depth<-totalDepth<-list(...)$depth
		else depth<-totalDepth<-max(nodeHeights(tree))
		if(lambda!=1){
			e<-lambdaTree(tree,lambda)$edge.length/tree$edge.length
			tree$edge.length<-tree$edge.length*e
			tree$maps<-mapply(function(x,y) x*y,tree$maps,e)
			tree$mapped.edge<-tree$mapped.edge*matrix(rep(e,ncol(tree$mapped.edge)),length(e),ncol(tree$mapped.edge))
		}
		if(depth!=max(nodeHeights(tree))){
			h<-max(nodeHeights(tree))
			s<-depth/h
			tree$edge.length<-tree$edge.length*s
			tree$maps<-lapply(tree$maps,"*",s)
			tree$mapped.edge<-tree$mapped.edge*s
		}
		return(tree)
	} else message("tree should be an object of class \"simmap\" or \"multiSimmap\"")
}

## function to get states at internal nodes from simmap style trees
## written by Liam J. Revell 2013, 2014, 2015, 2021
getStates<-function(tree,type=c("nodes","tips","both")){
	type<-type[1]
	if(inherits(tree,"multiPhylo")){
		tree<-unclass(tree)
		obj<-lapply(tree,getStates,type=type)
		nn<-names(obj[[1]])
		y<-sapply(obj,function(x,n) x[n],n=nn)
	} else if(inherits(tree,"phylo")){ 
		if(type%in%c("nodes","both")){
			a<-setNames(sapply(tree$maps,function(x) names(x)[1]),tree$edge[,1])
			a<-a[as.character(Ntip(tree)+1:tree$Nnode)]
		}
		if(type%in%c("tips","both")){
			b<-setNames(sapply(tree$maps,function(x) names(x)[length(x)]),tree$edge[,2])
			b<-setNames(b[as.character(1:Ntip(tree))],tree$tip.label)
		}
		y<-if(type=="both") c(a,b) else if(type=="nodes") a else b
	} else stop("tree should be an object of class \"phylo\" or \"multiPhylo\".")
	return(y)
}

## di2multi & multi2di for "contMap" & "densityMap" object classes

di2multi.contMap<-function(phy,...){
	phy$tree<-di2multi(phy$tree,...)
	phy
}

di2multi.densityMap<-function(phy,...){
	phy$tree<-di2multi(phy$tree,...)
	phy
}

multi2di.contMap<-function(phy,...){
	phy$tree<-multi2di(phy$tree,...)
	phy
}

multi2di.densityMap<-function(phy,...){
	phy$tree<-multi2di(phy$tree,...)
	phy
}

## multi2di for "simmap" object class

multi2di.simmap<-function(phy,...){
	obj<-multi2di(as.phylo(phy),...)
	M<-rbind(matchNodes(obj,phy),
		matchLabels(obj,phy))
	obj$maps<-vector(mode="list",length=nrow(obj$edge))
	for(i in 2:nrow(M)){
		if(!is.na(M[i,2])){
			obj$maps[[which(obj$edge[,2]==M[i,1])]]<-
				phy$maps[[which(phy$edge[,2]==M[i,2])]]
		} else {
			ii<-which(obj$edge[,2]==getParent(obj,M[i,1]))
			state<-names(obj$maps[[ii]])[length(obj$maps[[ii]])]
			obj$maps[[which(obj$edge[,2]==M[i,1])]]<-
				setNames(0,state)
		}
	}
	obj$node.states<-getStates(obj,"nodes")
	obj$states<-getStates(obj,"tips")
	obj$mapped.edge<-makeMappedEdge(obj$edge,obj$maps)
	class(obj)<-c("simmap",class(obj))
	obj
}

## di2multi & multi2di for "multiSimmap" object class

di2multi.multiSimmap<-function(phy,...){
	obj<-lapply(phy,di2multi,...)
	class(obj)<-c("multiSimmap","multiPhylo")
	obj
}
multi2di.multiSimmap<-function(phy,...){
	obj<-lapply(phy,multi2di,...)
	class(obj)<-c("multiSimmap","multiPhylo")
	obj
}

## function to rescale a tree according to an EB model
## written by Liam J. Revell 2017

ebTree<-function(tree,r){
	if(r!=0){
		H<-nodeHeights(tree)
		e<-(exp(r*H[,2])-exp(r*H[,1]))/r
		tree$edge.length<-e
	}
	tree
}

## function to expand clades in a plot by a given factor
## written by Liam J. Revell 2017
expand.clade<-function(tree,node,factor=5){
	cw<-reorder(tree)
	tips<-setNames(rep(1,Ntip(tree)),cw$tip.label)
	get.tips<-function(node,tree){
				dd<-getDescendants(tree,node)
				tree$tip.label[dd[dd<=Ntip(tree)]]
	}
	desc<-unlist(lapply(node,get.tips,tree=cw))
	for(i in 2:Ntip(cw)){
		tips[i]<-tips[i-1]+
			if(names(tips)[i]%in%desc){
				1 
			} else if(names(tips)[i-1]%in%desc){
				1
			} else 1/factor
	}
	obj<-list(tree=tree,tips=tips)
	class(obj)<-"expand.clade"
	obj
}

## S3 print method for the object class "expand.clade"
print.expand.clade<-function(x,...){
	cat("An object of class \"expand.clade\" consisting of:\n")
	cat(paste("(1) A phylogenetic tree (x$tree) with",Ntip(x$tree),
		"tips and\n	 ",x$tree$Nnode,"internal nodes.\n"))
	cat("(2) A vector (x$tips) containing the desired tip-spacing.\n\n")
}

## S3 plot method for the object class "expand.clade"
plot.expand.clade<-function(x,...){
	args<-list(...)
	args$tree<-x$tree
	args$tips<-x$tips
	if(inherits(args$tree,"simmap")) do.call(plotSimmap,args)
	else do.call(plotTree,args)
}

## function to add a geological or other temporal legend to a plotted tree
## written by Liam J. Revell 2017, 2019
geo.legend<-function(leg=NULL,colors=NULL,alpha=0.2,...){
	if(hasArg(cex)) cex<-list(...)$cex
	else cex<-par()$cex
	if(hasArg(plot)) plot<-list(...)$plot
	else plot<-TRUE
	if(hasArg(show.lines)) show.lines<-list(...)$show.lines
	else show.lines<-TRUE
	obj<-get("last_plot.phylo",envir=.PlotPhyloEnv)
	if(is.null(colors)){
		colors<-setNames(c(
			rgb(255,242,127,255,maxColorValue=255),
			rgb(255,230,25,255,maxColorValue=255),
			rgb(253,154,82,255,maxColorValue=255),
			rgb(127,198,78,255,maxColorValue=255),
			rgb(52,178,201,255,maxColorValue=255),
			rgb(129,43,146,255,maxColorValue=255),
			rgb(240,64,40,255,maxColorValue=255),
			rgb(103,165,153,255,maxColorValue=255),
			rgb(203,140,55,255,maxColorValue=255),
			rgb(179,225,182,255,maxColorValue=255),
			rgb(0,146,112,255,maxColorValue=255),
			rgb(127,160,86,255,maxColorValue=255),
			rgb(247,67,112,255,maxColorValue=255)),
			c("Quaternary","Neogene","Paleogene",
			"Cretaceous","Jurassic","Triassic",
			"Permian","Carboniferous","Devonian",
			"Silurian","Ordovician","Cambrian",
			"Precambrian"))
	}
	if(is.null(leg)){
		leg<-rbind(c(2.588,0),
			c(23.03,2.588),
			c(66.0,23.03),
			c(145.0,66.0),
			c(201.3,145.0),
			c(252.17,201.3),
			c(298.9,252.17),
			c(358.9,298.9),
			c(419.2,358.9),
			c(443.8,419.2),
			c(485.4,443.8),
			c(541.0,485.4),
			c(4600,541.0))
		rownames(leg)<-c("Quaternary","Neogene","Paleogene",
			"Cretaceous","Jurassic","Triassic",
			"Permian","Carboniferous","Devonian",
			"Silurian","Ordovician","Cambrian",
			"Precambrian")
		t.max<-max(obj$xx)
		ii<-which(leg[,2]<=t.max)
		leg<-leg[ii,]
		leg[max(ii),1]<-t.max
	}
	colors<-sapply(colors,make.transparent,alpha=alpha)
	if(plot){	
		y<-c(rep(0,2),rep(par()$usr[4],2))
		ylabel<--1/25*obj$Ntip
		if(obj$direction=="rightwards"){
			old.usr<-par()$usr
			h<-max(obj$xx)
			new.xlim<-c(h-par()$usr[1],h-par()$usr[2])
			par(usr=c(new.xlim,old.usr[3:4]))
		} else old.usr<-par()$usr
		for(i in 1:nrow(leg)){
			strh<-strheight(rownames(leg)[i])
			polygon(c(leg[i,1:2],leg[i,2:1]),y,
				col=colors[rownames(leg)[i]],border=NA)
			if(show.lines){
				lines(x=rep(leg[i,1],2),y=c(0,par()$usr[4]),
					lty="dotted",col="grey")
				lines(x=c(leg[i,1],mean(leg[i,])-0.8*cex*
					get.asp()*strheight(rownames(leg)[i])),
					y=c(0,ylabel),lty="dotted",col="grey")
				lines(x=c(leg[i,2],mean(leg[i,])+0.8*cex*
					get.asp()*strheight(rownames(leg)[i])),
					y=c(0,ylabel),lty="dotted",col="grey")
				lines(x=rep(mean(leg[i,])-0.8*cex*
					get.asp()*strheight(rownames(leg)[i]),2),
					y=c(ylabel,par()$usr[3]),lty="dotted",col="grey")
				lines(x=rep(mean(leg[i,])+0.8*cex*
					get.asp()*strheight(rownames(leg)[i]),2),
					y=c(ylabel,par()$usr[3]),lty="dotted",col="grey")
			}
			polygon(x=c(leg[i,1],
				mean(leg[i,])-0.8*cex*get.asp()*strh,
				mean(leg[i,])-0.8*cex*get.asp()*strh,
				mean(leg[i,])+0.8*cex*get.asp()*strh,
				mean(leg[i,])+0.8*cex*get.asp()*strh,
				leg[i,2]),y=c(0,ylabel,par()$usr[3],
				par()$usr[3],ylabel,0),
				col=colors[rownames(leg)[i]],border=NA)
			strh<-strh*get.asp()
			text(x=mean(leg[i,])+
				if(obj$direction=="leftwards") 0.12*strh else -0.12*strh,
				y=ylabel,labels=rownames(leg)[i],
				srt=90,adj=c(1,0.5),cex=cex)
		}
	}
	par(usr=old.usr)
	object<-list(leg=leg,colors=colors[1:nrow(leg)])
	class(object)<-"geo.legend"
	invisible(object)
}

print.geo.legend<-function(x,...){
	cat("A geological period legend:\n")
	colnames(x$leg)<-c("start","end")
	print(data.frame(x$leg,color=x$colors))
	cat("\n")
}

geo.palette<-function(){
	colors<-setNames(c(
		rgb(255,242,127,255,maxColorValue=255),
		rgb(255,230,25,255,maxColorValue=255),
		rgb(253,154,82,255,maxColorValue=255),
		rgb(127,198,78,255,maxColorValue=255),
		rgb(52,178,201,255,maxColorValue=255),
		rgb(129,43,146,255,maxColorValue=255),
		rgb(240,64,40,255,maxColorValue=255),
		rgb(103,165,153,255,maxColorValue=255),
		rgb(203,140,55,255,maxColorValue=255),
		rgb(179,225,182,255,maxColorValue=255),
		rgb(0,146,112,255,maxColorValue=255),
		rgb(127,160,86,255,maxColorValue=255),
		rgb(247,67,112,255,maxColorValue=255)),
		c("Quaternary","Neogene","Paleogene",
		"Cretaceous","Jurassic","Triassic",
		"Permian","Carboniferous","Devonian",
		"Silurian","Ordovician","Cambrian",
		"Precambrian"))
	leg<-rbind(c(2.588,0),
		c(23.03,2.588),
		c(66.0,23.03),
		c(145.0,66.0),
		c(201.3,145.0),
		c(252.17,201.3),
		c(298.9,252.17),
		c(358.9,298.9),
		c(419.2,358.9),
		c(443.8,419.2),
		c(485.4,443.8),
		c(541.0,485.4),
		c(4600,541.0))
	rownames(leg)<-c("Quaternary","Neogene","Paleogene",
		"Cretaceous","Jurassic","Triassic",
		"Permian","Carboniferous","Devonian",
		"Silurian","Ordovician","Cambrian",
		"Precambrian")
	colnames(leg)<-c("start","end")
	object<-list(period=leg,cols=colors)
	class(object)<-"geo.palette"
	object
}

print.geo.palette<-function(x,...){
	cat("A geological period color palette:\n")
	print(data.frame(x$period,color=x$cols))
	cat("\n")
}

## borrowed from mapplots
get.asp<-function(){
	pin<-par("pin")
	usr<-par("usr")
	asp<-(pin[2]/(usr[4]-usr[3]))/(pin[1]/(usr[2]-usr[1]))
	asp
}

# round.polygon<-function(x,y,col="transparent"){
	## just space holding for now	
# }

## draw a box around a clade
## written by Liam J. Revell 2017

cladebox<-function(tree,node,color=NULL,...){
	if(is.null(color)) color<-make.transparent("yellow",0.2)
	obj<-get("last_plot.phylo",envir=.PlotPhyloEnv)
	h<-max(nodeHeights(tree))
	parent<-tree$edge[which(tree$edge[,2]==node),1]
	x0<-max(c(obj$xx[node]+obj$xx[parent])/2,obj$xx[node]-0.05*h)
	x1<-obj$x.lim[2]
	dd<-getDescendants(tree,node)
	y0<-min(range(obj$yy[dd]))-0.5
	y1<-max(range(obj$yy[dd]))+0.5
	polygon(c(x0,x1,x1,x0),c(y0,y0,y1,y1),col=color,
		border=0)
}

## draw tip labels as linking lines to text
## written by Liam J. Revell 2017

linklabels<-function(text,tips,link.type=c("bent","curved","straight"),
	...){
	lastPP<-get("last_plot.phylo",envir=.PlotPhyloEnv)
	if(!(lastPP$direction%in%c("leftwards","rightwards")))
		stop("direction should be \"rightwards\" or \"leftwards\".")
	if(hasArg(cex)) cex<-list(...)$cex
	else cex<-1
	if(hasArg(col)) col<-list(...)$col
	else col<-"black"
	if(hasArg(lty)) lty<-list(...)$lty
	else lty<-"dashed"
	if(hasArg(lwd)) lwd<-list(...)$lwd
	else lwd<-1
	if(hasArg(link.offset)) link.offset<-list(...)$link.offset
	else link.offset<-0.1*max(lastPP$xx)
	if(hasArg(font)) font<-list(...)$font
	else font<-3
	link.type<-link.type[1]
	xpos<-lastPP$xx[tips]+strwidth("i")
	ypos<-lastPP$yy[tips]
	xmax<-rep(max(lastPP$xx),length(tips))+link.offset
	ylab<-seq(min(lastPP$yy),max(lastPP$yy),
		by=(max(lastPP$yy)-min(lastPP$yy))/(length(tips)-1))
	ylab<-ylab[rank(ypos)]
	text(xmax,ylab,gsub("_"," ",text),pos=4,font=font,cex=cex,
		offset=0)
	if(link.type=="curved"){
		for(i in 1:length(tips))
			drawCurve(c(xpos[i],xmax[i]),c(ypos[i],ylab[i]),
				scale=0.01*diff(range(lastPP$xx)),lty=lty,
				col=col,lwd=lwd)
	} else if(link.type=="bent"){
		tipmax<-max(lastPP$xx)
		for(i in 1:length(tips)){
			ff<-strwidth("W")
			segments(xpos[i],ypos[i],tipmax+link.offset/2,ypos[i],
				lty=lty,col=col,lwd=lwd)
			segments(tipmax+link.offset/2,ypos[i],tipmax+
				link.offset/2+ff,ylab[i],lty=lty,col=col,lwd=lwd)
			segments(tipmax+link.offset/2+ff,ylab[i],xmax[i],ylab[i],
				lty=lty,col=col,lwd=lwd)
		}
	} else if(link.type=="straight")
		segments(xpos,ypos,xmax,ylab,lty=lty,col=col)
}

## function to create curved clade labels for a fan tree
## written by Liam J. Revell 2017, 2022

arc.cladelabels<-function(tree=NULL,text,node=NULL,ln.offset=1.02,
	lab.offset=1.06,cex=1,orientation="curved",stretch=1,...){
	obj<-get("last_plot.phylo",envir=.PlotPhyloEnv)
	if(obj$type!="fan") stop("method works only for type=\"fan\"")
	h<-max(sqrt(obj$xx^2+obj$yy^2))
	if(hasArg(mark.node)) mark.node<-list(...)$mark.node
	else mark.node<-TRUE
	if(hasArg(interactive)) interactive<-list(...)$interactive
	else {
		if(is.null(node)) interactive<-TRUE
		else interactive<-FALSE
	}
	if(interactive) node<-getnode()
	if(hasArg(lwd)) lwd<-list(...)$lwd
	else lwd<-par()$lwd
	if(hasArg(col)) col<-list(...)$col
	else col<-par()$col
	if(hasArg(lend)) lend<-list(...)$lend
	else lend<-par()$lend
	if(hasArg(clockwise)) clockwise<-list(...)$clockwise
	else clockwise<-TRUE
	if(hasArg(n)) n<-list(...)$n
	else n<-0.05
	if(mark.node) points(obj$xx[node],obj$yy[node],pch=21,
		bg="red")
	if(is.null(tree)){
		tree<-list(edge=obj$edge,tip.label=1:obj$Ntip,
			Nnode=obj$Nnode)
		class(tree)<-"phylo"
	}
	d<-getDescendants(tree,node)
	d<-sort(d[d<=Ntip(tree)])
	deg<-atan(obj$yy[d]/obj$xx[d])*180/pi
	ii<-intersect(which(obj$yy[d]>=0),which(obj$xx[d]<0))
	deg[ii]<-180+deg[ii]
	ii<-intersect(which(obj$yy[d]<0),which(obj$xx[d]<0))
	deg[ii]<-180+deg[ii]
	ii<-intersect(which(obj$yy[d]<0),which(obj$xx[d]>=0))
	deg[ii]<-360+deg[ii]
	draw.arc(x=0,y=0,radius=ln.offset*h,deg1=min(deg),
		deg2=max(deg),lwd=lwd,col=col,lend=lend,n=n)
	if(orientation=="curved")
		arctext(text,radius=lab.offset*h,
			middle=mean(range(deg*pi/180)),cex=cex,
			clockwise=clockwise,stretch=stretch)
	else if(orientation=="horizontal"){
		x0<-lab.offset*cos(median(deg)*pi/180)*h
		y0<-lab.offset*sin(median(deg)*pi/180)*h
		text(x=x0,y=y0,label=text,
		adj=c(if(x0>=0) 0 else 1,if(y0>=0) 0 else 1),
		offset=0,cex=cex)
	}
}

## function to return a node index interactively from a plotted tree
## written by Liam J. Revell 2017
getnode<-function(...){
	if(hasArg(env)) env<-list(...)$env
	else env<-get("last_plot.phylo",envir=.PlotPhyloEnv)
	if(hasArg(show.pt)) show.pt<-list(...)$show.pt
	else show.pt<-FALSE
	xy<-unlist(locator(n=1))
	if(show.pt) points(xy[1],xy[2])
	d<-sqrt((xy[1]-env$xx)^2+(xy[2]-env$yy)^2)
	ii<-which(d==min(d))[1]
	ii
}

## function mostly to interactively label nodes by clicking
## written by Liam J. Revell 2017, 2020
labelnodes<-function(text,node=NULL,interactive=TRUE,
	shape=c("circle","ellipse","rect"),...){
	shape<-shape[1]
	if(hasArg(circle.exp)) circle.exp<-list(...)$circle.exp
	else circle.exp<-1.3
	if(hasArg(rect.exp)) rect.exp<-list(...)$rect.exp
	else rect.exp<-1.6
	if(hasArg(cex)) cex<-list(...)$cex
	else cex<-1
	if(hasArg(bg)) bg<-list(...)$bg
	else bg<-"white"
	obj<-get("last_plot.phylo",envir=.PlotPhyloEnv)
	h<-cex*strheight("A")
	w<-cex*strwidth(text)
	rad<-circle.exp*h*diff(par()$usr[1:2])/diff(par()$usr[3:4])
	if(is.null(node)){
		if(!interactive){
			cat("No nodes provided. Setting interactive mode to TRUE.\n")
			interactive<-TRUE
		}
		node<-vector(length=length(text))
	}
	for(i in 1:length(text)){
		if(interactive){
			cat(paste("Click on the node you would like to label ",
				text[i],".\n",sep=""))
			flush.console()
			ii<-getnode(env=obj)
			node[i]<-ii
		} else ii<-node[i]
		if(shape=="circle")
			draw.circle(obj$xx[ii],obj$yy[ii],rad,col=bg)
		else if(shape=="ellipse")
			draw.ellipse(obj$xx[ii],obj$yy[ii],0.8*w[i],h,
				col=bg)
		else if(shape=="rect")
			rect(xleft=obj$xx[ii]-0.5*rect.exp*w[i],
				ybottom=obj$yy[ii]-0.5*rect.exp*h,
				xright=obj$xx[ii]+0.5*rect.exp*w[i],
				ytop=obj$yy[ii]+0.5*rect.exp*h,col=bg,
				ljoin=1)
		text(obj$xx[ii],obj$yy[ii],label=text[i],cex=cex)
	}
	invisible(node)
}

## convert object of class "birthdeath" into birth & death rates
bd<-function(x){
	if(!inherits(x,"birthdeath")) stop("x should be an object of class 'birthdeath'")
	b<-x$para[2]/(1-x$para[1])
	d<-b-x$para[2]
	setNames(c(b,d),c("b","d"))
}

## compute AIC weights
aic.w<-function(aic){
	d.aic<-aic-min(aic)
	x<-exp(-1/2*d.aic)/sum(exp(-1/2*d.aic))
	class(x)<-"aic.w"
	x
}

print.aic.w<-function(x,...){
	if(hasArg(signif)) signif<-list(...)$signif
	else signif<-8
	print(round(unclass(x),signif))
}

## function to compute all paths towards the tips from a node
## written by	Liam J. Revell
node.paths<-function(tree,node){
	d<-Descendants(tree,node,"children")
	paths<-as.list(d)
	while(any(d>Ntip(tree))){
		jj<-1
		new.paths<-list()
		for(i in 1:length(paths)){
			if(paths[[i]][length(paths[[i]])]<=Ntip(tree)){ 
				new.paths[[jj]]<-paths[[i]]
				jj<-jj+1
			} else {
				ch<-Descendants(tree,paths[[i]][length(paths[[i]])],
					"children")
				for(j in 1:length(ch)){
					new.paths[[jj]]<-c(paths[[i]],ch[j])
					jj<-jj+1
				}
			}
		}
		paths<-new.paths
		d<-sapply(paths,function(x) x[length(x)])
	}
	paths
}

## function to compute a modification of Grafen's edge lengths
## written by Liam J. Revell 2016
modified.Grafen<-function(tree,power=2){
	max.np<-function(tree,node){
		np<-node.paths(tree,node)
		if(length(np)>0) max(sapply(np,length)) else 0
	}
	nn<-1:(Ntip(tree)+tree$Nnode)
	h<-sapply(nn,max.np,tree=tree)+1
	h<-(h/max(h))^power
	edge.length<-vector()
	for(i in 1:nrow(tree$edge)) 
		edge.length[i]<-diff(h[tree$edge[i,2:1]])
	tree$edge.length<-edge.length
	tree
}

## function to compute all rotations
## written by Liam J. Revell 2016
allRotations<-function(tree){
	if(!is.binary(tree)){
		was.binary<-FALSE
		if(is.null(tree$edge.length)){ 
			tree<-compute.brlen(tree)
			had.edge.lengths<-FALSE
		} else had.edge.lengths<-TRUE
		tree<-multi2di(tree)
	} else was.binary<-TRUE
	nodes<-1:tree$Nnode+Ntip(tree)
	trees<-vector(mode="list",length=2^length(nodes))
	ii<-2
	trees[[1]]<-tree
	for(i in 1:length(nodes)){
		N<-ii-1
		for(j in 1:N){
			trees[[ii]]<-rotate(trees[[j]],nodes[i])
			ii<-ii+1
		}
	}
	trees<-lapply(trees,untangle,"read.tree")
	if(!was.binary){
		trees<-lapply(trees,di2multi)
		if(!had.edge.lengths) trees<-lapply(trees,
			function(x){
				x$edge.length<-NULL
				x
			})
	}
	class(trees)<-"multiPhylo"
	trees
}

## function to rotate a multifurcation in all possible ways
## written by Liam J. Revell 2016
rotate.multi<-function(tree,node){
	kids<-Children(tree,node)
	if(length(kids)>2){
		ii<-sapply(kids,function(x,y) which(y==x),y=tree$edge[,2])
		jj<-permn(ii)
		foo<-function(j,i,t){
			t$edge[i,]<-t$edge[j,]
			if(!is.null(t$edge.length))
				t$edge.length[i]<-t$edge.length[j]
			untangle(t,"read.tree")
		}
		obj<-lapply(jj[2:length(jj)],foo,i=ii,t=tree)
		class(obj)<-"multiPhylo"
	} else obj<-untangle(rotate(tree,node),"read.tree")
	obj
}

## wrapper for bind.tree that takes objects of class "simmap"
## written by Liam J. Revell 2016
bind.tree.simmap<-function(x,y,where="root"){
	x<-reorder(x)
	y<-reorder(y)
	rootx<-x$edge[1,1]
	rooty<-y$edge[1,1]
	xy<-read.tree(text=write.tree(bind.tree(x,y,where)))
	Mx<-rbind(matchLabels(x,xy),matchNodes(x,xy,"distances"))
	My<-rbind(matchLabels(y,xy),matchNodes(y,xy,"distances"))
	if(where!="root"&&where<=Ntip(x))
		Mx[which(is.na(Mx[,2])),2]<-findMRCA(xy,y$tip.label)
	xy$maps<-vector(mode="list",length=nrow(xy$edge))
	ix<-sapply(Mx[-which(Mx[,1]==rootx),1],
		function(x,y) which(y==x),y=x$edge[,2])
	ixy<-sapply(Mx[-which(Mx[,1]==rootx),2],
		function(x,y) which(y==x),y=xy$edge[,2])
	xy$maps[ixy]<-x$maps[ix]
	iy<-sapply(My[-which(My[,1]==rooty),1],
		function(x,y) which(y==x),y=y$edge[,2])
	ixy<-sapply(My[-which(My[,1]==rooty),2],
		function(x,y) which(y==x),y=xy$edge[,2])
	xy$maps[ixy]<-y$maps[iy]
	xy$mapped.edge<-makeMappedEdge(xy$edge,xy$maps)
	ns<-c(setNames(getStates(xy,"tips"),1:Ntip(xy)),
		getStates(xy,"nodes"))
	xy$node.states<-cbind(ns[as.character(xy$edge[,1])],
		ns[as.character(xy$edge[,2])])
	xy$states<-getStates(xy,"tips")
	attr(xy,"class")<-c("simmap",class(xy))
	xy
}

## generic function to convert an object of class "simmap" to "phylo"
## written by Liam J. Revell 2016
as.phylo.simmap<-function(x,...){
	x$maps<-NULL
	x$mapped.edge<-NULL
	if(!is.null(x$node.states)) x$node.states<-NULL
	if(!is.null(x$states)) x$states<-NULL
	if(!is.null(x$Q)) x$Q<-NULL
	if(!is.null(x$logL)) x$logL<-NULL
	if(!is.null(attr(x,"map.order"))) attr(x,"map.order")<-NULL
	class(x)<-setdiff(class(x),"simmap")
	x
}

## generic function to convert an object of class "multiSimmap" to "multiPhylo"
## written by Liam J. Revell 2016
as.multiPhylo.multiSimmap<-function(x,...){
	obj<-lapply(x,as.phylo)
	class(obj)<-setdiff(class(x),"multiSimmap")
	obj
}

## generic function to convert object of class "phylo" to "multiPhylo"
## written by Liam J. Revell 2016
as.multiPhylo.phylo<-function(x,...){
	obj<-list(x)
	class(obj)<-"multiPhylo"
	obj
}

as.multiPhylo<-function(x,...){
	if (identical(class(x),"multiPhylo")) return(x)
	UseMethod("as.multiPhylo")
}

## get mapped states
## written by Liam J. Revell 2016
mapped.states<-function(tree,...){
	if(!(inherits(tree,"simmap")||inherits(tree,"multiSimmap")))
		stop("tree should be an object of class \"simmap\" or \"multiSimmap\".")
	else {
		if(inherits(tree,"simmap")){
			if(!is.null(tree$mapped.edge)) 
				obj<-sort(colnames(tree$mapped.edge))
			else 
				obj<-sort(unique(unlist(lapply(tree$maps,function(x) names(x)))))
		} else if(inherits(tree,"multiSimmap")) {
			obj<-sapply(tree,mapped.states,...)
		}
	}
	obj
}

## match labels between trees (equivalent to matchNodes)
## written by Liam J. Revell 2016
matchLabels<-function(tr1,tr2){
	foo<-function(x,y) if(length(obj<-which(y==x))>0) obj else NA
	M<-cbind(1:Ntip(tr1),sapply(tr1$tip.label,foo,y=tr2$tip.label))
	colnames(M)<-c("tr1","tr2")
	M
}

## compute the probability of states changes along edges of the tree
## written by Liam J. Revell 2015
edgeProbs<-function(trees){
	if(!inherits(trees,"multiSimmap")) stop("trees should be an object of class \"multiSimmap\".")
	SS<-sapply(trees,getStates,"tips")
	states<-sort(unique(as.vector(SS)))
	m<-length(states)
	TT<-sapply(states,function(x,y) sapply(y,paste,x,sep="->"),
		y=states)
	nn<-c(TT[upper.tri(TT)],TT[lower.tri(TT)])
	## this function computes for a given edge
	fn<-function(edge,trees,states){
		obj<-sapply(trees,function(x,e,s) 
		if(names(x$maps[[e]])[1]==
			s[1]&&names(x$maps[[e]])[length(x$maps[[e]])]==s[2]) TRUE
		else FALSE,e=edge,s=states)
		sum(obj)/length(obj)
	}
	edge.probs<-matrix(0,nrow(trees[[1]]$edge),m,
		dimnames=list(apply(trees[[1]]$edge,1,paste,collapse=","),nn))
	k<-1
	for(i in 1:m) for(j in 1:m){
		if(i!=j){ 
			edge.probs[,k]<-sapply(1:nrow(trees[[1]]$edge),fn,
				trees=trees,states=states[c(i,j)])
			k<-k+1
		}
	}
	edge.probs<-cbind(edge.probs,1-rowSums(edge.probs))
	colnames(edge.probs)[ncol(edge.probs)]<-"no change"
	edge.probs
}

## get a position in the tree interactively
## written by Liam J. Revell 2015, 2016, 2020
get.treepos<-function(message=TRUE,...){
	obj<-get("last_plot.phylo",envir=.PlotPhyloEnv)
	if(obj$type=="phylogram"&&obj$direction=="rightwards"){
		if(message){ 
			cat("Click on the tree position you want to capture...\n")
			flush.console()
		}
		if(hasArg(x)) x<-list(...)$x
		else x<-NULL
		if(hasArg(y)) y<-list(...)$y
		else y<-NULL
		if(is.null(x)||is.null(y)){
			x<-unlist(locator(1)) 	
			y<-x[2] 	
			x<-x[1]
		}
		d<-pos<-c()
		for(i in 1:nrow(obj$edge)){
			x0<-obj$xx[obj$edge[i,]]
			y0<-obj$yy[obj$edge[i,2]]
			if(x<x0[1]||x>x0[2]){
				d[i]<-min(dist(rbind(c(x,y),c(x0[1],y0))),
					dist(rbind(c(x,y),c(x0[2],y0))))
				pos[i]<-if(x>x0[2]) 0 else diff(obj$xx[obj$edge[i,]])
			} else {
				d[i]<-abs(y0-y)
				pos[i]<-obj$xx[obj$edge[i,2]]-x
			}
		}
		ii<-which(d==min(d))
		## check to make sure the root is not closer:
		root.d<-dist(rbind(c(x,y),c(obj$xx[obj$Ntip+1],obj$yy[obj$Ntip+1])))
		if(root.d<min(d)) return(list(where=obj$Ntip+1,pos=0))
		else list(where=obj$edge[ii,2],pos=pos[ii])
	} else stop("Does not work for the plotted tree type.")
}

## fastDist: uses fastHeight to compute patristic distance between a pair of species
fastDist<-function(tree,sp1,sp2){
	if(!inherits(tree,"phylo")) stop("tree should be an object of class \"phylo\".")
	if(is.null(tree$edge.length)) stop("tree should have edge lengths.")
	if(sp1==sp2) 0
	else fastHeight(tree,sp1,sp1)+fastHeight(tree,sp2,sp2)-
		2*fastHeight(tree,sp1,sp2)
}

## function reorders simmap tree
## written Liam Revell 2011, 2013, 2015, 2019
reorderSimmap<-function(tree,order="cladewise",index.only=FALSE,...){
	if(!inherits(tree,"phylo")) stop("tree should be an object of class \"phylo\".")
	ii<-reorder.phylo(tree,order,index.only=TRUE,...)
	if(!index.only){
		if(inherits(ii,"phylo")) ii<-whichorder(ii$edge[,2],
			tree$edge[,2]) ## bug workaround
		tree$edge<-tree$edge[ii,]
		tree$edge.length<-tree$edge.length[ii]
		if(!is.null(tree$maps)){
			tree$maps<-tree$maps[ii]
			tree$mapped.edge<-tree$mapped.edge[ii,,drop=FALSE]
		}
		attr(tree,"order")<-order
		return(tree)
	} else return(ii)
}

## S3 reorder method for objects of class "simmap"
reorder.simmap<-function(x,...) reorderSimmap(x,...)

# function whichorder
# written by Liam Revell 2011, 2013, 2015
whichorder<-function(x,y) sapply(x,function(x,y) which(x==y),y=y)

# function returns random state with probability given by y
# written by Liam J. Revell 2013, 2015
rstate<-function(y){
	if(length(y)==1) return(names(y)[1])
	else {
		p<-y/sum(y)
		if(any(p<0)){ 
			warning("Some probabilities (slightly?) < 0. Setting p < 0 to zero.")
			p[p<0]<-0
		}
		return(names(which(rmultinom(1,1,p)[,1]==1)))
	}
}

## mark the changes on a plotted "simmap" object
## written by Liam J. Revell 2015
markChanges<-function(tree,colors=NULL,cex=1,lwd=2,plot=TRUE){
	states<-sort(unique(getStates(tree)))
	if(is.null(colors)) colors<-setNames(palette()[1:length(states)],
		states)
	obj<-get("last_plot.phylo",envir=.PlotPhyloEnv)
	nc<-sapply(tree$maps,length)-1
	ii<-which(nc>0)
	nc<-nc[ii]
	xx<-yy<-vector()
	for(i in 1:length(ii)){
		for(j in 1:nc[i]){
			ss<-names(tree$maps[[ii[i]]])[j+1]
			mm<-tree$edge[ii[i],1]
			dd<-tree$edge[ii[i],2]
			x<-rep(obj$xx[mm]+cumsum(tree$maps[[ii[i]]])[j],2)
			y<-c(obj$yy[dd]-0.5*mean(strheight(LETTERS)*cex),
				obj$yy[dd]+0.5*mean(strheight(LETTERS)*cex))
			if(plot) lines(x,y,lwd=lwd,col=colors[ss],lend=2)
			xx<-c(xx,setNames(x[1],
				paste(names(tree$maps[[ii[i]]])[j:(j+1)],
				collapse="->")))
			yy<-c(yy,mean(y))
		}
	}
	XY<-cbind(xx,yy)
	colnames(XY)<-c("x","y")
	invisible(XY)
}

## function to label clades
## written by Liam J. Revell 2014, 2015, 2022
cladelabels<-function(tree=NULL,text,node,offset=NULL,wing.length=NULL,cex=1,
	orientation="vertical"){
	lastPP<-get("last_plot.phylo",envir=.PlotPhyloEnv)
	if(is.null(tree)){
		if(is.null(wing.length)) wing.length<-1
		if(is.null(offset)) offset<-8
		tree<-list(edge=lastPP$edge,
			tip.label=1:lastPP$Ntip,
			Nnode=lastPP$Nnode)
		H<-matrix(lastPP$xx[tree$edge],nrow(tree$edge),2)
		tree$edge.length<-H[,2]-H[,1]
		class(tree)<-"phylo"
	}
	if(is.null(offset)) offset<-0.5
	xx<-mapply(labelSubTree,node,text,
		MoreArgs=list(tree=tree,pp=lastPP,offset=offset,wl=wing.length,cex=cex,
		orientation=orientation),SIMPLIFY=FALSE)
}

## internal function used by cladelabels
## written by Liam J. Revell 2014, 2015, 2022
labelSubTree<-function(tree,nn,label,pp,offset,wl,cex,orientation){
	if(is.null(wl)) wl<-1
	tree<-reorder(tree)
	tips<-getDescendants(tree,nn)
	tips<-tips[tips<=Ntip(tree)]
	ec<-0.7 ## expansion constant
	sw<-pp$cex*max(strwidth(tree$tip.label[tips]))
	sh<-pp$cex*max(strheight(tree$tip.label))
	cw<-mean(strwidth(LETTERS)*cex)	
	h<-max(sapply(tips,function(x,tree)
		nodeHeights(tree)[which(tree$edge[,2]==x),2],
		tree=tree))+sw+offset*cw
	y<-range(pp$yy[tips])
	lines(c(h,h),y+ec*c(-sh,sh),col=par()$fg)
	lines(c(h-wl*cw,h),
		c(y[1]-ec*sh,y[1]-ec*sh),col=par()$fg)
	lines(c(h-wl*cw,h),
		c(y[2]+ec*sh,y[2]+ec*sh),col=par()$fg)
	text(h+cw,mean(y),
		label,srt=if(orientation=="horizontal") 0 else 90,
		adj=if(orientation=="horizontal") 0 else 0.5,cex=cex,
		col=par()$col.lab)
	list(xx=c(h,h),yy=y+ec*c(-sh,sh))
}

## get all the extant/extinct tip names
## written by Liam J. Revell 2012, 2015

getExtant<-function(tree,tol=1e-8){
	if(!inherits(tree,"phylo")) stop("tree should be object of class \"phylo\".")
	H<-nodeHeights(tree)
	tl<-max(H)
	x<-which(H[,2]>=(tl-tol))
	y<-tree$edge[x,2]
	y<-y[y<=Ntip(tree)]
	z<-tree$tip.label[y]
	return(z)
}

getExtinct<-function(tree,tol=1e-8) setdiff(tree$tip.label,getExtant(tree,tol))

# function splits tree at split
# written by Liam Revell 2011, 2014, 2015, 2020

splitTree<-function(tree,split){
	if(!inherits(tree,"phylo")) stop("tree should be an object of class \"phylo\".")
	if(split$node>Ntip(tree)){
		# first extract the clade given by shift$node
		tr2<-extract.clade(tree,node=split$node)
		tr2$root.edge<-tree$edge.length[which(tree$edge[,2]==split$node)]-split$bp
		#now remove tips in tr2 from tree
		tr1<-drop.clade(tree,tr2$tip.label)
		nn<-if(!is.null(tree$node.label)) c(tree$node.label,"NA") else "NA"
		tr1$tip.label[which(tr1$tip.label%in%nn)]<-"NA"
		tr1$edge.length[match(which(tr1$tip.label=="NA"),tr1$edge[,2])]<-split$bp
	} else {
		# first extract the clade given by shift$node
		tr2<-list(edge=matrix(c(2L,1L),1,2),tip.label=tree$tip.label[split$node],edge.length=tree$edge.length[which(tree$edge[,2]==split$node)]-split$bp,Nnode=1L)
		class(tr2)<-"phylo"
		# now remove tip in tr2 from tree
		tr1<-tree
		tr1$edge.length[match(which(tr1$tip.label==tr2$tip.label[1]),tr1$edge[,2])]<-split$bp
		tr1$tip.label[which(tr1$tip.label==tr2$tip.label[1])]<-"NA"
	}
	trees<-list(tr1,tr2)
	class(trees)<-"multiPhylo"
	trees
}

# function drops entire clade
# written by Liam Revell 2011, 2015, 2020

drop.clade<-function(tree,tip){
	## step 1, check to make sure tips form a monophyletic clade
	node<-getMRCA(tree,tip)
	desc<-getDescendants(tree,node)
	chk<-tree$tip.label[desc[desc<=Ntip(tree)]]
	if(!setequal(tip,chk)){
		cat("Caution: Species in tip do not form a monophyletic clade.\n")
		cat("				 Pruning all tips descended from ancestor.\n\n")
		tip<-chk
	}
	## step 2, find all edges in the clade & set them to zero length
	ee<-sapply(desc,function(node,edge) which(edge==node), 
		edge=tree$edge[,2])
	tree$edge.length[ee]<-0
	## step 3, bind a tip labeled "NA" to the node
	tree<-bind.tip(tree,"NA",0,node)
	## step 4, prune the other tips
	tree<-drop.tip(tree,tip)
	## step 5, return tree
	tree
}

## function to re-root a phylogeny along an edge
## written by Liam J. Revell 2011-2016, 2019, 2020

reroot<-function(tree,node.number,position=NULL,interactive=FALSE,...){
	if(!inherits(tree,"phylo")) stop("tree should be an object of class \"phylo\".")
	if(interactive){
		plotTree(tree,...)
		cat("Click where you would like re-root the plotted tree\n")
		flush.console()
		obj<-get.treepos(message=FALSE)
		node.number<-obj$where
		position<-tree$edge.length[which(tree$edge[,2]==node.number)]-obj$pos
	}
	if(node.number==(Ntip(tree)+1))
		cat("Note: you chose to re-root the tree at it's current root.\n")
	if(is.null(position)){
		if(node.number==(Ntip(tree)+1)) position=0
		else position<-tree$edge.length[which(tree$edge[,2]==node.number)]
	} else {
		if(node.number==(Ntip(tree)+1))
			cat("			A value of position != 0 has been reset to zero.\n")
	}
	if(hasArg(edgelabel)) edgelabel<-list(...)$edgelabel
	else edgelabel<-FALSE
	if(node.number!=(Ntip(tree)+1)){
		tt<-splitTree(tree,list(node=node.number,bp=position))
		p<-tt[[1]]
		d<-tt[[2]]
		tip<-if(length(which(p$tip.label=="NA"))>0) "NA" else p$tip.label[which(p$tip.label%in%tree$node.label)]
		p<-ape::root.phylo(p,outgroup=tip,resolve.root=TRUE,edgelabel=edgelabel)
		bb<-which(p$tip.label==tip)
		p$tip.label[bb]<-"NA"
		ee<-p$edge.length[which(p$edge[,2]==bb)]
		p$edge.length[which(p$edge[,2]==bb)]<-0
		cc<-p$edge[which(p$edge[,2]==bb),1]
		dd<-setdiff(p$edge[which(p$edge[,1]==cc),2],bb)
		p$edge.length[which(p$edge[,2]==dd)]<-p$edge.length[which(p$edge[,2]==dd)]+ee
		obj<-paste.tree(p,d)
	} else obj<-tree
	if(interactive) plotTree(obj,...)
	obj
}

## function to ladderize phylogeny with mapped discrete character
## written by Liam J. Revell 2014, 2015, 2019, 2023

ladderize.simmap<-function(tree,right=TRUE){
	if(!inherits(tree,"simmap")){
		if(!inherits(tree,"phylo")) 
			stop("tree should be an object of class \"phylo\".")
		else {
			cat("Do not detect a mapped character. Using ape::ladderize.\n")
			obj<-ladderize(tree,right=right)
		}
	} else {
		obj<-read.tree(text=write.tree(ladderize(tree,right=right)))
		rN<-Ntip(obj)+1
		T<-cbind(1:Ntip(obj),sapply(obj$tip.label,
			function(x,y) which(y==x),y=tree$tip.label))
		N<-matchNodes(obj,as.phylo(tree))
		M<-rbind(T,N[N[,1]!=rN,])
		ii<-sapply(M[,1],function(x,y) which(y==x),y=obj$edge[,2])
		jj<-sapply(M[,2],function(x,y) which(y==x),y=tree$edge[,2])
		obj$maps<-vector(length=nrow(tree$edge),mode="list")
		obj$mapped.edge<-matrix(NA,nrow(tree$edge),ncol(tree$mapped.edge),
			dimnames=list(apply(tree$edge,1,paste,collapse=","),
			colnames(tree$mapped.edge)))
		if(!is.null(tree$states)) 
			obj$states<-tree$states[sapply(obj$tip.label,
				function(x,y) which(y==x),y=tree$tip.label)]
		if(!is.null(tree$node.states)) obj$node.states<-matrix(NA,nrow(tree$edge),2)
		for(i in 1:length(ii)){
			obj$maps[[ii[i]]]<-tree$maps[[jj[i]]]
			obj$mapped.edge[ii[i],]<-tree$mapped.edge[jj[i],]
			if(!is.null(tree$node.states)) obj$node.states[ii[i],]<-
				tree$node.states[jj[i],]
		}
		class(obj)<-c("simmap","phylo")
	}
	obj
}

## for backward compatibility

repPhylo<-function(tree,times) rep(tree,times)

## S3 method rep for objects of class "phylo" and "multiPhylo"
## written by Liam J. Revell 2014, 2021

rep.phylo<-function(x,...){
	if(hasArg(times)) times<-list(...)$times
	else times<-(...)[[1]]
	for(i in 1:times) obj<-if(i==1) x else if(i==2) c(obj,x) else c(unclass(obj),list(x))
	class(obj)<-"multiPhylo"
	obj
}

rep.multiPhylo<-function(x,...){
	if(hasArg(times)) times<-list(...)$times
	else times<-(...)[[1]]
	for(i in 1:times) obj<-if(i==1) x else if(i>=2) c(obj,x)
	class(obj)<-"multiPhylo"
	obj
}

## function to drop one or more tips from a tree but retain all ancestral nodes as singletons
## written by Liam J. Revell 2014, 2015, 2023
drop.tip.singleton<-function(phy,tip,...){
	if(!inherits(phy,"singleton")&&!inherits(phy,"phylo")) 
		stop("phy should be an object of class \"singleton\".")
	N<-Ntip(phy)
	m<-length(tip)
	ii<-sapply(tip,function(x,y) which(y==x),y=phy$tip.label)
	phy$tip.label<-phy$tip.label[-ii]
	ii<-sapply(ii,function(x,y) which(y==x),y=phy$edge[,2])
	phy$edge<-phy$edge[-ii,]
	phy$edge.length<-phy$edge.length[-ii]
	phy$edge[phy$edge<=N]<-as.integer(rank(phy$edge[phy$edge<=N]))
	phy$edge[phy$edge>N]<-phy$edge[phy$edge>N]-m
	N<-N-m
	if(any(sapply(phy$edge[phy$edge[,2]>N,2],"%in%",phy$edge[,1])==FALSE)) internal<-TRUE
	else internal<-FALSE
	while(internal){
		ii<-which(sapply(phy$edge[,2],"%in%",c(1:N,phy$edge[,1]))==FALSE)[1]
		nn<-phy$edge[ii,2]
		phy$edge<-phy$edge[-ii,]
		phy$edge.length<-phy$edge.length[-ii]
		phy$edge[phy$edge>nn]<-phy$edge[phy$edge>nn]-1
		phy$Nnode<-phy$Nnode-length(ii)
		if(any(sapply(phy$edge[phy$edge[,2]>N,2],"%in%",phy$edge[,1])==FALSE)) internal<-TRUE
		else internal<-FALSE
	}
	class(phy)<-union("singleton",class(phy))
	phy
}

## S3 print method for object of class 'describe.simmap'
## written by Liam J. Revell 2014, 2015
print.describe.simmap<-function(x,...){
	if(inherits(x$tree,"multiPhylo")){
		cat(paste(length(x$tree),"trees with a mapped discrete character with states:\n",paste(colnames(x$ace),collapse=", "),"\n\n"))
		cat(paste("trees have",colMeans(x$count)["N"],"changes between states on average\n\n"))
		cat(paste("changes are of the following types:\n"))
		aa<-t(as.matrix(colMeans(x$count)[2:ncol(x$count)]))
		rownames(aa)<-"x->y"
		print(aa)
		cat(paste("\nmean total time spent in each state is:\n"))
		print(matrix(c(colMeans(x$times),colMeans(x$times[,1:ncol(x$times)]/x$times[,ncol(x$times)])),2,ncol(x$times),byrow=TRUE,
			dimnames=list(c("raw","prop"),c(colnames(x$times)))))
		cat("\n")
	} else if(inherits(x$tree,"phylo")){
		cat(paste("1 tree with a mapped discrete character with states:\n",paste(colnames(x$Tr),collapse=", "),"\n\n"))
		cat(paste("tree has",x$N,"changes between states\n\n"))
		cat(paste("changes are of the following types:\n"))
		print(x$Tr)
		cat(paste("\nmean total time spent in each state is:\n"))
		print(x$times)
		cat("\n")
	}
}

## S3 plot method for object of class 'describe.simmap'
## written by Liam J. Revell 2014, 2015, 2020
plot.describe.simmap<-function(x,...){
	if(hasArg(lwd)) lwd<-list(...)$lwd
	else lwd<-2
	if(hasArg(cex)) cex<-list(...)$cex
	else cex<-c(0.6,0.4)
	if(length(cex)==1) cex<-rep(cex,2)
	if(hasArg(type)) type<-list(...)$type
	else type<-"phylogram"
	if(hasArg(offset)) offset<-list(...)$offset
	else offset<-0
	if(hasArg(fsize)) fsize<-list(...)$fsize
	else fsize<-1
	if(inherits(x$tree,"multiPhylo")){
		states<-colnames(x$ace)
		if(hasArg(colors)) colors<-list(...)$colors
		else colors<-setNames(palette()[1:length(states)],states)
		plotTree(if(is.null(x$ref.tree)) x$tree[[1]] else x$ref.tree,lwd=lwd,
			offset=cex[2]+offset*fsize,...)
		nodelabels(pie=x$ace,piecol=colors[colnames(x$ace)],cex=cex[1])
		if(!is.null(x$tips)) tips<-x$tips 
		else tips<-to.matrix(getStates(x$tree[[1]],"tips"),seq=states) 
		tiplabels(pie=tips[if(is.null(x$ref.tree)) x$tree[[1]]$tip.label else 
			x$ref.tree$tip.label,],piecol=colors[colnames(tips)],cex=cex[2])
	} else if(inherits(x$tree,"phylo")){
		states<-colnames(x$Tr)
		if(hasArg(colors)) colors<-list(...)$colors
		else colors<-setNames(palette()[1:length(states)],states)
		plotSimmap(x$tree,lwd=lwd,colors=colors,type=type,offset=offset)
	}
}

## function finds the height of a given node
## written by Liam Revell 2014, 2015, 2016
nodeheight<-function(tree,node,...){
	if(hasArg(root.edge)) root.edge<-list(...)$root.edge
	else root.edge<-FALSE
	if(root.edge) ROOT<-if(!is.null(tree$root.edge)) tree$root.edge else 0
	else ROOT<-0 
	if(!inherits(tree,"phylo")) stop("tree should be an object of class \"phylo\".")
	if(node==(Ntip(tree)+1)) h<-0
	else {
		a<-setdiff(c(getAncestors(tree,node),node),Ntip(tree)+1)
		h<-sum(tree$edge.length[sapply(a,function(x,e) which(e==x),e=tree$edge[,2])])
	}
	h+ROOT
}

# fast pairwise MRCA function
# written by Liam Revell 2012, 2014, 2015
fastMRCA<-function(tree,sp1,sp2){
	if(!inherits(tree,"phylo")) stop("tree should be an object of class \"phylo\".")
	x<-match(sp1,tree$tip.label)
	y<-match(sp2,tree$tip.label)
	a<-c(x,getAncestors(tree,x))
	b<-c(y,getAncestors(tree,y))
	z<-a%in%b
	return(a[min(which(z))])
}

## function to find the height of the MRCA of sp1 & sp2
## written by Liam J. Revell 2014, 2015
fastHeight<-function(tree,sp1,sp2){
	if(!inherits(tree,"phylo")) stop("tree should be an object of class \"phylo\".")
	if(is.null(tree$edge.length)) stop("tree should have edge lengths.")
	sp1<-which(tree$tip.label==sp1)
	sp2<-which(tree$tip.label==sp2)
	a1<-c(sp1,getAncestors(tree,sp1))
	a2<-c(sp2,getAncestors(tree,sp2))
	a12<-intersect(a1,a2)
	if(length(a12)>1){ 
		a12<-a12[2:length(a12)-1]
		h<-sapply(a12,function(i,tree) tree$edge.length[which(tree$edge[,2]==i)],tree=tree)
		return(sum(h))
	} else return(0)
}

## function gets ancestor node numbers, to be used internally by 
## written by Liam J. Revell 2014
getAncestors<-function(tree,node,type=c("all","parent")){
	if(!inherits(tree,"phylo")) stop("tree should be an object of class \"phylo\".")
	type<-type[1]
	if(type=="all"){
		aa<-vector()
		rt<-Ntip(tree)+1
		currnode<-node
		while(currnode!=rt){
			currnode<-getAncestors(tree,currnode,"parent")
			aa<-c(aa,currnode)
		}
		return(aa)
	} else if(type=="parent"){
		aa<-tree$edge[which(tree$edge[,2]==node),1]
		return(aa)
	} else stop("do not recognize type")
}

## function for midpoint rooting
## written by Liam J. Revell 2014
## (being deprecated out in 2023)
midpoint_root<-function(tree){
	D<-cophenetic(tree)
	dd<-max(D)
	ii<-which(D==dd)[1]
	ii<-c(ceiling(ii/nrow(D)),ii%%nrow(D))
	if(ii[2]==0) ii[2]<-nrow(D)
	spp<-rownames(D)[ii]
	nn<-which(tree$tip.label==spp[2])
	tree<-reroot(tree,nn,tree$edge.length[which(tree$edge[,2]==nn)])
	aa<-getAncestors(tree,which(tree$tip.label==spp[1]))
	D<-c(0,dist.nodes(tree)[which(tree$tip.label==spp[1]),aa])
	names(D)[1]<-which(tree$tip.label==spp[1])
	i<-0
	while(D[i+1]<(dd/2)) i<-i+1
	tree<-reroot(tree,as.numeric(names(D)[i]),D[i+1]-dd/2)
	tree
}

midpoint.root<-function(tree,node.labels="support",...) 
	phangorn::midpoint(tree,node.labels="support",...)

# function computes phylogenetic variance-covariance matrix, including for internal nodes
# written by Liam J. Revell 2011, 2013, 2014, 2015
vcvPhylo<-function(tree,anc.nodes=TRUE,...){
	if(!inherits(tree,"phylo")) stop("tree should be an object of class \"phylo\".")
	# get & set optional arguments
	if(hasArg(internal)) internal<-list(...)$internal
	else internal<-anc.nodes
	if(internal!=anc.nodes){ 
		message(paste("arguments \"internal\" and \"anc.nodes\" are synonyms; setting internal =",anc.nodes))
		internal<-anc.nodes
	}
	if(hasArg(model)) model<-list(...)$model
	else model<-"BM"
	if(hasArg(tol)) tol<-list(...)$tol
	else tol<-1e-12
	if(model=="OU"){
		if(hasArg(alpha)) alpha<-list(...)$alpha
		else alpha<-0
	}
	if(model=="OU"&&alpha<tol) model<-"BM"
	if(model=="lambda"){
		if(hasArg(lambda)){ 
			lambda<-list(...)$lambda
			tree<-lambdaTree(tree,lambda)
		} else model<-"BM"
		model<-"BM"
	}
	if(model=="EB"){
		if(hasArg(r)){
			r<-list(...)$r
			tree<-ebTree(tree,r)
		} else model<-"BM"
	}
	# done settings
	n<-Ntip(tree)
	h<-nodeHeights(tree)[order(tree$edge[,2]),2]
	h<-c(h[1:n],0,h[(n+1):length(h)])
	M<-mrca(tree,full=internal)[c(1:n,internal*(n+2:tree$Nnode)),c(1:n,internal*(n+2:tree$Nnode))]
	C<-matrix(h[M],nrow(M),ncol(M))
	if(internal) rownames(C)<-colnames(C)<-c(tree$tip.label,n+2:tree$Nnode)
	else rownames(C)<-colnames(C)<-tree$tip.label
	if(model=="OU"){
		D<-dist.nodes(tree)
		rownames(D)[1:n]<-colnames(D)[1:n]<-tree$tip.label
		D<-D[rownames(C),colnames(C)]
		# not sure 
		C<-(1-exp(-2*alpha*C))*exp(-alpha*D)/(2*alpha) # Hansen (2007)
		# C<-(1-exp(-2*alpha*C))*exp(-2*alpha*(1-C))/(2*alpha) # Butler & King (2004)
	}
	return(C)
}

## simplified lambdaTree to be used internally by vcvPhylo
## written by Liam J. Revell 2014
lambdaTree<-function(tree,lambda){
	ii<-which(tree$edge[,2]>Ntip(tree))
	H1<-nodeHeights(tree)
	tree$edge.length[ii]<-lambda*tree$edge.length[ii]
	H2<-nodeHeights(tree)
	tree$edge.length[-ii]<-tree$edge.length[-ii]+		 H1[-ii,2]-H2[-ii,2]
	tree
}

## di2multi method for tree with mapped state
## written by Liam J. Revell 2013, 2015, 2016
di2multi.simmap<-function(phy,...){
	if(hasArg(tol)) tol<-list(...)$tol
	else tol<-1e-08
	if(!inherits(phy,"phylo")) stop("tree should be an object of class \"phylo\".")
	if(is.null(phy$maps)){
		cat("Warning: tree does not contain mapped state. Using di2multi.\n")
		return(di2multi(phy,tol))
	}
	N<-length(phy$tip.label)
	n<-length(intersect(which(phy$edge.length<tol),which(phy$edge[,2]>N)))
	if(n==0) return(phy)
	edge<-phy$edge
	edge[edge>N]<--edge[edge>N]+N
	edge.length<-phy$edge.length
	maps<-phy$maps
	Nnode<-phy$Nnode
	for(i in 1:n){
		ii<-intersect(which(edge.length<tol),which(edge[,2]<0))[1]
		node<-edge[ii,2]
		edge[edge==node]<-edge[ii,1]
		jj<-which(apply(edge,1,function(x) x[1]==x[2]))[1]
		edge<-edge[-jj,]
		edge.length<-edge.length[-jj]
		maps<-maps[-jj]	
		Nnode<-Nnode-1
	}
	nn<-sort(unique(edge[edge<0]),decreasing=TRUE)
	mm<-1:Nnode+N
	for(i in 1:length(edge)) if(edge[i]%in%nn) edge[i]<-mm[which(nn==edge[i])]
	mapped.edge<-makeMappedEdge(edge,maps)
	tt<-list(edge=edge,Nnode=Nnode,tip.label=phy$tip.label,edge.length=edge.length,
		maps=maps,mapped.edge=mapped.edge)
	class(tt)<-"phylo"
	if(!is.null(attr(phy,"order"))) attr(tt,"order")<-attr(phy,"order")
	if(!is.null(phy$node.states)) tt$node.states<-getStates(tt,"nodes")
	if(!is.null(phy$states)) tt$states<-getStates(tt,"tips")
	return(tt)
}

# returns the heights of each node
# written by Liam J. Revell 2011, 2012, 2013, 2015, 2016
# modified by Klaus Schliep 2017
nodeHeights<-function(tree,...){
		if(hasArg(root.edge)) root.edge<-list(...)$root.edge
		else root.edge<-FALSE
		if(root.edge) ROOT<-if(!is.null(tree$root.edge)) tree$root.edge else 0
		else ROOT<-0 
		nHeight <- function(tree){
				tree <- reorder(tree)
				edge <- tree$edge
				el <- tree$edge.length
				res <- numeric(max(tree$edge))
				for(i in seq_len(nrow(edge))) res[edge[i,2]] <- res[edge[i,1]] + el[i] 
				res
		}
		nh <- nHeight(tree)
		return(matrix(nh[tree$edge], ncol=2L)+ROOT)
}

## function drops all the leaves from the tree & collapses singleton nodes
## written by Liam J. Revell 2013, 2014, 2015
drop.leaves<-function(tree,...){
	if(!inherits(tree,"phylo")) stop("tree should be an object of class \"phylo\".")
	## optional arguments
	if(hasArg(keep.tip.labels)) keep.tip.labels<-list(...)$keep.tip.labels
	else keep.tip.labels<-FALSE
	## end optional arguments
	n<-Ntip(tree)
	edge<-tree$edge
	edge[edge>n]<--edge[edge>n]+n
	ii<-which(edge[,2]>0)
	edge<-edge[-ii,]
	if(!is.null(tree$edge.length)){
		edge.length<-tree$edge.length
		edge.length<-edge.length[-ii]
	}
	zz<-sapply(edge[,2],function(x,y) !(x%in%y),y=edge[,1])
	if(is.null(tree$node.label)) tree$node.label<-1:tree$Nnode+n
	nn<-matrix(tree$node.label[-edge],nrow(edge),ncol(edge))
	tip.label<-nn[zz,2]
	node.label<-c(nn[1,1],nn[!zz,2])
	edge[zz,2]<-1:sum(zz)
	Nnode<-length(unique(edge[edge<0]))
	rr<-cbind(sort(unique(edge[edge<0]),decreasing=TRUE),1:Nnode+sum(zz))
	for(i in 1:nrow(rr)) edge[edge==rr[i,1]]<-rr[i,2]
	tt<-list(edge=edge,Nnode=Nnode,tip.label=tip.label,edge.length=edge.length,node.label=node.label)
	class(tt)<-"phylo"
	tt<-collapse.singles(tt)
	if(keep.tip.labels){
		for(i in 1:length(tt$tip.label)){
			yy<-getDescendants(tree,node=which(tree$node.label==tt$tip.label[i])+n)
			tt$tip.label[i]<-paste(tree$tip.label[yy[yy<=n]],collapse=",")
		}
	}
	return(tt)
}

# function rounds the branch lengths of the tree & applies rounding to simmap tree
# written by Liam J. Revell 2012, 2013, 2015
roundBranches<-function(tree,digits=0){
	if(inherits(tree,"multiPhylo")){
		trees<-lapply(tree,roundBranches,digits=digits)
		class(trees)<-"multiPhylo"
		return(trees)
	} else if(inherits(tree,"phylo")) {
		tree$edge.length<-round(tree$edge.length,digits)
		if(!is.null(tree$maps)){
			for(i in 1:nrow(tree$edge)){
				temp<-tree$maps[[i]]/sum(tree$maps[[i]])
				tree$maps[[i]]<-temp*tree$edge.length[i]
			}
		}
		if(!is.null(tree$mapped.edge)){
			a<-vector()
			for(i in 1:nrow(tree$edge)) a<-c(a,names(tree$maps[[i]]))
			a<-unique(a)
			tree$mapped.edge<-matrix(data=0,length(tree$edge.length),length(a),dimnames=list(apply(tree$edge,1,function(x) paste(x,collapse=",")),state=a))
			for(i in 1:length(tree$maps)) for(j in 1:length(tree$maps[[i]])) tree$mapped.edge[i,names(tree$maps[[i]])[j]]<-tree$mapped.edge[i,names(tree$maps[[i]])[j]]+tree$maps[[i]][j]
		}
		return(tree)
	} else stop("tree should be an object of class \"phylo\" or \"multiPhylo\".")
}

# function to merge mapped states
# written by Liam J. Revell 2013, 2015, 2019
mergeMappedStates<-function(tree,old.states,new.state){
	if(inherits(tree,"multiSimmap")){
		tree<-unclass(tree)
		tree<-lapply(tree,mergeMappedStates,old.states=old.states,new.state=new.state)
		class(tree)<-c("multiSimmap","multiPhylo")
	} else if(inherits(tree,"simmap")) {
		maps<-tree$maps
		rr<-function(map,oo,nn){ 
			for(i in 1:length(map)) if(names(map)[i]%in%oo) names(map)[i]<-nn
			map
		}
		mm<-function(map){
			if(length(map)>1){
				new.map<-vector()
				j<-1
				new.map[j]<-map[1]
				names(new.map)[j]<-names(map)[1]
				for(i in 2:length(map)){
					if(names(map)[i]==names(map)[i-1]){ 
						new.map[j]<-map[i]+new.map[j]
						names(new.map)[j]<-names(map)[i]
					} else {
						j<-j+1
						new.map[j]<-map[i]
						names(new.map)[j]<-names(map)[i]
					}

				}
				map<-new.map
			}
			map
		}
		maps<-lapply(maps,rr,oo=old.states,nn=new.state)
		if(length(old.states)>1){ 
			maps<-lapply(maps,mm)
			mapped.edge<-tree$mapped.edge
			mapped.edge<-cbind(rowSums(mapped.edge[,colnames(mapped.edge)%in%old.states]),
				mapped.edge[,setdiff(colnames(mapped.edge),old.states)])
			colnames(mapped.edge)<-c(new.state,setdiff(colnames(tree$mapped.edge),old.states))
		} else {
			mapped.edge<-tree$mapped.edge
			colnames(mapped.edge)[which(colnames(mapped.edge)==old.states)]<-new.state
		}
		tree$maps<-maps
		tree$mapped.edge<-mapped.edge
	} else stop("tree should be an object of class \"simmap\" or \"multiSimmap\".")
	return(tree)
}

# function rotates a node or multiple nodes
# written by Liam J. Revell 2013, 2015
rotateNodes<-function(tree,nodes,polytom=c(1,2),...){
	if(!inherits(tree,"phylo")) stop("tree should be an object of class \"phylo\".")
	n<-Ntip(tree)
	if(nodes[1]=="all") nodes<-1:tree$Nnode+n
	for(i in 1:length(nodes)) tree<-rotate(tree,nodes[i],polytom)
	if(hasArg(reversible)) reversible<-list(...)$reversible
	else reversible<-TRUE
	if(reversible){ 
		ii<-which(tree$edge[,2]<=n)
		jj<-tree$edge[ii,2]
		tree$edge[ii,2]<-1:n
		tree$tip.label<-tree$tip.label[jj]
	}
	return(tree)
}

# function simulates random sampling from xbar, xvar, with sample sizes n
# written by Liam J. Revell 2012
sampleFrom<-function(xbar=0,xvar=1,n=1,randn=NULL,type="norm"){
	if(length(xvar)==1&&length(xbar)!=length(xvar)) xvar<-rep(xvar,length(xbar))
	if(!is.null(randn))
		for(i in 1:length(xbar)) n[i]<-floor(runif(n=1,min=randn[1],max=(randn[2]+1)))
	x<-vector()

	for(i in 1:length(xbar)){
		y<-rnorm(n=n[i],mean=xbar[i],sd=sqrt(xvar[i]))
	 		names(y)<-rep(names(xbar)[i],length(y))
	 		x<-c(x,y)
	}
	return(x)
}

# function adds a new tip to the tree
# written by Liam J. Revell 2012, 2013, 2014, 2015
bind.tip<-function(tree,tip.label,edge.length=NULL,where=NULL,position=0,interactive=FALSE,...){
	if(!inherits(tree,"phylo")) stop("tree should be an object of class \"phylo\".")
	use.edge.length<-if(is.null(tree$edge.length)) FALSE else TRUE
	if(use.edge.length==FALSE) tree<-compute.brlen(tree)
	if(interactive==TRUE){
		plotTree(tree,...)
		cat(paste("Click where you would like to bind the tip \"",tip.label,"\"\n",sep=""))
		flush.console()
		obj<-get.treepos(message=FALSE)
		where<-obj$where
		position<-obj$pos
	} else if(is.null(where)) where<-Ntip(tree)+1
	if(where<=Ntip(tree)&&position==0){
		pp<-1e-12
		if(tree$edge.length[which(tree$edge[,2]==where)]<=1e-12){
			tree$edge.length[which(tree$edge[,2]==where)]<-2e-12
			ff<-TRUE
		} else ff<-FALSE
	} else pp<-position
	if(is.null(edge.length)&&is.ultrametric(tree)){
		H<-nodeHeights(tree)
		if(where==(Ntip(tree)+1)) edge.length<-max(H)
		else edge.length<-max(H)-H[tree$edge[,2]==where,2]+position
	}
	tip<-list(edge=matrix(c(2,1),1,2),
		tip.label=tip.label,
		edge.length=edge.length,
		Nnode=1)
		class(tip)<-"phylo"
	obj<-bind.tree(tree,tip,where=where,position=pp)
	if(where<=Ntip(tree)&&position==0){
		nn<-obj$edge[which(obj$edge[,2]==which(obj$tip.label==tip$tip.label)),1]
		obj$edge.length[which(obj$edge[,2]==nn)]<-obj$edge.length[which(obj$edge[,2]==nn)]+1e-12
		obj$edge.length[which(obj$edge[,2]==which(obj$tip.label==tip$tip.label))]<-0
		obj$edge.length[which(obj$edge[,2]==which(obj$tip.label==tree$tip.label[where]))]<-0
	}
	root.time<-if(!is.null(obj$root.time)) obj$root.time else NULL
	obj<-untangle(obj,"read.tree")
	if(!is.null(root.time)) obj$root.time<-root.time
	if(interactive) plotTree(obj,...)
	if(!use.edge.length) obj$edge.length<-NULL
	obj
}

## function collapses the subtree descended from node to a star tree
## written by Liam J. Revell 2013, 2015, 2019
collapse.to.star<-function(tree,node){
	if(!inherits(tree,"phylo")) stop("tree should be an object of class \"phylo\".")
	if(node==(Ntip(tree)+1)){
		object<-list(edge=cbind(rep(Ntip(tree)+1,Ntip(tree)),1:Ntip(tree)),
			tip.label=tree$tip.label,Nnode=1)
		if(!is.null(tree$edge.length)) object$edge.length<-sapply(1:Ntip(tree),nodeheight,tree=tree)
		class(object)<-"phylo"
		tree<-object
	} else {
		nel<-if(is.null(tree$edge.length)) TRUE else FALSE
		if(nel) tree$edge.length<-rep(1,nrow(tree$edge))
		tt<-splitTree(tree,split=list(node=node,
			bp=tree$edge.length[which(tree$edge[,2]==node)]))
		ss<-starTree(species=tt[[2]]$tip.label,
			branch.lengths=diag(vcv(tt[[2]])))
		ss$root.edge<-0
		tree<-paste.tree(tt[[1]],ss)
		if(nel) tree$edge.length<-NULL
	}
	tree
}

## function returns the MRCA, or its height above the root, for a set of taxa (in tips)
## written by Liam Revell 2012, 2013, 2015, 2016
findMRCA<-function(tree,tips=NULL,type=c("node","height")){
	type<-type[1]
	if(!inherits(tree,"phylo")) stop("tree should be an object of class \"phylo\".")
	if(is.null(tips)){ 
		X<-mrca(tree)
		if(type=="height"){
			H<-nodeHeights(tree)
			X<-apply(X,c(1,2),function(x,y,z) y[which(z==x)[1]],y=H,z=tree$edge)
		}
		return(X)
		} else {
		node<-getMRCA(tree,tips)
		if (type == "node") return(node)
		else if(type=="height") return(nodeheight(tree,node))
	}
}

# function works like extract.clade in ape but will preserve a discrete character mapping
# written by Liam J. Revell 2013
extract.clade.simmap<-function(tree,node){
	if(!inherits(tree,"phylo")) stop("tree should be an object of class \"phylo\".")
	x<-getDescendants(tree,node)
	x<-x[x<=Ntip(tree)]
	drop.tip.simmap(tree,tree$tip.label[-x])
}

# function gets all subtrees that cannot be further subdivided into two clades of >= clade.size
# written by Liam J. Revell 2013, 2015
getCladesofSize<-function(tree,clade.size=2){
	if(!inherits(tree,"phylo")) stop("tree should be an object of class \"phylo\".")
	n<-Ntip(tree)
	nn<-1:(tree$Nnode+n)
	ndesc<-function(tree,node){
		x<-getDescendants(tree,node)
		sum(x<=Ntip(tree))
	}
	dd<-setNames(sapply(nn,ndesc,tree=tree),nn)
	aa<-n+1 # root
	nodes<-vector()
	while(length(aa)){
		bb<-lapply(aa,function(x,tree) tree$edge[which(tree$edge[,1]==x),2],tree=tree)
		cc<-lapply(bb,function(x) dd[as.character(x)])
		gg<-sapply(cc,function(x,cs) any(x<cs),cs=clade.size)
		nodes<-c(nodes,aa[gg])
		aa<-unlist(bb[!gg])
	}
	trees<-lapply(nodes,extract.clade,phy=tree)
	class(trees)<-"multiPhylo"
	return(trees)
}

# function counts transitions from a mapped history
# written by Liam J. Revell 2013, 2015
countSimmap<-function(tree,states=NULL,message=TRUE){
	if(inherits(tree,"multiPhylo")){
		ff<-function(zz){
 			XX<-countSimmap(zz,states,message)
			setNames(c(XX$N,as.vector(t(XX$Tr))),c("N",
			sapply(rownames(XX$Tr),paste,colnames(XX$Tr),sep=",")))
		}
		tree<-unclass(tree)
		XX<-t(sapply(tree,ff))	
		if(!message) return(XX)
		else return(list(Tr=XX,message=
			c("Column N is the total number of character changes on the tree",
			"Other columns give transitions x,y from x->y")))
	} else if(inherits(tree,"phylo")) {
		n<-sum(sapply(tree$maps,length))-nrow(tree$edge)
		if(is.null(states)) states<-colnames(tree$mapped.edge)
		m<-length(states)	
		TT<-matrix(NA,m,m,dimnames=list(states,states))
		gg<-function(map,a,b){
			if(length(map)==1) zz<-0
			else {
				zz<-0; i<-2
				while(i<=length(map)){
					if(names(map)[i]==b&&names(map)[i-1]==a) zz<-zz+1
				i<-i+1
				}
			} 
			return(zz)
		}
		for(i in 1:m) for(j in 1:m)
			if(i==j) TT[i,j]<-0
			else TT[i,j]<-sum(sapply(tree$maps,gg,a=states[i],b=states[j]))
		if(!message) return(list(N=n,Tr=TT))
		else return(list(N=n,Tr=TT,message=c(
			"N is the total number of character changes on the tree",
			"Tr gives the number of transitions from row state->column state")))
	}
}

# function to match nodes between trees
# written by Liam J. Revell 2012, 2013, 2015
matchNodes<-function(tr1,tr2,method=c("descendants","distances"),...){
	if(!inherits(tr1,"phylo")||!inherits(tr1,"phylo")) stop("tr1 & tr2 should both be objects of class \"phylo\".")
	if(hasArg(quiet)) quiet<-list(...)$quiet
	else quiet<-FALSE
	method<-method[1]
	method<-matchType(method,c("descendants","distances"))
	if(method=="descendants"){
		desc.tr1<-lapply(1:tr1$Nnode+length(tr1$tip),function(x) extract.clade(tr1,x)$tip.label)
		names(desc.tr1)<-1:tr1$Nnode+length(tr1$tip)
		desc.tr2<-lapply(1:tr2$Nnode+length(tr2$tip),function(x) extract.clade(tr2,x)$tip.label)
		names(desc.tr2)<-1:tr2$Nnode+length(tr2$tip)
		Nodes<-matrix(NA,length(desc.tr1),2,dimnames=list(NULL,c("tr1","tr2")))
		for(i in 1:length(desc.tr1)){
			Nodes[i,1]<-as.numeric(names(desc.tr1)[i])
			for(j in 1:length(desc.tr2))
				if(all(desc.tr1[[i]]%in%desc.tr2[[j]])&&all(desc.tr2[[j]]%in%desc.tr1[[i]]))
					Nodes[i,2]<-as.numeric(names(desc.tr2)[j])
		}
	} else if(method=="distances"){
		if(hasArg(tol)) tol<-list(...)$tol
		else tol<-1e-6
		if(hasArg(corr)) corr<-list(...)$corr
		else corr<-FALSE
		if(corr) tr1$edge.length<-tr1$edge.length/max(nodeHeights(tr1))
		if(corr) tr2$edge.length<-tr2$edge.length/max(nodeHeights(tr2))
		D1<-dist.nodes(tr1)[1:length(tr1$tip),1:tr1$Nnode+length(tr1$tip)]
		D2<-dist.nodes(tr2)[1:length(tr2$tip),1:tr2$Nnode+length(tr2$tip)]
		rownames(D1)<-tr1$tip.label
		rownames(D2)<-tr2$tip.label
		common.tips<-intersect(tr1$tip.label,tr2$tip.label)
		D1<-D1[common.tips,]
		D2<-D2[common.tips,]
		Nodes<-matrix(NA,tr1$Nnode,2,dimnames=list(NULL,c("tr1","tr2")))
		for(i in 1:tr1$Nnode){
			if(corr) z<-apply(D2,2,function(X,y) cor(X,y),y=D1[,i])
			else z<-apply(D2,2,function(X,y) 1-sum(abs(X-y)),y=D1[,i])
			Nodes[i,1]<-as.numeric(colnames(D1)[i])
			if(any(z>=(1-tol))){
				a<-as.numeric(names(which(z>=(1-tol))))
				if(length(a)==1) Nodes[i,2]<-a
				else {
					Nodes[i,2]<-a[1]
					if(!quiet) warning("polytomy detected; some node matches may be arbitrary")
				}
			}
		}
	}
	return(Nodes)
}

# function applies the branch lengths of a reference tree to a second tree, including mappings

# written by Liam J. Revell 2012, 2015
applyBranchLengths<-function(tree,edge.length){
	if(inherits(tree,"multiPhylo")){
		trees<-lapply(tree,applyBranchLengths,edge.length=edge.length)
		class(trees)<-"multiPhylo"
		return(trees)
	} else if(inherits(tree,"phylo")) {
		tree$edge.length<-edge.length
		if(!is.null(tree$maps)){
			for(i in 1:nrow(tree$edge)){
				temp<-tree$maps[[i]]/sum(tree$maps[[i]])
				tree$maps[[i]]<-temp*tree$edge.length[i]
			}
		}
		if(!is.null(tree$mapped.edge)){
			a<-vector()
			for(i in 1:nrow(tree$edge)) a<-c(a,names(tree$maps[[i]]))
			a<-unique(a)
			tree$mapped.edge<-matrix(data=0,length(tree$edge.length),length(a),dimnames=list(apply(tree$edge,1,function(x) paste(x,collapse=",")),state=a))
			for(i in 1:length(tree$maps)) for(j in 1:length(tree$maps[[i]])) tree$mapped.edge[i,names(tree$maps[[i]])[j]]<-tree$mapped.edge[i,names(tree$maps[[i]])[j]]+tree$maps[[i]][j]
		}
		return(tree)
	}
}

# function to compute phylogenetic VCV using joint Pagel's lambda
# written by Liam Revell 2011, 2024

phyl.vcv<-function(X,C,lambda){
	if(!is.null(rownames(X))) C<-C[rownames(X),rownames(X)] ## sort by rownames of X (if present)
	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 for joint estimation of lambda for multiple traits
# written by Liam Revell 2011/2012
likMlambda<-function(lambda,X,C){
	# compute R, conditioned on lambda
	temp<-phyl.vcv(X,C,lambda);
	C<-temp$C; R<-temp$R; a<-temp$alpha
	# 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
	kronRC<-kronecker(R,C)
	logL<--t(y-D%*%a)%*%solve(kronRC,y-D%*%a)/2-n*m*log(2*pi)/2-determinant(kronRC,logarithm=TRUE)$modulus/2
	return(logL)
}

# function matches data to tree
# written by Liam J. Revell 2011, 2015
matchDatatoTree<-function(tree,x,name){
	if(!inherits(tree,"phylo")) stop("tree should be an object of class \"phylo\".")
	if(is.matrix(x)) x<-x[,1]
	if(is.null(names(x))){
		if(length(x)==Ntip(tree)){
			print(paste(name,"has no names; assuming x is in the same order as tree$tip.label"))
			names(x)<-tree$tip.label
		} else

			stop(paste(name,"has no names and is a different length than tree$tip.label"))

	}
	if(any(is.na(match(names(x),tree$tip.label)))){
		print(paste("some species in",name,"are missing from tree, dropping missing taxa from",name))
		x<-x[intersect(tree$tip.label,names(x))]
	}
	return(x)
}

# function matches tree to data
# written by Liam J. Revell 2011, 2015
matchTreetoData<-function(tree,x,name){
	if(!inherits(tree,"phylo")) stop("tree should be an object of class \"phylo\".")
	if(any(is.na(match(tree$tip.label,names(x))))){
		print(paste("some species in tree are missing from",name,", dropping missing taxa from the tree"))
		tree<-drop.tip(tree,setdiff(tree$tip.label,names(x)))
	}
	if(any(is.na(x))){
		print(paste("some data in",name,"given as 'NA', dropping corresponding species from tree"))
		tree<-drop.tip(tree,names(which(is.na(x))))
	}
	return(tree)
}

# function finds the maximum value of Pagel's lambda
# written by Liam J. Revell 2011
maxLambda<-function(tree){
	if(!inherits(tree,"phylo")) stop("tree should be an object of class \"phylo\".")
	if(is.ultrametric(tree)){
		H<-nodeHeights(tree)
		return(max(H[,2])/max(H[,1]))
	} else return(1)
}

# function reorders the columns of mapped.edge from a set of simmap trees

# written by Liam J. Revell 2013, 2015
orderMappedEdge<-function(trees,ordering=NULL){
	if(!inherits(trees,"phylo")&&!inherits(trees,"multiPhylo")) 
		stop("trees should be an object of class \"phylo\" or \"multiPhylo\".")
	f1<-function(tree,ordering){
		mapped.edge<-matrix(0,nrow(tree$mapped.edge),length(ordering),
			dimnames=list(rownames(tree$mapped.edge),ordering))
		mapped.edge[,colnames(tree$mapped.edge)]<-tree$mapped.edge
		tree$mapped.edge<-mapped.edge
		return(tree)
	}
	f2<-function(tree) colnames(tree$mapped.edge)
	if(inherits(trees,"phylo")) states<-colnames(trees$mapped.edge)
	else if(inherits(trees,"multiPhylo")) states<-unique(as.vector(sapply(trees,f2)))
	if(length(ordering)>1) if(length(intersect(states,ordering))<length(states)){
		warning("not all states represented in input ordering; setting to default")
		ordering<-NULL
	}
	if(is.null(ordering)) ordering<-"alphabetical"
	if(length(ordering)==1){
		ordering<-matchType(ordering,c("alphabetical","numerical"))
		if(ordering=="alphabetical") ordering<-sort(states)
		else if(ordering=="numerical") ordering<-as.character(sort(as.numeric(states)))
	}
	if(inherits(trees,"phylo")) trees<-f1(trees,ordering)
	else { 
		trees<-lapply(trees,f1,ordering=ordering)
		class(trees)<-"multiPhylo"
	}
	return(trees)
}

# function gets sister node numbers or names
# written by Liam J. Revell 2013, 2015
getSisters<-function(tree,node,mode=c("number","label")){
	if(!inherits(tree,"phylo")) stop("tree should be an object of class \"phylo\".")
	mode<-mode[1]
	if(is.character(node)) node<-match(node,c(tree$tip.label,tree$node.label))
	sisters<-tree$edge[which(tree$edge[,1]==tree$edge[which(tree$edge[,2]==node),1]),2]
	sisters<-setdiff(sisters,node)
	if(mode=="number") return(sisters)
	else if(mode=="label"){
		result<-list()
		n<-Ntip(tree)
		if(is.null(tree$node.label)&&any(sisters>n)) result$nodes<-sisters[which(sisters>n)] 
		else if(any(sisters>n)) result$nodes<-tree$node.label[sisters[which(sisters>n)]-n]
		if(any(sisters<=n)) result$tips<-tree$tip.label[sisters[which(sisters<=n)]]
		return(result)
	}
}

# gets descendant node numbers
# written by Liam Revell 2012, 2013, 2014
getDescendants<-function(tree,node,curr=NULL){
	if(!inherits(tree,"phylo")) stop("tree should be an object of class \"phylo\".")
	if(is.null(curr)) curr<-vector()
	daughters<-tree$edge[which(tree$edge[,1]==node),2]
	curr<-c(curr,daughters)
	if(length(curr)==0&&node<=Ntip(tree)) curr<-node
	w<-which(daughters>Ntip(tree))
	if(length(w)>0) for(i in 1:length(w)) 
		curr<-getDescendants(tree,daughters[w[i]],curr)
	return(curr)
}

# function computes vcv for each state, and stores in array
# written by Liam J. Revell 2011, 2012, 2016
multiC<-function(tree,internal=FALSE){
	if(!inherits(tree,"phylo")) stop("tree should be an object of class \"phylo\".")
	if(!inherits(tree,"simmap")) stop("tree should be an object of class \"simmap\".")
	m<-ncol(tree$mapped.edge)
	# compute separate C for each state
	mC<-list()
	for(i in 1:m){
		mtree<-list(edge=tree$edge,
			Nnode=tree$Nnode,
			tip.label=tree$tip.label,
			edge.length=tree$mapped.edge[,i])
		class(mtree)<-"phylo"
		mC[[i]]<-if(internal) vcvPhylo(mtree,internal=TRUE) else vcv.phylo(mtree)
	}
	names(mC)<-colnames(tree$mapped.edge)
	mC
}

# function pastes subtree onto tip
# written by Liam Revell 2011, 2015
paste.tree<-function(tr1,tr2){
	if(!inherits(tr1,"phylo")||!inherits(tr2,"phylo")) stop("tr1 & tr2 should be objects of class \"phylo\".")
	if(length(tr2$tip)>1){ 
		temp<-tr2$root.edge; tr2$root.edge<-NULL
		tr1$edge.length[match(which(tr1$tip.label=="NA"),tr1$edge[,2])]<-tr1$edge.length[match(which(tr1$tip.label=="NA"),tr1$edge[,2])]+temp
	}
	tr.bound<-bind.tree(tr1,tr2,where=which(tr1$tip.label=="NA"))	
	return(tr.bound)
}

# match type
# written by Liam J. Revell 2012
matchType<-function(type,types){
	for(i in 1:length(types))
		if(all(strsplit(type,split="")[[1]]==strsplit(types[i],split="")[[1]][1:length(strsplit(type,split="")[[1]])]))
			type=types[i]
	return(type)
}
	
## function 'untangles' (or attempts to untangle) a tree with crossing branches
## written by Liam J. Revell 2013, 2015, 2020, 2021
untangle<-function(tree,method=c("reorder","read.tree")){
	if(inherits(tree,"multiPhylo")){
		tree<-lapply(tree,untangle,method=method)
		class(tree)<-"multiPhylo"
	} else {
		if(!inherits(tree,"phylo")) stop("tree should be an object of class \"phylo\".")
		obj<-attributes(tree)
		method<-method[1]
		if(method=="reorder") tree<-reorder(reorder(tree,"pruningwise"))
		else if(method=="read.tree"){
			tip.label<-tree$tip.label
			tree$tip.label<-1:Ntip(tree)
			if(inherits(tree,"simmap")) 
				tree<-read.simmap(text=capture.output(write.simmap(tree,file=stdout())))
			else tree<-if(Ntip(tree)>1) read.tree(text=write.tree(tree)) else read.newick(text=write.tree(tree))
			tree$tip.label<-tip.label[as.numeric(tree$tip.label)]
		}
		ii<-!names(obj)%in%names(attributes(tree))
		attributes(tree)<-c(attributes(tree),obj[ii])
	}
	tree
}
liamrevell/phytools documentation built on March 4, 2024, 3:27 a.m.