3 ## Nonparametric Rate Smoothing Method by Sanderson
5 ## Copyright 2003 Gangolf Jobb and Korbinian Strimmer
7 ## This file is part of the R-package `ape'.
8 ## See the file ../COPYING for licensing issues.
11 function(lowerNodes,upperNodes,edgeLengths,minEdgeLength,tipLabels)
14 as.integer(lowerNodes),
15 as.integer(upperNodes),
16 as.double(edgeLengths),
17 as.double(minEdgeLength),
18 as.integer(length(edgeLengths)),
19 as.character(tipLabels),
20 as.integer(length(tipLabels)),
45 result=double(getNEdges()),
60 function(params,scale)
65 result=double(getNEdges()),
70 function(params,scale)
75 result=double(getNEdges()),
83 result=double(getNFreeParams()),
89 prepareTree <- function(phy, minEdgeLength = 1e-06)
91 len <- phy$edge.length
92 if (length(len) > 2048) stop("Only 2048 branches in tree allowed!")
93 low <- phy$edge[, 1] # edges in the tree
95 setTree(low, upp, len, minEdgeLength, phy$tip.labels)
98 optimTree <- function(phy, expo = 2) # call prepareTree first
100 dur <- rep(log(0.5), getNFreeParams() ) # start value
101 objL <- function(d) objFuncLogScale(d, expo)
102 opt <- optim(dur, objL, method = "BFGS")
106 ### this is just for testing purposes, to get the tree we are
107 ### actually using when there are many small branch lengths
108 phylogram <- function(phy, ...)
110 if (class(phy) != "phylo") stop("object \"phy\" is not of class \"phylo\"")
112 ## added by EP for the new coding of "phylo" (2006-10-04):
113 phy <- new2old.phylo(phy)
116 prepareTree(phy, ...)
117 ##opt <- optimTree(phy, ...)
120 newTree$edge.length <- getEdgeLengths()
128 chronogram <- function(phy, scale = 1, expo = 2, minEdgeLength = 1e-06)
130 if (class(phy) != "phylo") stop("object \"phy\" is not of class \"phylo\"")
132 ## added by EP for the new coding of "phylo" (2006-10-04):
133 phy <- new2old.phylo(phy)
136 prepareTree(phy, minEdgeLength = minEdgeLength)
137 opt <- optimTree(phy, expo = expo)
140 newTree$edge.length <- getDurations(opt$par, scale)
146 ratogram <- function(phy, scale = 1, expo = 2, minEdgeLength = 1e-06)
148 if (class(phy) != "phylo")
149 stop("object \"phy\" is not of class \"phylo\"")
151 ## added by EP for the new coding of "phylo" (2006-10-04):
152 phy <- new2old.phylo(phy)
155 prepareTree(phy, minEdgeLength = minEdgeLength)
156 opt <- optimTree(phy, expo = expo)
159 newTree$edge.length <- getRates(opt$par, scale)
165 NPRS.criterion <- function(phy, chrono, expo = 2, minEdgeLength = 1e-06)
167 if (!is.ultrametric(chrono))
168 stop("tree \"chrono\" is not ultrametric (clock-like)")
170 ## added by EP for the new coding of "phylo" (2006-10-04):
171 phy <- new2old.phylo(phy)
172 chrono <- new2old.phylo(chrono)
175 prepareTree(chrono, minEdgeLength = minEdgeLength)
176 parms <- getExternalParams()
177 prepareTree(phy, minEdgeLength = minEdgeLength)
178 objFuncLogScale(parms, expo)