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