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