/archive/source.archive/bm.and.ou/trait.rjmcmc.09132010a.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, probRegimes=0.10, probMergeSplit=0.05, probRoot=0.1, heat=1, fileBase="result") 
  5{
  6	if(length(model)>1)stop(paste("Please specify either ", sQuote("OU"), " or ", sQuote("BM")))
  7	
  8	require(ouch)
  9	require(geiger)
 10	
 11### prepare objects for rjMCMC
 12	currModel		<- model								# initialize family of models explored
 13
 14	dataList		<- prepare.ouchdata(phy, data)			# check tree and data; convert into S4
 15
 16	ape.tre			<- currTree <- dataList$ape.tre			# S3 phylogeny
 17	orig.dat		<- dataList$orig.dat					# S3 data	
 18	ouch.tre		<- dataList$ouch.tre					# S4 OUCH phylogeny
 19	ouch.dat		<- dataList$ouch.dat					# S4 OUCH data
 20	
 21	nn				<- Ntip(ape.tre)						# max number of rates (one per branch)
 22	
 23	currRates.c		<- rep(1,nn)							# initialize array for recording which branch belongs to which rate category
 24    currRates		<- rep(fitContinuous(ape.tre,orig.dat)[[1]]$beta,nn)							# initialize array for rate estimates branchwise
 25	currRoot		<- mean(orig.dat)
 26    delta.rates		<- rep(0,nn)
 27	currAlpha		<- 0.001								# initialize alpha
 28	currTheta		<- 0									# initialize trait optima under OU
 29	delta.theta		<- rep(0,nn)							# initialize local 'rates'
 30	currRegimes		<- rep(1,nrow(ouch.dat))				# initial set of selective regimes
 31	
 32	nRateProp = 0											# proposal count: update rate parameter(s)
 33	nRcatProp = 0											# proposal count: update rate categories
 34	nRootProp = 0
 35	nAlphaProp = 0											# proposal count: update alpha parameter
 36	nRegimesProp = 0										# proposal count: update selective regimes
 37	nThetaProp = 0
 38	nRateOK = 0												# proposal count: successful rate changes
 39	nRcatOK = 0												# proposal count: successful rate categories updates
 40	nRootOK = 0
 41	nAlphaOK = 0											# proposal count: successful alpha updates
 42	nRegimesOK = 0											# proposal count: successful regime updates
 43	nThetaOK = 0
 44	l <- array(dim=c(ngen,4))								# array for storing output: size is number of MCMC generations
 45    l[] <- NA												# initialize array for storing lnLs
 46	
 47	cOK=c(													# list of proposal tallies
 48		  nRateOK, 
 49		  nRcatOK, 
 50		  nRootOK, 
 51		  nAlphaOK, 
 52		  nRegimesOK, 
 53		  nThetaOK
 54		  )		
 55	
 56	errorLog=file(paste(getwd(),paste(currModel, "rjTraitErrors.log",sep="."),sep="/"),open='w+')
 57	generate.error.message(i=NULL, modCurr=NULL, modNew=NULL, lnR=NULL, errorLog=errorLog, init=TRUE)
 58	runLog=file(paste(getwd(),paste(currModel, "rjTrait.log",sep="."),sep="/"),open='w+')
 59	generate.log(bundled.parms=NULL, currModel, file=runLog, init=TRUE)
 60	
 61### Begin rjMCMC
 62    for (i in 1:ngen) {
 63		error.flag=FALSE
 64        lnLikelihoodRatio <- lnProposalRatio <- lnPriorRatio <- 0
 65		startparms=list(rates=currRates, cats=currRates.c, root=currRoot, alph=currAlpha, reg=currRegimes, theta=currTheta)
 66		
 67## OU IMPLEMENTATION ## 
 68		if(currModel=="OU") {
 69			if (runif(1) < (2 * probMergeSplit)) {
 70				if(runif(1)<0.5) { # adjust regimes
 71					nr=split.or.merge(currTheta, currRegimes, decide.s.or.m(currRegimes, regimes=TRUE)->s.m, regimes=TRUE)
 72					newRegimes=nr$new.categories
 73					newTheta=nr$new.values
 74					newRates=currRates
 75					newRates.c=currRates.c
 76					nRegimesProp=nRegimesProp+1
 77				} else { # adjust rate categories
 78					nr=split.or.merge(currRates, currRates.c, decide.s.or.m(currRates.c)->s.m)
 79					newRates=nr$new.values
 80					newRates.c=nr$new.categories
 81					newRegimes=currRegimes
 82					newTheta=currTheta
 83					nRcatProp=nRcatProp+1
 84				}
 85				newAlpha=currAlpha
 86				lnProposalRatio=nr$lnl.prop
 87				lnPriorRatio=nr$lnl.prior
 88			} else { # neither split nor merge
 89				if(runif(1)<probAlpha) { # adjust alpha
 90					newAlpha=adjust.rates(currAlpha)
 91					newRates=currRates
 92					newRates.c=currRates.c
 93					newRegimes=currRegimes
 94					newTheta=currTheta
 95					nAlphaProp=nAlphaProp+1
 96				} else { 
 97					if(runif(1)<0.5) {# adjust rates
 98						newRates=adjust.rates(currRates)
 99						newRates.c=fix.categories(newRates)
100						newAlpha=currAlpha
101						newRegimes=currRegimes
102						newTheta=currTheta
103						nRateProp=nRateProp+1
104					} else { # adjust theta
105						newTheta=sapply(currTheta, adjust.value)
106						newAlpha=currAlpha
107						newRates=currRates
108						newRates.c=currRates.c
109						newRegimes=currRegimes
110						nThetaProp=nThetaProp+1
111					}
112				}
113			}
114			
115			currTree=apply.BMM.to.apetree(ape.tre, currRates)
116			newTree=apply.BMM.to.apetree(ape.tre, newRates)
117			modCurr=new.ou.lik.fn(currRates, currAlpha, currRegimes, currTheta, currTree, ouch.dat)
118			modNew=try(new.ou.lik.fn(newRates, newAlpha, newRegimes, newTheta, newTree, ouch.dat),silent=TRUE)
119		} else if(currModel=="BM"){
120## BM IMPLEMENTATION ##
121			if (runif(1) < (2 * probMergeSplit)) {
122				nr=split.or.merge(currRates, currRates.c, decide.s.or.m(currRates.c)->s.m)
123				newRates=nr$new.values
124				newRates.c=nr$new.categories
125				newRoot=currRoot
126				nRcatProp=nRcatProp+1
127				lnProposalRatio=nr$lnl.prop
128				lnPriorRatio=nr$lnl.prior
129			} else { 
130				if(runif(1)<probRoot) { # adjust root
131					newRoot=adjust.value(currRoot)
132					newRates=currRates
133					newRates.c=currRates.c
134					nRootProp=nRootProp+1
135				} else { # adjust rates
136					newRates=adjust.rates(currRates)
137					newRates.c=fix.categories(newRates)
138					newRoot=currRoot
139					nRateProp=nRateProp+1
140				}
141			}
142			
143			currTree=apply.BMM.to.apetree(ape.tre, currRates)
144			newTree=apply.BMM.to.apetree(ape.tre, newRates)
145			modCurr=new.bm.lik.fn(currRates, currRoot, currTree, ouch.dat)
146			modNew=try(new.bm.lik.fn(newRates, newRoot, newTree, ouch.dat),silent=TRUE)
147		}
148	
149		if(!inherits(modNew, "try-error")) {
150			lnLikelihoodRatio = modNew$lnL - modCurr$lnL
151		} else {
152			lnLikelihoodRatio = -Inf
153			modNew$lnL=NaN
154		}
155		r=assess.lnR((heat * lnLikelihoodRatio + heat * lnPriorRatio + lnProposalRatio)->lnR)
156		
157		if (runif(1) <= r$r) {			# adopt proposal
158			if(currModel=="OU") {
159				currAlpha <- newAlpha
160				currRegimes <- newRegimes
161				currTheta <- newTheta
162			} else {
163				currRoot <- newRoot
164			}
165			currRates <- newRates
166			currRates.c <- newRates.c
167			curr.lnL <- modNew$lnL
168			currTree <- newTree
169		} else {						# deny proposal			
170			curr.lnL <- modCurr$lnL 
171		}
172		if(r$error) generate.error.message(i, modCurr, modNew, lnR, errorLog)
173		
174		l[i,]<-c(curr.lnL, max(currRates.c), currAlpha, max(currRegimes))
175		endparms=list(rates=currRates, cats=currRates.c, root=currRoot, alph=currAlpha, reg=currRegimes, theta=currTheta)
176		cOK=tally.mcmc.parms(startparms, endparms, cOK)
177		
178		bundled.parms=list(gen=i, mrate=exp(mean(log(currRates))), cats=max(currRates.c), root=currRoot, alpha=currAlpha, reg=max(currRegimes), theta=mean(currTheta), lnL=curr.lnL)
179		if(i%%sampleFreq==0) generate.log(bundled.parms, currModel, file=runLog)
180	}
181	
182# End rjMCMC
183	close(errorLog)
184	close(runLog)
185	summarize.run(cOK, (cProp<-c(nRateProp, nRcatProp, nRootProp, nAlphaProp, nRegimesProp, nThetaProp)), currModel)	
186	if(currModel=="BM") {
187		return(list(r.cat=l[,2], lnL=l[,1]))
188	} else {
189		return(list(r.cat=l[,2], r.alpha=l[,3], r.regimes=l[,4], lnL=l[,1]))
190	}
191}
192
193
194## AUXILIARY FUNCTIONS
195
196
197new.ou.lik.fn <- function(rates, alpha, regimes, theta, ape.tre, ouch.dat) { # from OUCH 2.7-1:::hansen()
198# mod 09.13.2010
199	ouch.tre=ape2ouch(ape.tre, scale=FALSE)
200	ouch.dat$regimes=as.factor(regimes)
201	alpha.ouch=as.matrix(sqrt(alpha))
202	beta=ouch:::regime.spec(ouch.tre, ouch.dat['regimes'])
203	dat=do.call(c,lapply(ouch.dat['trait'],function(y)y[ouch.tre@term]))
204	n <- length(dat)
205	otree.df=as(ouch.tre,"data.frame")
206	ordering=otree.df[n:nrow(otree.df),'labels']
207	ev <- eigen(alpha.ouch,symmetric=TRUE)
208	w <- .Call(ouch:::ouch_weights,object=ouch.tre,lambda=ev$values,S=ev$vectors,beta=beta)
209
210	v <- geiger:::ouMatrix(vcv.phylo(ape.tre), alpha)
211	v.match <- match(row.names(v),ordering)
212	v <- v[v.match,v.match]
213		## GEIGER call does not assume ultrametricity (and thus allows different 'rates' under BM)
214		# -old-	v <- .Call(ouch:::ouch_covar,object=ouch.tre,lambda=ev$values,S=ev$vectors,sigma.sq=1)
215	
216	y <- matrix(theta,ncol=1)
217	e <- w%*%y-dat
218	dim(y) <- dim(y)[1]
219	dim(e) <- n
220	q <- e%*%solve(v,e)
221	det.v <- determinant(v,logarithm=TRUE)
222	if (det.v$sign!=1)stop("ou.lik.fn error: non-positive determinant",call.=FALSE)
223	dev <- n*log(2*pi)+as.numeric(det.v$modulus)+q[1,1]
224	list(
225		 theta=theta,
226		 weight=w,
227		 vcov=v,
228		 resids=e,
229		 lnL=-0.5*dev
230		 )
231}
232
233new.bm.lik.fn <- function(rates, root, ape.tre, ouch.dat) { # from OUCH 2.7-1:::brown()	
234# mod 09.13.2010	
235	tree=ape2ouch(ape.tre, scale=FALSE)
236	data=ouch.dat['trait']
237	
238	if (!is(tree, "ouchtree")) {stop(sQuote("tree"), " must be an object of class ", sQuote("ouchtree"))}
239    if (is.data.frame(data)) {
240        nm <- rownames(data)
241        data <- lapply(as.list(data), function(x) {names(x) <- nm; x})
242    }
243    if (is.numeric(data)) {
244        nm <- deparse(substitute(data))[1]
245        data <- list(data)
246        names(data) <- nm
247    }
248    if (is.list(data)) {
249        if (any(sapply(data, class) != "numeric") || any(sapply(data, length) != tree@nnodes)) {stop(sQuote("data"), " vector(s) must be numeric, with one entry per node of the tree")}
250        if (any(sapply(data, function(x) (is.null(names(x))) || (!setequal(names(x), tree@nodes))))) {stop(sQuote("data"), " vector names (or data-frame row names) must match node names of ", sQuote("tree"))}
251        for (xx in data) {
252            no.dats <- which(is.na(xx[tree@nodes[tree@term]]))
253            if (length(no.dats) > 0) {stop("missing data on terminal node(s): ", paste(tree@nodes[tree@term[no.dats]], collapse = ","))}
254        }
255    }
256    nterm <- tree@nterm
257    nchar <- length(data)
258    if (is.null(names(data))) names(data) <- paste("char", 1:nchar, sep = "")
259    dat.tmp <- lapply(data, function(x) x[tree@nodes])
260    dat <- lapply(dat.tmp, function(y) y[tree@term])
261	dat <- dat[[1]]
262	nterm <- Ntip(ape.tre)
263	nchar <- 1
264	w <- matrix(data=1,nrow=nterm,ncol=1)
265	b <- tree@branch.times
266	y <- matrix(root)
267	e <- w%*%y-dat
268	q <- t(e)%*%solve(b,e)
269	v <- q/nterm
270	dev <- nchar*nterm*(1+log(2*pi))+nchar*log(det(b))+nterm*log(det(v))
271	list(
272		 root=root,
273		 lnL = -0.5*dev
274		 )	
275}
276
277rtpois<-function(N, lambda, k){
278	p=ppois(k, lambda, lower.tail=FALSE)
279	out=qpois(runif(N, min=0, max=1-p), lambda)
280	out
281}
282
283get.descendants.of.node <- function(node, phy) {
284	descendants.of.node <- function(node, phy) {
285		current.node=node[node%in%phy$edge[,1]]
286		if (length(current.node)>0) {
287			desc=c(sapply(current.node, function(x) {return(phy$edge[which(phy$edge[,1]==x),2])}))
288			write(desc,file="TEMP",append=TRUE)
289			descendants.of.node(desc,phy)
290		}	
291	}
292	descendants.of.node(node, phy)
293	if(file.exists("TEMP")) {
294		out=unlist(strsplit(readLines("TEMP")," "))
295		unlink("TEMP")
296		return(sort(unique(as.numeric(out))))
297	}
298}	
299
300get.ancestor.of.node<-function(node, phy) {
301	anc=phy$edge[which(phy$edge[,2]==node),1]
302	anc
303}
304
305assess.lnR <- function(lnR) {
306	if(is.na(lnR) || abs(lnR)==Inf) {
307		error=TRUE
308		r=0
309	} else {
310		if(lnR < -100) {
311			r=0 
312		} else if(lnR > 0) {
313			r=1 
314		} else {
315			r=exp(lnR)
316		}
317		error=FALSE
318	}
319	return(list(r=r, error=error))
320}
321
322decide.s.or.m <- function(categories, regimes=FALSE) {
323	if (max(categories) == 1) {
324		return("split")
325	}
326	else if (max(categories) == length(categories)) {
327		return("merge")
328	}
329#	else if (regimes && (((length(categories)+1)/2)==max(categories))) {
330# limits number of selective regimes to number of tips
331#		return("merge")
332#	}
333	else if (runif(1) < 0.5) {
334		return("split")
335	}
336	else return("merge")
337}
338	
339split.or.merge <- function(values, categories, task=c("split","merge"), regimes=FALSE) 
340{	
341# updated to adjust for an additional theta when regimes are split (or one fewer theta if regimes are merged)
342# implement Drummon and Suchard 2010 method, adapting local relaxed clocks for local theta and BMMrates - draw changes in BMMrates from truncated Poisson
343	cc=categories
344
345	if(regimes) vv=values[cc] else vv=values
346	new.vv=vv
347	new.cc=cc
348	nn=length(vv)
349	
350	if(task=="merge") {
351		ncat <- max(cc)
352		m.cc <- sample(1:ncat)[1:2] # grab two categories
353		new.cc[new.cc == m.cc[1]] <- m.cc[2]	# merge selected categories
354		new.cc <- fix.categories(new.cc)
355		or1 <- vv[cc == m.cc[1]][1] # value affected
356		or2 <- vv[cc == m.cc[2]][1]	# value affected
357		n1 <- sum(cc == m.cc[1]) # no. entries affected
358		n2 <- sum(cc == m.cc[2]) # no. entries affected
359		nr <- (n1 * or1 + n2 * or2)/(n1 + n2) # find new value for merged category
360		new.vv[cc == m.cc[1] | cc == m.cc[2]] <- nr # new value assigned to merged classes
361		numSplittableBeforeMerge <- sum(table(cc) > 1)
362		numSplittableAfterMerge <- sum(table(new.cc) > 1)
363		if (max(cc) == nn) {
364			probMerge = 1
365		} else {
366			probMerge = 0.5
367		}
368		
369		if (max(new.cc) == 1) {
370			probSplit = 1
371		} else {
372			probSplit = 0.5
373		}
374		factor = probSplit/probMerge		
375# AS WRITTEN by LJ HARMON
376		lnProposalRatio = log(factor) + log(nn) + log(ncat) - log(numSplittableAfterMerge) - log(2^(n1 + n2) - 2) - log(nr/mean(vv)) - log(n1 + n2)
377		lnPriorRatio = log(Stirling2(nn, ncat)) - log(Stirling2(nn, ncat - 1))
378#	lnHastings = log((2N-2-K+1)/K) # where N is tips, K is number of local parms in tree
379	}
380	if(task=="split") {
381		ncat <- max(cc)
382		while (1) { # find non-unique category
383			rcat <- round(runif(1) * ncat + 0.5) # class to split
384			nc <- sum(cc == rcat) 
385			if (nc > 1) 
386			break
387		}
388		while (1) { # 
389			new <- round(runif(nc) * 2 + 0.5)
390			if (length(table(new)) == 2) 
391			break
392		}
393		new.cc[cc == rcat][new == 2] <- ncat + 1
394		or <- vv[cc == rcat][1] # old rate
395		n1 <- sum(new == 1)	# no. elements affected
396		n2 <- sum(new == 2)	# no. elements affected
397		lims <- c(-0.5 * n1 * or, 0.5 * n2 * or)
398		lims <- lims[order(lims)]
399		u <- runif(1, min = lims[1], max = lims[2])
400		nr1 <- or + u/n1 # new rate 1
401		nr2 <- or - u/n2 # new rate 2
402		new.vv[new.cc == rcat] <- nr1 # assign new rate 1
403		new.vv[new.cc == ncat + 1] <- nr2 # assign new rate 2
404		new.cc <- fix.categories(new.cc)
405		numSplittableBeforeSplit <- sum(table(cc) > 1)
406		numSplittableAfterSplit <- sum(table(new.cc) > 1)
407		if (max(cc) == 1) {
408			probSplit = 1
409		} else {
410			probSplit = 0.5
411		}
412		
413		if (max(new.cc) == nn) {
414			probMerge = 1
415		} else {
416			probMerge = 0.5
417		}
418		
419		factor = probMerge/probSplit
420		
421# AS WRITTEN by LJ HARMON
422		lnProposalRatio = log(factor) + log(numSplittableBeforeSplit) + log(2^(n1 + n2) - 2) + log(or/mean(vv)) + log(n1 + n2) - log(nn) - log(ncat + 1)
423		lnPriorRatio = log(Stirling2(nn, ncat)) - log(Stirling2(nn, ncat + 1))
424#		lnHastings=log((K+1)/(2N-2-K)) # where N is tips, K is number of local parms in tree
425
426	}
427	if(regimes) new.vv=unique(new.vv[order(new.cc)])
428	list(new.values=new.vv, new.categories=new.cc, lnl.prop=lnProposalRatio, lnl.prior=lnPriorRatio)
429}
430
431
432fix.categories <- function(x)
433{
434	f<-factor(x)
435	n<-nlevels(f)
436	o<-numeric(n)
437	for(i in 1:n) o[i]<-which(f==levels(f)[i])[1]
438	nf<-ordered(f, levels(f)[order(o)])	
439	nx<-as.numeric(nf)
440	nx
441}
442
443choose.one <- function(bit.vector)
444{
445	bb=bit.vector
446	s=sample(1:length(bb),1)
447	bb[s]=1-bb[s]
448	bb
449}
450
451fix.heritable.values<-function(bit.vector, values, phy) {
452# assuming most recently changed value in correspondence with bit.vector has been changed on entry
453	bb=bit.vector
454	vv=values
455	names(vv)<-names(bb)<-phy$edge[,2]
456	new.bb=choose.one(bb)
457#	new.vv=vv (do adjust value on the element that was just changed (if 0->1)
458	shifts=sort(as.numeric(names(which(bb>0))))	
459	desc=lapply(shifts, function(x)get.descendants.of.node(x,phy))
460# fix if descendants of one set are included in another
461	new.vv=match()
462}
463
464assign.edgelengths<-function(ouchtree) {
465	ouchtree=as(ouchtree, "data.frame")
466	ouchtree$edges=NA
467	for(i in 2:nrow(ouchtree)) {
468		ouchtree$edges[i]=ouchtree$times[i]-(ouchtree[match(ouchtree[i,'ancestors'], ouchtree[,'nodes']),'times'])
469	}
470	ouchtree$edges[1]=0
471	ouchtree
472}
473
474apply.BMM.to.ouchtree<-function(ouchtree, rates) {
475## FIXME: ensure that the tree maintains structure before and after this function is called
476## VERY SLOW 
477## NOT CURRENTLY implemented
478	ot=as(ouchtree, "data.frame")
479	ot=assign.edgelengths(ot)
480	ot$edges.new=c(0, rates)*ot$edges
481	ot=ot[order(ot$times),]
482	ot$times.new=NA
483	ot$times.new[1]=0
484	for(i in 2:nrow(ot)) {
485		ot$times.new[i]=ot[match(ot[i,'ancestors'], ot[,'nodes']),'times.new']+ot[i,'edges.new']
486	}
487	ot$times=ot$times.new
488	return(with(ot, ouchtree(nodes=nodes, ancestors=ancestors, times=times, labels=labels)))
489}
490
491apply.BMM.to.apetree<-function(apetree, rates) {
492## FASTER than apply.BMM.to.ouchtree
493	apetree$edge.length=apetree$edge.length*rates
494	return(apetree)
495}
496
497adjust.value <- function(value) {
498	rr=value
499	rch <- runif(1, min = -0.5, max = 0.5)
500	rr <- rr + rch
501	rr
502}
503
504adjust.rates <- function(values) {
505	vv=values
506	ncat <- max(fix.categories(vv)->cc)
507	rcat <- round(runif(1) * (ncat) + 0.5)
508	r <- log(vv[cc == rcat][1])
509	rch <- runif(1, min = -2, max = 2)
510	nr <- r + rch
511	vv[cc == rcat] <- exp(nr)
512	vv
513}
514
515prepare.ouchdata <- function(phy, data) {
516	td <- treedata(phy, data, sort = T)	
517	ouch.tre <- ape2ouch(td$phy->ape.tre, scale=FALSE)		
518	ouch.dat=data.frame(as.numeric(ouch.tre@nodes), as.character(ouch.tre@nodelabels))
519	names(ouch.dat)=c("nodes","labels")
520	ouch.dat[ouch.dat[,2]=="",2]=NA
521	ouch.dat$trait=NA
522	mM=match(ouch.dat$labels,names(data))
523	for(i in 1:nrow(ouch.dat)){
524		if(!is.na(mM[i])){
525			ouch.dat$trait[i]=data[mM[i]]	
526		}
527	}
528	return(list(ouch.tre=ouch.tre, ouch.dat=ouch.dat, ape.tre=td$phy, orig.dat=td$data[,1]))
529}
530
531summarize.run<-function(cOK, cProp, currModel) {
532	names(cProp)<-names(cOK)<-c("rates", "rcats", "root", "alpha", "sRegimes", "theta")
533	names(cProp)=paste("prop.",names(cProp),sep="")
534	names(cOK)=paste("okay.",names(cOK),sep="")
535	cat("\n")
536	if(currModel=="BM") {
537		print(cProp[1:3])
538		print(cOK[1:3])
539	} else {
540		print(cProp[c(1,2,4,5,6)])
541		print(cOK[c(1,2,4,5,6)])
542	}
543}
544
545generate.log<-function(bundled.parms, currModel, file, init=FALSE) {
546#	bundled.parms=list(gen=i, mrate=mean(currRates), cats=max(currRates.c), root=currRoot, alpha=currAlpha, reg=max(currRegimes), theta=mean(currTheta), lnL=curr.lnL)
547	res=bundled.parms
548
549	if(init) {
550		if(currModel=="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")
551	} else {
552		if(currModel=="BM") {
553			msg<-paste(res$gen, sQuote(currModel), sprintf("%.3f", res$mrate), res$cats, sprintf("%.3f", res$root), sprintf("%.3f", res$lnL),sep="\t")
554		} else {
555			msg<-paste(res$gen, sQuote(currModel), sprintf("%.3f", res$mrate), res$cats, sprintf("%.4f", res$alpha), res$reg, sprintf("%.3f", res$theta), sprintf("%.3f", res$lnL),sep="\t")
556		}
557	}
558	write(msg, file=file, append=TRUE)
559}
560
561generate.error.message<-function(i, modCurr, modNew, lnR, errorLog, init=FALSE) {
562	if(init) {
563		initmsg = write(paste("gen", "curr.lnL", "new.lnL", "lnR", sep="\t"), file=errorLog)
564	} else {
565		write(paste(i, sprintf("%.3f", modCurr$lnL), sprintf("%.3f", modNew$lnL), sprintf("%.3f", lnR), sep="\t"), file=errorLog)
566	}
567}
568
569tally.mcmc.parms<-function(startparms, endparms, cOK) {
570	index=array(dim=(length(startparms)-1))
571	for(i in 1:(length(startparms)-1)) {
572		index[i]=all(startparms[[i]]==endparms[[i]])
573	}
574	if(!all(index)) {
575		cOK[max(which(index==FALSE))]=cOK[max(which(index==FALSE))]+1
576	}
577	if(length(startparms$theta) == length(endparms$theta)) {
578		if(!all(sort(startparms$theta)==sort(endparms$theta))) {
579			cOK[6]=cOK[6]+1
580		}
581	}
582	cOK
583}
584	
585Stirling2 <- function(n,m)
586{
587    ## Purpose:  Stirling Numbers of the 2-nd kind
588    ## 		S^{(m)}_n = number of ways of partitioning a set of
589    ##                      $n$ elements into $m$ non-empty subsets
590    ## Author: Martin Maechler, Date:  May 28 1992, 23:42
591    ## ----------------------------------------------------------------
592    ## Abramowitz/Stegun: 24,1,4 (p. 824-5 ; Table 24.4, p.835)
593    ## Closed Form : p.824 "C."
594    ## ----------------------------------------------------------------
595
596    if (0 > m || m > n) stop("'m' must be in 0..n !")
597    k <- 0:m
598    sig <- rep(c(1,-1)*(-1)^m, length= m+1)# 1 for m=0; -1 1 (m=1)
599    ## The following gives rounding errors for (25,5) :
600    ## r <- sum( sig * k^n /(gamma(k+1)*gamma(m+1-k)) )
601    ga <- gamma(k+1)
602    round(sum( sig * k^n /(ga * rev(ga))))
603}
604