# function performs least-squares phylogeny inference by nni # written by Liam J. Revell 2011 optim.phylo.ls<-function(D,stree=NULL,set.neg.to.zero=TRUE,fixed=FALSE,tol=1e-10,collapse=TRUE){ # load phangorn package if(!require(phangorn)) stop("must first install \"phangorn\" package.") # change D to a matrix (if actually an object of class "dist") if(class(D)=="dist") D<-as.matrix(D) # compute the number of species n<-nrow(D) if(is.null(stree)) stree<-rtree(n=n,tip.label=rownames(D),br=NULL,rooted=F) # random starting tree else if(class(stree)!="phylo"){ message("starting tree must be an object of class \"phylo.\" using random starting tree.") stree<-rtree(n=n,tip.label=rownames(D),br=NULL,rooted=F) # random starting tree } if(!is.binary.tree(stree)) stree<-multi2di(stree) if(is.rooted(stree)) stree<-unroot(stree) # get ls branch lengths for stree best.tree<-ls.tree(stree,D) Q<-attr(best.tree,"Q-score") bestQ<-0 # to start the loop # for search Nnni<-0 # loop while Q is not improved by nni while(bestQ-Qi) { rnames[k]<-paste(c(i,j),collapse=","); k<-k+1 } X<-matrix(0,n*(n-1)/2,nrow(tree$edge),dimnames=list(rnames,cnames)) anc.nodes<-compute.ancestor.nodes(tree) anc.nodes<-anc.nodes[,c(n+1:tree$Nnode,1:n)] for(i in 1:(n-1)){ for(j in (i+1):n){ nodes.i<-names(anc.nodes[i,anc.nodes[i,]==1]) for(k in 1:length(nodes.i)) X[paste(c(i,j),collapse=","),match(as.numeric(nodes.i[k]),tree$edge[,2])]<-X[paste(c(i,j),collapse=","),match(as.numeric(nodes.i[k]),tree$edge[,2])]+1 nodes.j<-names(anc.nodes[j,anc.nodes[j,]==1]) for(k in 1:length(nodes.j)) X[paste(c(i,j),collapse=","),match(as.numeric(nodes.j[k]),tree$edge[,2])]<-X[paste(c(i,j),collapse=","),match(as.numeric(nodes.j[k]),tree$edge[,2])]+1 } } X[X>1]<-0 return(X) } # function computes the ancestor node numbers for each tip number # written by Liam J. Revell 2011 compute.ancestor.nodes<-function(tree){ n<-length(tree$tip) m<-tree$Nnode X<-matrix(0,n,n+m,dimnames=list(1:n,1:(n+m))) for(i in 1:n){ currnode<-i while(currnode!=(n+1)){ X[i,currnode]<-1 currnode<-tree$edge[match(currnode,tree$edge[,2]),1] } X[i,currnode]<-1 } return(X) }