PageRenderTime 50ms CodeModel.GetById 17ms app.highlight 26ms RepoModel.GetById 1ms app.codeStats 1ms

/in.progress/auteurExtended/R/phylogenetic.utilities.R

http://github.com/eastman/auteur
R | 215 lines | 157 code | 39 blank | 19 comment | 19 complexity | a7ace8c4c1277ec7bbbc2ce7cadd2c70 MD5 | raw file
  1oumat<-function (vmat, alpha, sigma) 
  2{
  3	tips=unique(dim(vmat))
  4	
  5	out <- .Call("oumat",	
  6				 TIMES=list(vmat=as.double(array(vmat))),
  7				 ALPHA=as.double(alpha^2),
  8				 SIGMA=as.double(sigma), 
  9				 TIPS=as.integer(tips),
 10				 PACKAGE = "auteurExtended")
 11	return(matrix(out$oumat,ncol=tips))
 12}
 13
 14
 15
 16ouweights <- function(phy, epochs, object, ancestors, alpha){
 17	if(!all(c("anc","des","time","regimes")%in%names(object)))stop("Cannot decipher supplied 'object'.")
 18	ntip = Ntip(phy)
 19	object=object[order(object$des),]
 20	optima = sort(unique(as.numeric(object$regimes)))
 21    nreg = length(optima)
 22	
 23	out <- .Call("ouweights",	
 24				 PHYLO=list(
 25							epochs=epochs,
 26							ancestors=ancestors,
 27							regimes=as.double(object$regimes),
 28							optima=as.double(optima),
 29							alpha=as.double(alpha^2),
 30							tips=as.integer(ntip)),
 31				 PACKAGE = "auteurExtended")
 32	
 33	W=matrix(out$weights,ncol=nreg)
 34	opt=optima
 35	return(list(W=W, optima=opt))
 36} 
 37
 38vmat <- function(phy){
 39	n=Ntip(phy)
 40	out <- .Call("vmat", tree=list(
 41									ROOT = as.integer(n+1),
 42									MAXNODE = as.integer(max(phy$edge[,1])),
 43									ENDOFCLADE = as.integer(dim(phy$edge)[1]),
 44									ANC = as.integer(phy$edge[,1]),
 45									DES = as.integer(phy$edge[,2]),
 46									EDGES = as.double(c(phy$edge.length,0)),
 47									VCV = as.double(array(matrix(0, n, n)))),
 48									PACKAGE = "auteurExtended")
 49	v=matrix(out$VCV,nrow=n,byrow=FALSE)
 50	rownames(v)<-colnames(v)<-phy$tip.label
 51	return(v)
 52} 
 53
 54		
 55
 56ultrametricize=function(phy, tol=1e-8, trim=c("min","max","mean")){
 57	paths=diag(vmat(phy))
 58	trim=switch(match.arg(trim),
 59				min = min(paths),
 60				max = max(paths),
 61				mean = mean(paths))
 62	if(diff(range(paths))<=tol) {
 63		for(i in 1:length(paths)){
 64			e=phy$edge.length[which(phy$edge[,2]==i)->bl]
 65			phy$edge.length[bl]=e+(trim-paths[i])
 66		}
 67		return(phy)
 68	} else {
 69		stop("Difference in path lengths is larger than supplied tolerance")
 70	}
 71}
 72
 73#general phylogenetic utility for finding the edge labels (from phy$edge[,2]) associated with the most exclusive monophyletic grouping containing at least supplied 'tips' where 'tips' are elements of the phy$tip.label vector
 74#author: JM EASTMAN 2010
 75
 76collect.nodes <-
 77function(tips,phy,strict=FALSE){
 78	if(length(tips)==1) stop("Must supply at least two tips.")
 79	tt=phy$tip.label
 80	tt=lapply(tips, function(x) which(phy$tip.label==x))
 81	nodes=lapply(tt, function(x) get.ancestors.of.node(x, phy))
 82	all.nodes=nodes[[1]]
 83	for(i in 2:length(nodes)) {
 84		all.nodes=intersect(all.nodes,nodes[[i]])
 85	}
 86	base.node=max(all.nodes)
 87	if(strict) {
 88		u=unlist(nodes)
 89		return(c(sort(unique(u[u>=base.node])), unlist(tt)))
 90	} else {
 91		return(c(base.node,get.descendants.of.node(base.node,phy)))
 92	}
 93}
 94
 95#general phylogenetic utility for finding the phy$tip.labels associated with the most exclusive monophyletic grouping containing at least supplied 'tips' where 'tips' are elements of the phy$tip.label vector
 96#author: JM EASTMAN 2010
 97
 98find.relatives <-
 99function(tips,phy){
100	nn=collect.nodes(tips, phy, strict=FALSE)
101	return(phy$tip.label[nn[nn<=Ntip(phy)]])
102}
103
104#general phylogenetic utility for returning first ancestor (as a numeric referent) of the supplied node 
105#author: JM EASTMAN 2010
106
107get.ancestor.of.node <-
108function(node, phy) {
109	return(phy$edge[which(phy$edge[,2]==node),1])
110}
111
112
113
114#general phylogenetic utility for returning all ancestors (listed as given in phy$edge[,2]) of a node 
115#author: JM EASTMAN 2010
116
117get.ancestors.of.node <-
118function(node, phy) {
119	a=c()
120	if(node==(Ntip(phy)+1->root)) return(NULL)
121	f=get.ancestor.of.node(node, phy)
122	a=c(a,f)
123	if(f>root) a=c(a, get.ancestors.of.node(f, phy))
124	return(a)
125}
126
127
128#general phylogenetic utility for returning the first (usually, unless a polytomy exists) two descendants the supplied node 
129#author: JM EASTMAN 2010
130
131get.desc.of.node <-
132function(node, phy) {
133	return(phy$edge[which(phy$edge[,1]==node),2])
134}
135
136
137#general phylogenetic utility for returning all descendants (listed as given in phy$edge[,2]) of a node (excluding the supplied node) 
138#author: JM EASTMAN 2010, based on GEIGER:::node.leaves()
139
140get.descendants.of.node <-
141function (node, phy, tips=FALSE) 
142{
143    node <- as.numeric(node)
144    n <- Ntip(phy)
145	if(node <= n) return(NULL)
146    l <- c()
147    d <- get.desc.of.node(node, phy)
148    for (j in d) {
149		l <- c(l, j)
150        if (j > n) {
151			l <- c(l, get.descendants.of.node(j, phy))
152		}
153    }
154    if(tips) return(l[l<=n]) else return(l)
155}
156
157
158get.epochs<-function(phy) 
159# ordered: tip to root
160{
161	ancestors=lapply(1:Ntip(phy), function(x) c(x,rev(sort(get.ancestors.of.node(x,phy)))))
162	phylo=get.times(phy)
163	tt=phylo$time[order(phylo$des)]
164	e=lapply(ancestors, function(x) tt[x])
165	names(e)=phy$tip.label
166	return(e)
167}
168
169get.times<-function(phy) 
170{
171	df=data.frame(phy$edge)
172	names(df)=c("anc","des")
173	anc=lapply(df$des,get.ancestors.of.node,phy)
174	order.edges=c(0,phy$edge.length)[order(c(Ntip(phy)+1,df$des))]
175	df$time=0
176	for(n in 1:length(anc)){
177		df$time[n]=sum(order.edges[anc[[n]]])+order.edges[df$des[n]]
178	}
179	
180	df=rbind(c(0,Ntip(phy)+1,0),df)
181	
182	return(df)
183}
184
185#general phylogenetic utility for determining whether a node is the root of the phylogeny
186#author: JM EASTMAN 2010
187
188is.root <-
189function(node,phy) {
190	if(node==phy$edge[1,1]) return(TRUE) else return(FALSE)
191}
192
193#general phylogenetic utility for returning all ancestors of nodes in a phylogeny
194#author: JM EASTMAN 2011
195
196get.ancestors <-
197function(phy) {
198	r=lapply(1:Ntip(phy), function(x) c(x,rev(sort(get.ancestors.of.node(x,phy)))))
199	r
200}
201
202
203#general phylogenetic utility for pruning the most recent speciation event from a phylogeny
204#author: JM EASTMAN 2010
205
206prunelastsplit<-function(phy) {
207	require(ape)
208	nn=as.numeric(names(xx<-branching.times(phy))[which(xx==min(xx))])
209	for(node in 1:length(nn)) {	
210		sp=sample(phy$edge[which(phy$edge[,1]==nn[node]),2],1)
211		obj=drop.tip(phy,phy$tip.label[sp])
212	}
213	return(obj)
214}
215