# function plots posterior density of mapped states from stochastic mapping # written by Liam J. Revell 2012-2013 densityMap<-function(trees,res=100,fsize=NULL,ftype=NULL,lwd=3,check=FALSE,legend=NULL,outline=FALSE,type="phylogram",direction="rightwards",plot=TRUE,...){ if(hasArg(mar)) mar<-list(...)$mar else mar<-rep(0.3,4) tol<-1e-10 if(class(trees)!="multiPhylo") stop("trees not 'multiPhylo' object; just use plotSimmap") h<-sapply(unclass(trees),function(x) max(nodeHeights(x))) steps<-0:res/res*max(h) trees<-rescaleSimmap(trees,totalDepth=max(h)) if(check){ X<-matrix(FALSE,length(trees),length(trees)) for(i in 1:length(trees)) X[i,]<-sapply(trees,all.equal.phylo,current=trees[[i]]) if(!all(X)) stop("some of the trees don't match in topology or relative branch lengths") } tree<-trees[[1]] H<-nodeHeights(tree) message("sorry - this might take a while; please be patient") trees<-unclass(trees) for(i in 1:nrow(tree$edge)){ YY<-cbind(c(H[i,1],steps[intersect(which(steps>H[i,1]),which(stepsH[i,1]),which(steps1){ for(k in 2:length(trees[[j]]$maps[[i]])){ XX[k,1]<-XX[k-1,2] XX[k,2]<-XX[k,1]+trees[[j]]$maps[[i]][k] } } for(k in 1:nrow(YY)){ lower<-which(XX[,1]<=YY[k,1]); lower<-lower[length(lower)] upper<-which(XX[,2]>=(YY[k,2]-tol))[1]; AA<-0 names(lower)<-names(upper)<-NULL for(l in lower:upper) AA<-AA+(min(XX[l,2],YY[k,2])-max(XX[l,1],YY[k,1]))/(YY[k,2]-YY[k,1])*as.numeric(rownames(XX)[l]) ZZ[k]<-ZZ[k]+AA/length(trees) } } tree$maps[[i]]<-YY[,2]-YY[,1] names(tree$maps[[i]])<-round(ZZ*1000) } cols<-rainbow(1001,start=0.7,end=0); names(cols)<-0:1000 tree$mapped.edge<-makeMappedEdge(tree$edge,tree$maps) tree$mapped.edge<-tree$mapped.edge[,order(as.numeric(colnames(tree$mapped.edge)))] x<-list(tree=tree,cols=cols); class(x)<-"densityMap" if(plot) plot.densityMap(x,fsize=fsize,ftype=ftype,lwd=lwd,legend=legend,outline=outline,type=type,mar=mar,direction=direction) invisible(x) } # function # written by Liam J. Revell 2012-2013 plot.densityMap<-function(x,...){ if(class(x)=="densityMap"){ tree<-x$tree cols<-x$cols } else stop("x should be an object of class 'densityMap'") H<-nodeHeights(tree) # get & set optional arguments if(hasArg(legend)) legend<-list(...)$legend else legend<-NULL if(hasArg(fsize)) fsize<-list(...)$fsize else fsize<-NULL if(hasArg(ftype)) ftype<-list(...)$ftype else ftype<-NULL if(hasArg(outline)) outline<-list(...)$outline else outline<-FALSE if(hasArg(lwd)) lwd<-list(...)$lwd else lwd<-3 if(hasArg(leg.txt)) leg.txt<-list(...)$leg.txt else leg.txt<-c("0","PP(state=1)","1") if(hasArg(type)) type<-list(...)$type else type<-"phylogram" if(hasArg(mar)) mar<-list(...)$mar else mar<-rep(0.3,4) if(hasArg(direction)) direction<-list(...)$direction else direction<-"rightwards" if(is.null(legend)) legend<-0.5*max(H) if(is.null(fsize)) fsize<-c(1,1) if(length(fsize)==1) fsize<-rep(fsize,2) if(is.null(ftype)) ftype<-c("i","reg") if(length(ftype)==1) ftype<-c(ftype,"reg") # done optional arguments if(legend){ if(legend>max(H)){ message("legend scale cannot be longer than total tree length; resetting") legend<-0.5*max(H) } } if(type=="phylogram"){ N<-length(tree$tip.label) if(legend) ylim<-c(1-0.12*(N-1),N) else ylim<-NULL if(outline){ par(col="white") plotTree(tree,fsize=fsize[1],lwd=lwd+2,offset=0.2*lwd/3+0.2/3,ftype=ftype[1],ylim=ylim,mar=mar,direction=direction) par(col="black") add<-TRUE } else add<-FALSE plotSimmap(tree,cols,pts=FALSE,lwd=lwd,fsize=fsize[1],mar=mar,ftype=ftype[1],add=add,ylim=ylim,direction=direction) if(legend){ ff<-function(dd){ if(!("."%in%dd)) dig<-0 else dig<-length(dd)-which(dd==".") dig } dig<-max(sapply(strsplit(leg.txt[c(1,3)],split=""),ff)) add.color.bar(legend,cols,title=leg.txt[2],lims<-as.numeric(leg.txt[c(1,3)]),digits=dig,prompt=FALSE,x=0,y=1-0.08*(N-1),lwd=lwd,fsize=fsize[2]) } } else if(type=="fan"){ if(outline){ par(col="white") invisible(capture.output(plotTree(tree,type="fan",lwd=lwd+2,mar=mar,fsize=fsize[1],ftype=ftype[1]))) par(col="black") } invisible(capture.output(plotSimmap(tree,cols,lwd=lwd,mar=mar,fsize=fsize[1],add=outline,ftype=ftype[1],type="fan"))) if(legend){ ff<-function(dd){ if(!("."%in%dd)) dig<-0 else dig<-length(dd)-which(dd==".") dig } dig<-max(sapply(strsplit(leg.txt[c(1,3)],split=""),ff)) add.color.bar(legend,cols,title=leg.txt[2],lims<-as.numeric(leg.txt[c(1,3)]),digits=dig,prompt=FALSE,x=0.9*par()$usr[1],y=0.9*par()$usr[3],lwd=lwd,fsize=fsize[2]) } } }