/archive/source.archive/bm.and.ou/trait.rjmcmc.10042010.R

  1## Written by LJ HARMON & AH: 2008
  2## Modified by JM EASTMAN: 09.2010
  3
  4rjMCMC.trait<-function (phy, data, ngen=1000, sampleFreq=100, model=c("OU","BM"), probAlpha=0.05, probMergeSplit=0.05, probRoot=0.1, theta.rates.ratio=1, heat=1, fileBase="result") 
  5{
  6	if(length(model)>1)stop(paste("Please specify either ", sQuote("OU"), " or ", sQuote("BM")))
  7## MORE ERROR CHECKING NEEDED HERE ## 
  8	require(ouch)
  9	require(geiger)
 10	
 11### prepare objects for rjMCMC
 12	cur.model		<- model								
 13
 14	dataList		<- prepare.ouchdata(phy, data)			
 15	ape.tre			<- cur.tre <- dataList$ape.tre			
 16	orig.dat		<- dataList$orig.dat						
 17	ouch.tre		<- dataList$ouch.tre					
 18	ouch.dat		<- dataList$ouch.dat						
 19	nn				<- length(ape.tre$edge.length)	
 20	node.des		<- sapply(unique(c(ape.tre$edge[1,1],ape.tre$edge[,2])), function(x) get.descendants.of.node(x, ape.tre))
 21	names(node.des) <- c(ape.tre$edge[1,1], unique(ape.tre$edge[,2]))
 22	ll <- list()										
 23
 24# initialize parameters
 25	init.rate		<- generate.starting.point(data,ape.tre,node.des,theta=FALSE)
 26    cur.rates		<- init.rate$values			
 27	cur.delta.rates	<- init.rate$delta
 28	cur.root		<- adjust.value(mean(orig.dat))
 29	cur.alpha		<- adjust.rate(0.001)
 30	init.theta		<- generate.starting.point(data,ape.tre,node.des,theta=TRUE)
 31	cur.theta		<- init.theta$values							
 32	cur.delta.theta	<- init.theta$delta	
 33	cur.root.theta	<- assign.root.theta(cur.theta,ape.tre)
 34	cur.regimes		<- as.factor(c(cur.root.theta,cur.theta)) 
 35	cur.tre			<- apply.BMM.to.apetree(ape.tre, cur.rates)
 36	if(model=="OU")	{
 37		mod.cur = new.ou.lik.fn(cur.rates, cur.alpha, cur.regimes, cur.theta, cur.tre, ouch.dat)
 38	} else {
 39		mod.cur = new.bm.lik.fn(cur.rates, cur.root, cur.tre, orig.dat)
 40	}
 41
 42
 43# proposal counts
 44	nRateProp =		0	
 45	nRateSwapProp = 0									
 46	nRcatProp =		0										
 47	nRootProp =		0
 48	nAlphaProp =	0										
 49	nTcatProp =		0										
 50	nThetaProp =	0
 51	nThetaSwapProp= 0
 52	nRateOK =		0											
 53	nRateSwapOK =	0
 54	nRcatOK =		0											
 55	nRootOK =		0
 56	nAlphaOK =		0										
 57	nTcatOK =		0
 58	nThetaOK =		0
 59	nThetaSwapOK =	0
 60
 61	cOK=c(												
 62		nRateOK,
 63		nRateSwapOK,
 64		nRcatOK, 
 65		nRootOK, 
 66		nAlphaOK,
 67		nThetaOK,
 68		nThetaSwapOK,
 69		nTcatOK 
 70	)		
 71	
 72	tickerFreq=ceiling(ngen/30)
 73	
 74# file handling
 75	errorLog=file(paste(getwd(),paste(cur.model, fileBase, "rjTraitErrors.log",sep="."),sep="/"),open='w+')
 76	generate.error.message(i=NULL, mod.cur=NULL, mod.new=NULL, lnR=NULL, errorLog=errorLog, init=TRUE)
 77	runLog=file(paste(getwd(),paste(cur.model, fileBase, "rjTrait.log",sep="."),sep="/"),open='w+')
 78	generate.log(bundled.parms=NULL, cur.model, file=runLog, init=TRUE)
 79
 80### Begin rjMCMC
 81    for (i in 1:ngen) {
 82        lnLikelihoodRatio <- lnHastingsRatio <- lnPriorRatio <- 0
 83		startparms = c(nRateProp, nRateSwapProp, nRcatProp, nRootProp, nAlphaProp, nThetaProp, nThetaSwapProp, nTcatProp)
 84	
 85## OU IMPLEMENTATION ## 
 86		if(cur.model=="OU") {
 87			if (runif(1) < (2 * probMergeSplit)) {
 88				if(runif(1)<(theta.rates.ratio/(theta.rates.ratio+1))) {	# adjust theta categories
 89					nr=split.or.merge(cur.delta.theta, cur.theta, ape.tre, node.des, theta=TRUE)
 90					new.alpha=cur.alpha
 91					new.theta=nr$new.values ##
 92					new.delta.theta=nr$new.delta ##
 93					new.root.theta<-nrt<-assign.root.theta(new.theta, ape.tre) ## 
 94					new.regimes=as.factor(c(nrt,new.theta)) ##
 95					new.rates=cur.rates 
 96					new.delta.rates=cur.delta.rates 
 97					nTcatProp=nTcatProp+1
 98				} else {							# adjust rate categories
 99					nr=split.or.merge(cur.delta.rates, cur.rates, ape.tre, node.des, theta=FALSE)
100					new.alpha=cur.alpha
101					new.root.theta=cur.root.theta
102					new.theta=cur.theta
103					new.delta.theta=cur.delta.theta
104					new.regimes=cur.regimes
105					new.rates=nr$new.values ##
106					new.delta.rates=nr$new.delta ##
107					nRcatProp=nRcatProp+1
108				}
109				lnHastingsRatio=nr$lnHastingsRatio
110				lnPriorRatio=nr$lnPriorRatio
111			} else { # neither split nor merge
112				if(runif(1)<probAlpha) {			# adjust alpha
113					new.alpha=adjust.rate(cur.alpha) ##
114					new.root.theta=cur.root.theta
115					new.theta=cur.theta
116					new.delta.theta=cur.delta.theta
117					new.regimes=cur.regimes 
118					new.rates=cur.rates
119					new.delta.rates=cur.delta.rates
120					nAlphaProp=nAlphaProp+1
121				} else { 
122					if(runif(1)>(theta.rates.ratio/(theta.rates.ratio+1))) {
123						if(runif(1)<0.2 | sum(cur.delta.rates)==0) {
124							new.alpha=cur.alpha		# adjust rates
125							new.root.theta=cur.root.theta
126							new.theta=cur.theta
127							new.delta.theta=cur.delta.theta
128							new.regimes=cur.regimes 
129							new.rates=adjust.rate(cur.rates) ##
130							new.delta.rates=cur.delta.rates
131							nRateProp=nRateProp+1
132						} else {					# swap rates
133							nr=switcheroo(cur.delta.rates, cur.rates, ape.tre, node.des, theta=FALSE)
134							new.alpha=cur.alpha		
135							new.root.theta=cur.root.theta
136							new.theta=cur.theta
137							new.delta.theta=cur.delta.theta
138							new.regimes=cur.regimes 
139							new.rates=nr$new.values ##
140							new.delta.rates=nr$new.delta ##
141							nRateSwapProp=nRateSwapProp+1
142						}
143					} else {	
144						if(runif(1)<0.2 | sum(cur.delta.theta)==0) {			# adjust theta
145							new.alpha=cur.alpha 
146							new.theta=adjust.value(cur.theta) ##
147							new.delta.theta=cur.delta.theta
148							new.root.theta<-nrt<-assign.root.theta(new.theta, ape.tre) ## 
149							new.regimes=as.factor(c(nrt,new.theta)) ##
150							new.rates=cur.rates
151							new.delta.rates=cur.delta.rates
152							nThetaProp=nThetaProp+1	
153						} else {					# swap theta
154							nr=switcheroo(cur.delta.theta, cur.theta, ape.tre, node.des, theta=TRUE)
155							new.alpha=cur.alpha 
156							new.theta=nr$new.values ##
157							new.delta.theta=nr$new.delta ##
158							new.root.theta<-nrt<-assign.root.theta(new.theta, ape.tre) ## 
159							new.regimes=as.factor(c(nrt,new.theta)) ##
160							new.rates=cur.rates
161							new.delta.rates=cur.delta.rates
162							nThetaSwapProp=nThetaSwapProp+1	
163						}
164					}
165				}
166			}
167			
168			new.tre=apply.BMM.to.apetree(ape.tre, new.rates)
169			mod.new=new.ou.lik.fn(new.rates, new.alpha, new.regimes, new.theta, new.tre, ouch.dat)
170		} else if(cur.model=="BM"){
171## BM IMPLEMENTATION ##
172			if (runif(1) < (2 * probMergeSplit)) {							# adjust rate categories
173				nr=split.or.merge(cur.delta.rates, cur.rates, ape.tre, node.des, theta=FALSE)
174				new.rates=nr$new.values
175				new.delta.rates=nr$new.delta
176				new.root=cur.root
177				nRcatProp=nRcatProp+1
178				lnHastingsRatio=nr$lnHastingsRatio
179				lnPriorRatio=nr$lnPriorRatio
180			} else { 
181				if(runif(1)<probRoot | sum(cur.delta.rates)==0) {		# adjust root
182					new.root=adjust.value(cur.root)
183					new.rates=cur.rates
184					new.delta.rates=cur.delta.rates
185					nRootProp=nRootProp+1
186				} else {													# adjust rates
187					if(runif(1)<0.05) {
188						new.rates=adjust.rate(cur.rates)
189						new.delta.rates=cur.delta.rates
190						new.root=cur.root
191						nRateProp=nRateProp+1
192					} else {												# swap rate classes
193						nr=switcheroo(cur.delta.rates, cur.rates, ape.tre, node.des, theta=FALSE)
194						new.rates=nr$new.values ##
195						new.delta.rates=nr$new.delta ##
196						new.root=cur.root
197						nRateSwapProp=nRateSwapProp+1
198					}
199				}
200			}
201			
202			new.tre=apply.BMM.to.apetree(ape.tre, new.rates)
203			mod.new=try(new.bm.lik.fn(new.rates, new.root, new.tre, orig.dat),silent==TRUE)
204		}
205		if(!inherits(mod.new, "try-error")) {
206			lnLikelihoodRatio = mod.new$lnL - mod.cur$lnL
207		} else {
208			lnLikelihoodRatio = -Inf
209			mod.new$lnL=NaN
210			failures="./failed.model.parameters/"
211			failed.trees=paste(failures, "trees", sep="/")
212			lapply(list(failures, failed.trees), function(x) if(!file.exists(x)) dir.create(x))
213				
214			pdf(file=paste(failed.trees, paste(i,"failed.model.pdf",sep="."),sep="/"))
215				plot(new.tre)
216			dev.off()
217			write.table(i,file=paste(failures, "failed.rates.txt",sep="/"),append=TRUE,col=FALSE,row=FALSE,quote=FALSE)
218			write.table(matrix(new.rates,nrow=1),file=paste(failures,"failed.rates.txt",sep="/"),append=TRUE,col=FALSE,row=FALSE,quote=FALSE)
219			write.table(matrix(c(i, new.root),nrow=1),file=paste(failures,"failed.root.txt",sep="/"),append=TRUE,col=FALSE,row=FALSE,quote=FALSE)
220		}
221		
222	# compare likelihoods
223		endparms = c(nRateProp, nRateSwapProp, nRcatProp, nRootProp, nAlphaProp, nThetaProp, nThetaSwapProp, nTcatProp)
224		
225		r=assess.lnR((heat * lnLikelihoodRatio + heat * lnPriorRatio + lnHastingsRatio)->lnR)
226
227		if (runif(1) <= r$r) {			# adopt proposal
228			if(cur.model=="OU") {
229				decision="adopt"
230				cur.alpha <- new.alpha
231				cur.root.theta <- new.root.theta
232				cur.regimes <- new.regimes
233				cur.theta <- new.theta
234				cur.delta.theta <- new.delta.theta
235			} else {
236				decision="adopt"
237				cur.root <- new.root
238			}
239			cur.rates <- new.rates
240			cur.delta.rates <- new.delta.rates
241			curr.lnL <- mod.new$lnL
242			cur.tre <- new.tre
243			mod.cur <- mod.new
244			cOK <- determine.accepted.proposal(startparms, endparms, cOK)
245		} else {						# deny proposal	
246			decision="reject"
247			curr.lnL <- mod.cur$lnL 
248		}
249		
250		if(is.infinite(curr.lnL)) stop("starting point has exceptionally poor likelihood")
251
252		if(r$error) generate.error.message(i, mod.cur, mod.new, lnR, errorLog)
253		
254		if(i%%tickerFreq==0) {
255			if(i==tickerFreq) cat("|",rep(" ",9),toupper("generations complete"),rep(" ",9),"|","\n")
256			cat(". ")
257		}
258		
259		if(i%%sampleFreq==0) {
260			bundled.parms=list(gen=i, mrate=exp(mean(log(cur.rates))), cats=sum(cur.delta.rates)+1, root=cur.root, alpha=cur.alpha, reg=(sum(cur.delta.theta)+1), theta=mean(cur.theta), lnL=curr.lnL)
261			generate.log(bundled.parms, cur.model, file=runLog)
262			ll$lnL=rbind(ll$lnL, curr.lnL)
263			ll$rates=rbind(ll$rates, cur.rates) 
264			if(cur.model=="OU") {
265				ll$theta=rbind(ll$theta, cur.theta)
266				ll$alpha=rbind(ll$alpha, cur.alpha)
267				ll$root=rbind(ll$root, cur.root.theta)
268			}
269			if(cur.model=="BM") {
270				ll$root=rbind(ll$root, cur.root)
271			}
272			if(i==ngen) {
273				if(cur.model=="OU") dimnames(ll$theta)[[2]]=ape.tre$edge[,2]
274				dimnames(ll$rates)[[2]]=ape.tre$edge[,2]
275			}
276		}
277	}
278	
279# End rjMCMC
280	close(errorLog)
281	close(runLog)
282	summarize.run(cOK, endparms, cur.model)	
283	return(list(results=ll,data=orig.dat,tree=ape.tre))
284}
285
286
287## AUXILIARY FUNCTIONS
288
289
290new.ou.lik.fn <- function(rates, alpha, regimes, theta, ape.tre, ouch.dat) { # from OUCH 2.7-1:::hansen()
291# mod 09.13.2010 JM Eastman
292	ouch.tre=ape2ouch(ape.tre, scale=FALSE)
293	ouch.dat$regimes=as.factor(regimes)
294	theta.tmp=as.numeric(levels(as.factor(theta)))
295	theta=theta.tmp[order(match(theta.tmp,theta))]
296	alpha.ouch=as.matrix(sqrt(alpha))
297	beta=ouch:::regime.spec(ouch.tre, ouch.dat['regimes'])
298	dat=do.call(c,lapply(ouch.dat['trait'],function(y)y[ouch.tre@term]))
299	names(dat)=ouch.dat$labels[Ntip(ape.tre):nrow(ouch.dat)]
300	n <- length(dat)
301	otree.df=as(ouch.tre,"data.frame")
302	ordering=otree.df[n:nrow(otree.df),'labels']
303	ev <- eigen(alpha.ouch,symmetric=TRUE)
304	w <- .Call(ouch:::ouch_weights,object=ouch.tre,lambda=ev$values,S=ev$vectors,beta=beta)
305	v <- geiger:::ouMatrix(vcv.phylo(ape.tre), alpha)
306	v.match <- match(row.names(v),ordering)
307	v <- v[v.match,v.match]
308		## GEIGER call does not assume ultrametricity (and thus allows different 'rates' under BM)
309		# -old-	v <- .Call(ouch:::ouch_covar,object=ouch.tre,lambda=ev$values,S=ev$vectors,sigma.sq=1)
310	
311	y <- matrix(theta,ncol=1)
312	e <- w%*%y-dat
313	dim(y) <- dim(y)[1]
314	dim(e) <- n
315	q <- e%*%solve(v,e)
316	det.v <- determinant(v,logarithm=TRUE)
317	if (det.v$sign!=1)stop("ou.lik.fn error: non-positive determinant",call.=FALSE)
318	dev <- n*log(2*pi)+as.numeric(det.v$modulus)+q[1,1]
319	list(
320		 theta=theta,
321		 weight=w,
322		 vcov=v,
323		 resids=e,
324		 lnL=-0.5*dev
325		 )
326}
327
328new.bm.lik.fn <- function(rates, root, ape.tre, orig.dat) { # from OUCH 2.7-1:::brown()	
329# mod 09.16.2010 JM Eastman
330
331	tree=ape.tre
332	data=orig.dat
333
334	nterm=Ntip(tree)
335	nchar=length(data)
336	dat=data[match(tree$tip.label,names(data))]
337
338	nchar <- 1
339	w <- matrix(data=1,nrow=nterm,ncol=1)
340	n <- length(dat)
341
342	b <- vcv.phylo(ape.tre)
343
344	y <- matrix(root)
345	e <- w%*%y-dat
346	q <- t(e)%*%solve(b,e)
347	v <- q/nterm
348	dev <- nchar*nterm*(1+log(2*pi))+nchar*log(det(b))+nterm*log(det(v))
349	list(
350		 root=root,
351		 lnL = -0.5*dev
352		 )	
353}
354
355split.or.merge <- function(cur.delta, cur.values, phy, node.des, theta=FALSE) { 
356# mod 09.17.2010 JM Eastman
357	bb=cur.delta
358	vv=cur.values
359	names(vv)<-names(bb)<-phy$edge[,2]
360	new.bb=choose.one(bb)
361	new.vv=vv
362	
363	s=names(new.bb[bb!=new.bb])
364	all.shifts=as.numeric(names(bb[bb>0]))
365	all.D=node.des[[which(names(node.des)==s)]]
366	if(length(all.D)!=0) {
367		untouchable=unlist(lapply(all.shifts[all.shifts>s], function(x)node.des[[which(names(node.des)==x)]]))
368		remain.unchanged=union(all.shifts, untouchable)
369	} else {
370		untouchable=NULL
371		remain.unchanged=list()
372	}
373	
374	marker=match(s, names(new.vv))
375	nn=length(vv)
376	K=sum(bb)+1
377	N=Ntip(phy)
378	
379	if(sum(new.bb)>sum(bb)) {			# add transition: SPLIT
380#		print("s")
381		if(theta) {							# assign new theta to current node
382			new.vv[marker] <- nr <- adjust.value(new.vv[marker])
383		} else {							# assign new rate to current node
384			new.vv[marker] <- nr <- adjust.rate(new.vv[marker])
385		}
386		if(length(remain.unchanged)==0) {	# assign new value to all descendants
387			new.vv[match(all.D,names(new.vv))] = nr
388		} else {							# assign new value to all 'open' descendants 
389			new.vv[match(all.D[!(all.D%in%remain.unchanged)],names(new.vv))] = nr
390		}
391		lnHastingsRatio = log((K+1)/(2*N-2-K)) ### from Drummond and Suchard 2010: where N is tips, K is number of local parms in tree
392		lnPriorRatio = log(ppois(K+1,log(2),lower.tail=FALSE)/ppois(K,log(2),lower.tail=FALSE))
393	} else {							# drop transition: MERGE
394#		print("m")
395		anc = get.ancestor.of.node(s, phy)
396		if(!is.root(anc, phy)) {			# substitute ancestral value for current
397			new.vv[match(s, names(new.vv))] <- nr <- new.vv[match(anc, names(new.vv))] 
398		} else {
399			sister.tmp=get.desc.of.node(anc,phy)
400			sister=sister.tmp[sister.tmp!=s]
401			if(bb[match(sister, names(bb))]==0) {
402				nr=as.numeric(vv[match(sister, names(vv))])
403			} else {
404				if(theta) {						# assign new theta to first descendant of root
405					new.vv[marker] <- nr <- adjust.value(new.vv[marker])
406				} else {						# assign new rate to first descendant of root
407					new.vv[marker] <- nr <- adjust.rate(new.vv[marker])
408				}
409			}
410		}
411		if(length(remain.unchanged)==0) {	# assign new value to all descendants
412			new.vv[match(all.D, names(new.vv))] = nr
413		} else {							# assign new value to all 'open' descendants
414			new.vv[match(all.D[!(all.D%in%remain.unchanged)],names(new.vv))] = nr
415		}
416		lnHastingsRatio = log((2*N-2-K+1)/K) ### from Drummond and Suchard 2010: where N is tips, K is number of local parms in tree
417		lnPriorRatio = log(ppois(K-1,log(2),lower.tail=FALSE)/ppois(K,log(2),lower.tail=FALSE))
418		
419	}
420	if(((length(table(new.vv))-1)-sum(new.bb))!=0) fail=TRUE else fail=FALSE
421	return(list(new.delta=new.bb, new.values=new.vv, lnHastingsRatio=lnHastingsRatio, lnPriorRatio=lnPriorRatio, fail=fail))
422}
423
424switcheroo <- function(cur.delta, cur.values, phy, node.des, theta=FALSE) { 
425	bb=cur.delta
426	vv=cur.values
427	new.bb=sample(bb)
428	names(vv)<-names(bb)<-names(new.bb)<-phy$edge[,2]
429	new.vv=vv
430	
431	old.shifts=as.numeric(names(bb[bb!=0]))
432	old.values=vv[match(old.shifts, names(vv))]
433	
434	anc=lapply(old.shifts, function(x)get.ancestor.of.node(x,phy))
435	if(!is.root(min(unlist(anc)),phy)) {
436	# find most ancestral value in tree
437		anc.value=unname(vv[match(min(unlist(anc)),names(vv))]) 
438	} else {
439		if(theta) anc.value=mean(sapply(old.values, function(x) adjust.value(x))) else anc.value=mean(sapply(old.values, function(x) adjust.rate(x)))
440	}
441	all.desc.anc=node.des[[which(names(node.des)==min(unlist(anc)))]]
442	new.vv[match(all.desc.anc,names(new.vv))] = anc.value
443
444	new.shifts=sort(as.numeric(names(new.bb[new.bb!=0])))
445	new.tips=new.shifts[new.shifts<=Ntip(phy)]
446	new.ints=new.shifts[new.shifts>Ntip(phy)]
447	new.shifts=c(new.ints, new.tips)
448	
449	new.values=as.numeric(old.values)[sample(1:length(old.values))]
450	for(z in 1:length(new.shifts)) {
451		node=new.shifts[z]
452		new.vv[match(node, names(new.vv))]=new.values[z]
453		new.vv[match(node.des[[which(names(node.des)==node)]],names(new.vv))]=new.values[z]
454	} 
455	return(list(new.delta=new.bb, new.values=new.vv))
456}
457
458generate.starting.point <- function(data, phy, node.des, theta=FALSE) { 
459	nn=length(phy$edge.length)
460	ntip=Ntip(phy)
461	nshifts=rtpois(1,log(ntip),nn)
462	if(nshifts!=0) bb=sample(c(rep(1,nshifts),rep(0,nn-nshifts)),replace=FALSE) else bb=rep(0,nn)
463	names(bb)=phy$edge[,2]
464	shifts=as.numeric(names(bb)[bb==1])
465	if(theta) {
466		min.max=c(adjust.value(min(data)), adjust.value(max(data)))
467		values=runif(sum(bb)+1, min=min.max[min.max==min(min.max)], max=min.max[min.max==max(min.max)])
468	} else {
469		init.rate=fitContinuous(phy,data)[[1]]$beta
470		min.max=c(adjust.rate(init.rate), adjust.rate(init.rate))
471		values=runif(sum(bb)+1, min=min.max[min.max==min(min.max)], max=min.max[min.max==max(min.max)])
472	}
473	internal.shifts<-tip.shifts<-numeric(0)
474	internal.shifts=sort(shifts[shifts>ntip])
475	tip.shifts=shifts[shifts<=ntip]
476	
477	if(length(internal.shifts)==0 & length(tip.shifts)==0) {
478		vv=rep(values, nn)
479	} else {
480		vv=bb
481		vv[]=values[length(values)]
482		i=0
483		if(length(internal.shifts)!=0) {
484			for(i in 1:length(internal.shifts)) {
485				d=node.des[which(names(node.des)==internal.shifts[i])]
486#		d=get.descendants.of.node(internal.shifts[i],phy)
487				vv[match(c(internal.shifts[i], d), names(vv))]=values[i]
488			}
489		}
490		if(length(tip.shifts)!=0) {
491			for(j in 1:length(tip.shifts)) {
492				vv[match(tip.shifts[j], names(vv))]=values[j+i]
493			}
494		}
495	}
496	return(list(delta=unname(bb), values=unname(vv)))
497}
498
499assign.root.theta <- function(cur.theta, phy) {
500	root=phy$edge[1,1]
501	desc=phy$edge[which(phy$edge[,1]==root),2]
502	root.theta=cur.theta[sample(which(phy$edge[,2]==desc),1)]
503	root.theta
504}
505
506get.descendants.of.node <- function(node, phy) {
507	storage=list()
508	if(node%in%phy$edge[,1]) {
509		desc=c(sapply(node, function(x) {return(phy$edge[which(phy$edge[,1]==x),2])}))
510		while(any(desc%in%phy$edge[,1])) {
511			desc.true=desc[which(desc%in%phy$edge[,1])]
512			desc.desc=lapply(desc.true, function(x) return(get.desc.of.node(x, phy)))
513			
514			storage=list(storage,union(desc,desc.desc))
515			
516			desc=unlist(desc.desc)
517		}
518		storage=sort(unique(unlist(union(desc,storage))))
519	}
520	return(storage)
521}
522
523get.desc.of.node <- function(node, phy) {
524	return(phy$edge[which(phy$edge[,1]==node),2])
525}
526
527get.ancestor.of.node<-function(node, phy) {
528	anc=phy$edge[which(phy$edge[,2]==node),1]
529	anc
530}
531
532choose.one <- function(cur.delta)
533{
534	bb=cur.delta
535	s=sample(1:length(bb),1)
536	bb[s]=1-bb[s]
537	bb
538}
539
540is.root<-function(node,phy) {
541	if(node==phy$edge[1,1]) return(TRUE) else return(FALSE)
542}
543
544assess.lnR <- function(lnR) {
545	if(is.na(lnR) || abs(lnR)==Inf) {
546		error=TRUE
547		r=0
548	} else {
549		if(lnR < -20) {
550			r=0 
551		} else if(lnR > 0) {
552			r=1 
553		} else {
554			r=exp(lnR)
555		}
556		error=FALSE
557	}
558	return(list(r=r, error=error))
559}
560
561apply.BMM.to.apetree<-function(apetree, rates) {
562	apetree$edge.length=apetree$edge.length*rates
563	return(apetree)
564}
565
566adjust.value <- function(value) {
567	vv=value
568	vv.levels=levels(factor(vv))
569	rch <- runif(length(vv.levels), min = -1, max = 1)
570	for(i in 1:length(vv.levels)){
571		vv[which(vv==vv.levels[i])]=vv[which(vv==vv.levels[i])]+rch[i]
572	}
573	vv
574}
575
576adjust.rate <- function(rate) {
577	vv=rate
578	v=log(vv)
579	vv.levels=levels(factor(vv))
580	rch <- runif(length(vv.levels), min = -1, max = 1)
581	for(i in 1:length(vv.levels)){
582		v[which(vv==vv.levels[i])]=v[which(vv==vv.levels[i])]+rch[i]
583	}
584	v=exp(v)
585	v
586}
587	
588prepare.ouchdata <- function(phy, data) {
589	td <- treedata(phy, data, sort = T)	
590	ouch.tre <- ape2ouch(td$phy->ape.tre, scale=FALSE)		
591	ouch.dat=data.frame(as.numeric(ouch.tre@nodes), as.character(ouch.tre@nodelabels))
592	names(ouch.dat)=c("nodes","labels")
593	ouch.dat[ouch.dat[,2]=="",2]=NA
594	ouch.dat$trait=NA
595	mM=match(ouch.dat$labels,names(data))
596	for(i in 1:nrow(ouch.dat)){
597		if(!is.na(mM[i])){
598			ouch.dat$trait[i]=data[mM[i]]	
599		}
600	}
601	return(list(ouch.tre=ouch.tre, ouch.dat=ouch.dat, ape.tre=td$phy, orig.dat=td$data[,1]))
602}
603
604summarize.run<-function(cOK, endparms, cur.model) {
605#	end = c(nRateProp, nRateSwapProp, nRcatProp, nRootProp, nAlphaProp, nThetaProp, nThetaSwapProp, nTcatProp)
606
607	names(endparms)<-names(cOK)<-c("rates", "rate.swap", "rate.catg", "root", "alpha", "theta", "theta.swap", "theta.catg")
608	names(endparms)=paste("pr.",names(endparms),sep="")
609	names(cOK)=paste("ok.",names(cOK),sep="")
610	cat("\n")
611	if(cur.model=="BM") {
612		print(endparms[paste("pr.", c("rates", "rate.swap", "rate.catg", "root"),sep="")])
613		print(cOK[paste("ok.", c("rates", "rate.swap", "rate.catg", "root"),sep="")])
614	} else {
615		print(endparms[paste("pr.", c("rates", "rate.swap", "rate.catg", "alpha", "theta", "theta.swap", "theta.catg"),sep="")])
616		print(cOK[paste("ok.", c("rates", "rate.swap", "rate.catg", "alpha", "theta", "theta.swap", "theta.catg"),sep="")])
617	}
618}
619
620generate.log<-function(bundled.parms, cur.model, file, init=FALSE) {
621#	bundled.parms=list(gen=i, mrate=mean(cur.rates), cats=sum(cur.delta.rates)+1, root=cur.root, alpha=cur.alpha, reg=sum(cur.delta.theta)+1, theta=mean(cur.theta), lnL=curr.lnL)
622	res=bundled.parms
623
624	if(init) {
625		if(cur.model=="BM") msg=paste("gen", "model", "mean.rate", "rates", "root", "lnL", sep="\t") else msg=paste("gen", "model", "mean.rate", "rates", "alpha", "regimes", "mean.theta", "lnL", sep="\t")
626	} else {
627		if(cur.model=="BM") {
628			msg<-paste(res$gen, sQuote(cur.model), sprintf("%.3f", res$mrate), res$cats, sprintf("%.3f", res$root), sprintf("%.3f", res$lnL),sep="\t")
629		} else {
630			msg<-paste(res$gen, sQuote(cur.model), sprintf("%.3f", res$mrate), res$cats, sprintf("%.4f", res$alpha), res$reg, sprintf("%.3f", res$theta), sprintf("%.3f", res$lnL),sep="\t")
631		}
632	}
633	write(msg, file=file, append=TRUE)
634}
635
636generate.error.message<-function(i, mod.cur, mod.new, lnR, errorLog, init=FALSE) {
637	if(init) {
638		initmsg = write(paste("gen", "curr.lnL", "new.lnL", "lnR", sep="\t"), file=errorLog)
639	} else {
640		write(paste(i, sprintf("%.3f", mod.cur$lnL), sprintf("%.3f", mod.new$lnL), sprintf("%.3f", lnR), sep="\t"), file=errorLog)
641	}
642}
643
644determine.accepted.proposal<-function(startparms, endparms, cOK) {
645	which(startparms!=endparms)->marker
646	cOK[marker]=cOK[marker]+1
647	cOK
648}
649
650rtpois<-function(N, lambda, k) {
651	p=ppois(k, lambda, lower.tail=FALSE)
652	out=qpois(runif(N, min=0, max=1-p), lambda)
653	out
654}
655				 
656
657				 
658