/working/auteur.extended/R/generate.starting.point.R

 1#utility for generating a list of indicator variables and values for each branch in a phylogeny, using a random model complexity as a starting point (or by constraining model complexity)
 2#author: JM EASTMAN 2010
 3
 4generate.starting.point <-
 5function(data, phy, node.des, model, K=FALSE, jumpsize) { 
 6	nn=length(phy$edge.length)
 7	ntip=Ntip(phy)
 8	if(!K) nshifts=rtpois(1,log(ntip),nn) else nshifts=K-1
 9	if(nshifts!=0) bb=sample(c(rep(1,nshifts),rep(0,nn-nshifts)),replace=FALSE) else bb=rep(0,nn)
10	names(bb)=phy$edge[,2]
11	shifts=as.numeric(names(bb)[bb==1])
12	if(model=="OU") {
13		min.max=c(adjustvalue(min(data), jumpsize), adjustvalue(max(data), jumpsize))
14		values=runif(sum(bb)+1, min=min.max[min.max==min(min.max)], max=min.max[min.max==max(min.max)])
15	} else if(model=="BM") {
16		init.rate=fit.continuous(phy,data)
17		min.max=c(adjustrate(init.rate, jumpsize), adjustrate(init.rate, jumpsize))
18		values=runif(sum(bb)+1, min=min.max[min.max==min(min.max)], max=min.max[min.max==max(min.max)])
19	}
20	internal.shifts<-tip.shifts<-numeric(0)
21	internal.shifts=sort(shifts[shifts>ntip])
22	tip.shifts=shifts[shifts<=ntip]
23	
24	if(length(internal.shifts)==0 & length(tip.shifts)==0) {
25		vv=rep(values, nn)
26	} else {
27		# assign randomly chosen values for each branch, respecting phylogenetic relationship
28		vv=bb
29		vv[]=values[length(values)]
30		i=0
31		if(length(internal.shifts)!=0) {
32			for(i in 1:length(internal.shifts)) {
33				d=node.des[which(names(node.des)==internal.shifts[i])]
34				vv[match(c(internal.shifts[i], unlist(d)), names(vv))]=values[i]
35			}
36		}
37		if(length(tip.shifts)!=0) {
38			for(j in 1:length(tip.shifts)) {
39				vv[match(tip.shifts[j], names(vv))]=values[j+i]
40			}
41		}
42	}
43	return(list(delta=unname(bb), values=unname(vv)))
44}
45