R/exhaustiveMP.R

In phytools: Phylogenetic Tools for Comparative Biology (and Other Things)

# function does branch & bound or exhaustive MP tree search
# uses "phangorn (Schliep 2011) & "ape" (Paradis & Strimmer 2004)
# data is a phyDat object; method can be "branch.and.bound" or "exhaustive"
# written by Liam J. Revell 2011, 2013

exhaustiveMP<-function(data,tree=NULL,method="branch.and.bound"){
	if(method=="branch.and.bound"){
		if(length(data)>15) stop("branch and bound only allowed for n<=15")
		if(is.null(tree)){
			if(attr(data,"type")=="DNA"){
				print("no input tree; starting with NJ tree")
				tree<-NJ(dist.dna(as.DNAbin(data)))
			} else {
				print("no input tree; using random starting tree")
				tree<-rtree(n=length(data),tip.label=names(data),br=NULL,rooted=FALSE)
			}
		}
		trees<-branch.and.bound(data,tree)
	} else if(method=="exhaustive"){
		if(length(data)>10) stop("exhaustive search only allowed for n<=10")
		if(!is.null(tree)) print("starting tree not necessary for exhaustive search")
		trees<-exhaustive.search(data)
	}
	if(length(trees)==1) trees<-trees[[1]]
	return(trees)
}

# function performs branch & bound search
# written by Liam J. Revell 2011

branch.and.bound<-function(data,tree){
	# first, compute the parsimony score on the input tree
	bound<-parsimony(tree,data)
	# now pick three species to start
################## this is new and I don't know if it helps
	if(is.null(tree$edge.length)){
		print("starting with 3 species chosen at random")
		new<-list(stree(n=3,tip.label=sample(tree$tip.label,3)))
	} else {
		print("starting with 3 species chosen to maximize distance") 
		mdSp<-names(sort(colSums(cophenetic(tree)),decreasing=TRUE))[1:3]
		mdSp<-c("Tarsier","Chimp","J.Macaque")
		print(mdSp)
		new<-list(stree(n=3,tip.label=mdSp)) 
	}
################## ends here
	class(new)<-"multiPhylo"
	added<-new[[1]]$tip.label; remaining<-setdiff(tree$tip.label,added)
	# branch & bound
	while(length(remaining)>0){
		old<-new; new<-list()
		new.tip<-sample(remaining,1)
		pscores<-vector()
		for(i in 1:length(old)){			
			temp<-add.everywhere(old[[i]],new.tip)
			score<-parsimony(temp,data)
			new<-c(new,temp[score<=bound])
			pscores<-c(pscores,score[score<=bound])
		}
		added<-c(added,new.tip)
		print(paste(length(added),"species added;",length(new),"trees retained",collapse=""))
		remaining<-setdiff(tree$tip.label,added)
	}
	# ok, done, now sort what needs to be returned
	trees<-new[pscores==min(pscores)]
	for(i in 1:length(trees)) attr(trees[[i]],"pscore")<-min(pscores)
	return(trees) # return all mp trees
}

# function does exhaustive tree search
# written by Liam J. Revell 2011

exhaustive.search<-function(data){
	all.trees<-allTrees(n=length(data),tip.label=names(data),rooted=FALSE)
	print(paste("searching",length(all.trees),"trees",collapse=""))
	all.trees = .uncompressTipLabel(all.trees)
	pscores<-parsimony(all.trees,data)
	minscore<-min(pscores); trees<-all.trees[pscores==minscore]
	for(i in 1:length(trees)) attr(trees[[i]],"pscore")<-min(pscores)
	return(trees)
}