1 ## plot.phylo.R (2013-01-11)
5 ## Copyright 2002-2013 Emmanuel Paradis
7 ## This file is part of the R-package `ape'.
8 ## See the file ../COPYING for licensing issues.
11 function(x, type = "phylogram", use.edge.length = TRUE,
12 node.pos = NULL, show.tip.label = TRUE,
13 show.node.label = FALSE, edge.color = "black",
14 edge.width = 1, edge.lty = 1, font = 3, cex = par("cex"),
15 adj = NULL, srt = 0, no.margin = FALSE, root.edge = FALSE,
16 label.offset = 0, underscore = FALSE, x.lim = NULL,
17 y.lim = NULL, direction = "rightwards", lab4ut = "horizontal",
18 tip.color = "black", plot = TRUE, rotate.tree = 0,
21 Ntip <- length(x$tip.label)
23 warning("found less than 2 tips in the tree")
26 if (any(tabulate(x$edge[, 1]) == 1))
27 stop("there are single (non-splitting) nodes in your tree; you may need to use collapse.singles()")
29 .nodeHeight <- function(Ntip, Nnode, edge, Nedge, yy)
30 .C("node_height", as.integer(Ntip), as.integer(Nnode),
31 as.integer(edge[, 1]), as.integer(edge[, 2]),
32 as.integer(Nedge), as.double(yy),
33 DUP = FALSE, PACKAGE = "ape")[[6]]
35 .nodeDepth <- function(Ntip, Nnode, edge, Nedge)
36 .C("node_depth", as.integer(Ntip), as.integer(Nnode),
37 as.integer(edge[, 1]), as.integer(edge[, 2]),
38 as.integer(Nedge), double(Ntip + Nnode),
39 DUP = FALSE, PACKAGE = "ape")[[6]]
41 .nodeDepthEdgelength <- function(Ntip, Nnode, edge, Nedge, edge.length)
42 .C("node_depth_edgelength", as.integer(Ntip),
43 as.integer(Nnode), as.integer(edge[, 1]),
44 as.integer(edge[, 2]), as.integer(Nedge),
45 as.double(edge.length), double(Ntip + Nnode),
46 DUP = FALSE, PACKAGE = "ape")[[7]]
48 Nedge <- dim(x$edge)[1]
51 type <- match.arg(type, c("phylogram", "cladogram", "fan",
52 "unrooted", "radial"))
53 direction <- match.arg(direction, c("rightwards", "leftwards",
54 "upwards", "downwards"))
55 if (is.null(x$edge.length)) use.edge.length <- FALSE
57 ## the order of the last two conditions is important:
58 if (type %in% c("unrooted", "radial") || !use.edge.length ||
59 is.null(x$root.edge) || !x$root.edge) root.edge <- FALSE
60 if (type == "fan" && root.edge) {
61 warning("drawing root edge with type = 'fan' is not yet supported")
65 phyloORclado <- type %in% c("phylogram", "cladogram")
66 horizontal <- direction %in% c("rightwards", "leftwards")
67 xe <- x$edge # to save
69 ## we first compute the y-coordinates of the tips.
70 phyOrder <- attr(x, "order")
71 ## make sure the tree is in cladewise order:
72 if (is.null(phyOrder) || phyOrder != "cladewise") {
73 x <- reorder(x) # fix from Klaus Schliep (2007-06-16)
74 if (!identical(x$edge, xe)) {
75 ## modified from Li-San Wang's fix (2007-01-23):
76 ereorder <- match(x$edge[, 2], xe[, 2])
77 if (length(edge.color) > 1) {
78 edge.color <- rep(edge.color, length.out = Nedge)
79 edge.color <- edge.color[ereorder]
81 if (length(edge.width) > 1) {
82 edge.width <- rep(edge.width, length.out = Nedge)
83 edge.width <- edge.width[ereorder]
85 if (length(edge.lty) > 1) {
86 edge.lty <- rep(edge.lty, length.out = Nedge)
87 edge.lty <- edge.lty[ereorder]
91 ### By contrats to ape (< 2.4), the arguments edge.color, etc., are
92 ### not elongated before being passed to segments(), except if needed
94 yy <- numeric(Ntip + Nnode)
95 TIPS <- x$edge[x$edge[, 2] <= Ntip, 2]
98 ## 'z' is the tree in pruningwise order used in calls to .C
99 z <- reorder(x, order = "pruningwise")
102 if (is.null(node.pos)) {
104 if (type == "cladogram" && !use.edge.length) node.pos <- 2
107 yy <- .nodeHeight(Ntip, Nnode, z$edge, Nedge, yy)
109 ## node_height_clado requires the number of descendants
110 ## for each node, so we compute `xx' at the same time
111 ans <- .C("node_height_clado", as.integer(Ntip),
112 as.integer(Nnode), as.integer(z$edge[, 1]),
113 as.integer(z$edge[, 2]), as.integer(Nedge),
114 double(Ntip + Nnode), as.double(yy),
115 DUP = FALSE, PACKAGE = "ape")
119 if (!use.edge.length) {
120 if (node.pos != 2) xx <- .nodeDepth(Ntip, Nnode, z$edge, Nedge) - 1
123 xx <- .nodeDepthEdgelength(Ntip, Nnode, z$edge, Nedge, z$edge.length)
127 rotate.tree <- twopi * rotate.tree/360
128 switch(type, "fan" = {
129 ## if the tips are not in the same order in tip.label
130 ## and in edge[, 2], we must reorder the angles: we
131 ## use `xx' to store temporarily the angles
132 TIPS <- x$edge[which(x$edge[, 2] <= Ntip), 2]
133 xx <- seq(0, twopi * (1 - 1/Ntip) - twopi * open.angle/360,
135 theta <- double(Ntip)
137 theta <- c(theta, numeric(Nnode))
138 theta <- .nodeHeight(Ntip, Nnode, z$edge, Nedge, theta)
139 if (use.edge.length) {
140 r <- .nodeDepthEdgelength(Ntip, Nnode, z$edge, Nedge, z$edge.length)
142 r <- .nodeDepth(Ntip, Nnode, z$edge, Nedge)
145 theta <- theta + rotate.tree
149 nb.sp <- .nodeDepth(Ntip, Nnode, z$edge, Nedge)
150 XY <- if (use.edge.length)
151 unrooted.xy(Ntip, Nnode, z$edge, z$edge.length, nb.sp, rotate.tree)
153 unrooted.xy(Ntip, Nnode, z$edge, rep(1, Nedge), nb.sp, rotate.tree)
154 ## rescale so that we have only positive values
155 xx <- XY$M[, 1] - min(XY$M[, 1])
156 yy <- XY$M[, 2] - min(XY$M[, 2])
158 X <- .nodeDepth(Ntip, Nnode, z$edge, Nedge)
162 ## angle (1st compute the angles for the tips):
163 yy <- c((1:Ntip)*twopi/Ntip, rep(0, Nnode))
164 Y <- .nodeHeight(Ntip, Nnode, z$edge, Nedge, yy)
165 xx <- X * cos(Y + rotate.tree)
166 yy <- X * sin(Y + rotate.tree)
172 xx <- tmp - min(tmp) + 1
175 if (direction == "rightwards") xx <- xx + x$root.edge
176 if (direction == "upwards") yy <- yy + x$root.edge
179 if (no.margin) par(mai = rep(0, 4))
180 if (is.null(x.lim)) {
184 pin1 <- par("pin")[1] # width of the device in inches
185 strWi <- strwidth(x$tip.label, "inches") # id. for the tip labels
186 ## 1.04 comes from that we are using a regular axis system
187 ## with 4% on both sides of the range of x:
188 xx.tips <- xx[1:Ntip] * 1.04
189 ## 'alp' is the conversion coefficient from
190 ## user coordinates to inches:
191 alp <- try(uniroot(function(a) max(a*xx.tips + strWi) - pin1,
192 c(0, 1e6))$root, silent = TRUE)
193 ## if the above fails, give 1/3 of the device for the tip labels:
194 if (is.character(alp)) tmp <- max(xx.tips)*1.5 else {
195 tmp <- if (show.tip.label) max(xx.tips + strWi/alp) else max(xx.tips)
198 } else x.lim <- c(1, Ntip)
199 } else switch(type, "fan" = {
200 if (show.tip.label) {
201 offset <- max(nchar(x$tip.label) * 0.018 * max(yy) * cex)
202 x.lim <- c(min(xx) - offset, max(xx) + offset)
203 } else x.lim <- c(min(xx), max(xx))
205 if (show.tip.label) {
206 offset <- max(nchar(x$tip.label) * 0.018 * max(yy) * cex)
207 x.lim <- c(0 - offset, max(xx) + offset)
208 } else x.lim <- c(0, max(xx))
210 if (show.tip.label) {
211 offset <- max(nchar(x$tip.label) * 0.03 * cex)
212 x.lim <- c(-1 - offset, 1 + offset)
213 } else x.lim <- c(-1, 1)
215 } else if (length(x.lim) == 1) {
217 if (phyloORclado && !horizontal) x.lim[1] <- 1
218 if (type %in% c("fan", "unrooted") && show.tip.label)
219 x.lim[1] <- -max(nchar(x$tip.label) * 0.018 * max(yy) * cex)
220 if (type == "radial")
222 if (show.tip.label) -1 - max(nchar(x$tip.label) * 0.03 * cex)
226 if (phyloORclado && direction == "leftwards") xx <- x.lim[2] - xx
227 if (is.null(y.lim)) {
229 if (horizontal) y.lim <- c(1, Ntip) else {
231 pin2 <- par("pin")[2] # height of the device in inches
232 strWi <- strwidth(x$tip.label, "inches")
233 ## 1.04 comes from that we are using a regular axis system
234 ## with 4% on both sides of the range of x:
235 yy.tips <- yy[1:Ntip] * 1.04
236 ## 'alp' is the conversion coefficient from
237 ## user coordinates to inches:
238 alp <- try(uniroot(function(a) max(a*yy.tips + strWi) - pin2,
239 c(0, 1e6))$root, silent = TRUE)
240 ## if the above fails, give 1/3 of the device for the tip labels:
241 if (is.character(alp)) tmp <- max(yy.tips)*1.5 else {
242 tmp <- if (show.tip.label) max(yy.tips + strWi/alp) else max(yy.tips)
246 } else switch(type, "fan" = {
247 if (show.tip.label) {
248 offset <- max(nchar(x$tip.label) * 0.018 * max(yy) * cex)
249 y.lim <- c(min(yy) - offset, max(yy) + offset)
250 } else y.lim <- c(min(yy), max(yy))
252 if (show.tip.label) {
253 offset <- max(nchar(x$tip.label) * 0.018 * max(yy) * cex)
254 y.lim <- c(0 - offset, max(yy) + offset)
255 } else y.lim <- c(0, max(yy))
257 if (show.tip.label) {
258 offset <- max(nchar(x$tip.label) * 0.03 * cex)
259 y.lim <- c(-1 - offset, 1 + offset)
260 } else y.lim <- c(-1, 1)
262 } else if (length(y.lim) == 1) {
264 if (phyloORclado && horizontal) y.lim[1] <- 1
265 if (type %in% c("fan", "unrooted") && show.tip.label)
266 y.lim[1] <- -max(nchar(x$tip.label) * 0.018 * max(yy) * cex)
267 if (type == "radial")
268 y.lim[1] <- if (show.tip.label) -1 - max(nchar(x$tip.label) * 0.018 * max(yy) * cex) else -1
271 if (phyloORclado && direction == "downwards") yy <- max(yy) - yy
272 if (phyloORclado && root.edge) {
273 if (direction == "leftwards") x.lim[2] <- x.lim[2] + x$root.edge
274 if (direction == "downwards") y.lim[2] <- y.lim[2] + x$root.edge
276 asp <- if (type %in% c("fan", "radial", "unrooted")) 1 else NA # fixes by Klaus Schliep (2008-03-28 and 2010-08-12)
277 plot(0, type = "n", xlim = x.lim, ylim = y.lim, ann = FALSE, axes = FALSE, asp = asp, ...)
281 adj <- if (phyloORclado && direction == "leftwards") 1 else 0
282 if (phyloORclado && show.tip.label) {
283 MAXSTRING <- max(strwidth(x$tip.label, cex = cex))
285 if (direction == "rightwards") {
286 lox <- label.offset + MAXSTRING * 1.05 * adj
288 if (direction == "leftwards") {
289 lox <- -label.offset - MAXSTRING * 1.05 * (1 - adj)
290 ##xx <- xx + MAXSTRING
294 MAXSTRING <- MAXSTRING * 1.09 * (psr[4] - psr[3])/(psr[2] - psr[1])
295 loy <- label.offset + MAXSTRING * 1.05 * adj
298 if (direction == "downwards") {
300 ##yy <- yy + MAXSTRING
305 if (type == "phylogram") {
306 phylogram.plot(x$edge, Ntip, Nnode, xx, yy,
307 horizontal, edge.color, edge.width, edge.lty)
310 ereorder <- match(z$edge[, 2], x$edge[, 2])
311 if (length(edge.color) > 1) {
312 edge.color <- rep(edge.color, length.out = Nedge)
313 edge.color <- edge.color[ereorder]
315 if (length(edge.width) > 1) {
316 edge.width <- rep(edge.width, length.out = Nedge)
317 edge.width <- edge.width[ereorder]
319 if (length(edge.lty) > 1) {
320 edge.lty <- rep(edge.lty, length.out = Nedge)
321 edge.lty <- edge.lty[ereorder]
323 circular.plot(z$edge, Ntip, Nnode, xx, yy, theta,
324 r, edge.color, edge.width, edge.lty)
326 cladogram.plot(x$edge, xx, yy, edge.color, edge.width, edge.lty)
330 "rightwards" = segments(0, yy[ROOT], x$root.edge, yy[ROOT]),
331 "leftwards" = segments(xx[ROOT], yy[ROOT], xx[ROOT] + x$root.edge, yy[ROOT]),
332 "upwards" = segments(xx[ROOT], 0, xx[ROOT], x$root.edge),
333 "downwards" = segments(xx[ROOT], yy[ROOT], xx[ROOT], yy[ROOT] + x$root.edge))
334 if (show.tip.label) {
335 if (is.expression(x$tip.label)) underscore <- TRUE
336 if (!underscore) x$tip.label <- gsub("_", " ", x$tip.label)
339 text(xx[1:Ntip] + lox, yy[1:Ntip] + loy, x$tip.label, adj = adj,
340 font = font, srt = srt, cex = cex, col = tip.color)
342 if (type == "unrooted") {
343 if (lab4ut == "horizontal") {
344 y.adj <- x.adj <- numeric(Ntip)
345 sel <- abs(XY$axe) > 0.75 * pi
346 x.adj[sel] <- -strwidth(x$tip.label)[sel] * 1.05
347 sel <- abs(XY$axe) > pi/4 & abs(XY$axe) < 0.75 * pi
348 x.adj[sel] <- -strwidth(x$tip.label)[sel] * (2 * abs(XY$axe)[sel] / pi - 0.5)
349 sel <- XY$axe > pi / 4 & XY$axe < 0.75 * pi
350 y.adj[sel] <- strheight(x$tip.label)[sel] / 2
351 sel <- XY$axe < -pi / 4 & XY$axe > -0.75 * pi
352 y.adj[sel] <- -strheight(x$tip.label)[sel] * 0.75
353 text(xx[1:Ntip] + x.adj * cex, yy[1:Ntip] + y.adj * cex,
354 x$tip.label, adj = c(adj, 0), font = font,
355 srt = srt, cex = cex, col = tip.color)
356 } else { # if lab4ut == "axial"
357 adj <- abs(XY$axe) > pi/2
358 srt <- 180 * XY$axe / pi
359 srt[adj] <- srt[adj] - 180
360 adj <- as.numeric(adj)
361 xx.tips <- xx[1:Ntip]
362 yy.tips <- yy[1:Ntip]
364 xx.tips <- xx.tips + label.offset * cos(XY$axe)
365 yy.tips <- yy.tips + label.offset * sin(XY$axe)
367 ## `srt' takes only a single value, so can't vectorize this:
368 ## (and need to 'elongate' these vectors:)
369 font <- rep(font, length.out = Ntip)
370 tip.color <- rep(tip.color, length.out = Ntip)
371 cex <- rep(cex, length.out = Ntip)
373 text(xx.tips[i], yy.tips[i], cex = cex[i],
374 x$tip.label[i], adj = adj[i], font = font[i],
375 srt = srt[i], col = tip.color[i])
378 if (type %in% c("fan", "radial")) {
379 xx.tips <- xx[1:Ntip]
380 yy.tips <- yy[1:Ntip]
381 angle <- atan2(yy.tips, xx.tips) # in radians
383 xx.tips <- xx.tips + label.offset * cos(angle)
384 yy.tips <- yy.tips + label.offset * sin(angle)
387 angle <- angle * 180/pi # switch to degrees
388 angle[s] <- angle[s] + 180
390 ## `srt' takes only a single value, so can't vectorize this:
391 ## (and need to 'elongate' these vectors:)
392 font <- rep(font, length.out = Ntip)
393 tip.color <- rep(tip.color, length.out = Ntip)
394 cex <- rep(cex, length.out = Ntip)
396 text(xx.tips[i], yy.tips[i], x$tip.label[i], font = font[i],
397 cex = cex[i], srt = angle[i], adj = adj[i],
402 text(xx[ROOT:length(xx)] + label.offset, yy[ROOT:length(yy)],
403 x$node.label, adj = adj, font = font, srt = srt, cex = cex)
405 L <- list(type = type, use.edge.length = use.edge.length,
406 node.pos = node.pos, show.tip.label = show.tip.label,
407 show.node.label = show.node.label, font = font,
408 cex = cex, adj = adj, srt = srt, no.margin = no.margin,
409 label.offset = label.offset, x.lim = x.lim, y.lim = y.lim,
410 direction = direction, tip.color = tip.color,
411 Ntip = Ntip, Nnode = Nnode)
412 assign("last_plot.phylo", c(L, list(edge = xe, xx = xx, yy = yy)),
413 envir = .PlotPhyloEnv)
417 phylogram.plot <- function(edge, Ntip, Nnode, xx, yy, horizontal,
418 edge.color, edge.width, edge.lty)
420 nodes <- (Ntip + 1):(Ntip + Nnode)
426 ## un trait vertical a chaque noeud...
428 y0v <- y1v <- numeric(Nnode)
429 ## store the index of each node in the 1st column of edge:
430 NodeInEdge1 <- vector("list", Nnode)
433 j <- NodeInEdge1[[ii]] <- which(edge[, 1] == i)
434 tmp <- range(yy[edge[j, 2]])
438 ## ... et un trait horizontal partant de chaque tip et chaque noeud
444 nc <- length(edge.color)
445 nw <- length(edge.width)
446 nl <- length(edge.lty)
448 if (nc + nw + nl == 3) {
449 color.v <- edge.color
450 width.v <- edge.width
453 Nedge <- dim(edge)[1]
454 edge.color <- rep(edge.color, length.out = Nedge)
455 edge.width <- rep(edge.width, length.out = Nedge)
456 edge.lty <- rep(edge.lty, length.out = Nedge)
457 DF <- data.frame(edge.color, edge.width, edge.lty, stringsAsFactors = FALSE)
458 color.v <- rep("black", Nnode)
459 width.v <- rep(1, Nnode)
460 lty.v <- rep(1, Nnode)
462 br <- NodeInEdge1[[i]]
463 if (length(br) > 2) {
464 x <- unique(DF[br, 1])
465 if (length(x) == 1) color.v[i] <- x
466 x <- unique(DF[br, 2])
467 if (length(x) == 1) width.v[i] <- x
468 x <- unique(DF[br, 3])
469 if (length(x) == 1) lty.v[i] <- x
473 if (any(DF[A, ] != DF[B, ])) {
474 color.v[i] <- edge.color[B]
475 width.v[i] <- edge.width[B]
476 lty.v[i] <- edge.lty[B]
478 y0v <- c(y0v, y0v[i])
479 y1v <- c(y1v, yy[i + Ntip])
480 x0v <- c(x0v, x0v[i])
481 color.v <- c(color.v, edge.color[A])
482 width.v <- c(width.v, edge.width[A])
483 lty.v <- c(lty.v, edge.lty[A])
485 y0v[i] <- yy[i + Ntip]
487 color.v[i] <- edge.color[A]
488 width.v[i] <- edge.width[A]
489 lty.v[i] <- edge.lty[A]
496 segments(x0h, y0h, x1h, y0h, col = edge.color, lwd = edge.width, lty = edge.lty) # draws horizontal lines
497 segments(x0v, y0v, x0v, y1v, col = color.v, lwd = width.v, lty = lty.v) # draws vertical lines
499 segments(y0h, x0h, y0h, x1h, col = edge.color, lwd = edge.width, lty = edge.lty) # draws vertical lines
500 segments(y0v, x0v, y1v, x0v, col = color.v, lwd = width.v, lty = lty.v) # draws horizontal lines
504 cladogram.plot <- function(edge, xx, yy, edge.color, edge.width, edge.lty)
505 segments(xx[edge[, 1]], yy[edge[, 1]], xx[edge[, 2]], yy[edge[, 2]],
506 col = edge.color, lwd = edge.width, lty = edge.lty)
508 circular.plot <- function(edge, Ntip, Nnode, xx, yy, theta,
509 r, edge.color, edge.width, edge.lty)
510 ### 'edge' must be in pruningwise order
514 theta0 <- theta[edge[, 2]]
515 costheta0 <- cos(theta0)
516 sintheta0 <- sin(theta0)
523 segments(x0, y0, x1, y1, col = edge.color, lwd = edge.width, lty = edge.lty)
525 tmp <- which(diff(edge[, 1]) != 0)
526 start <- c(1, tmp + 1)
527 Nedge <- dim(edge)[1]
530 ## function dispatching the features to the arcs
531 foo <- function(edge.feat, default) {
532 if (length(edge.feat) == 1) return(rep(edge.feat, Nnode))
534 edge.feat <- rep(edge.feat, length.out = Nedge)
535 feat.arc <- rep(default, Nnode)
537 tmp <- edge.feat[start[k]]
538 if (tmp == edge.feat[end[k]]) feat.arc[k] <- tmp
543 co <- foo(edge.color, "black")
544 lw <- foo(edge.width, 1)
545 ly <- foo(edge.lty, 1)
550 X <- rep(r[edge[i, 1]], 100)
551 Y <- seq(theta[edge[i, 2]], theta[edge[j, 2]], length.out = 100)
552 lines(X*cos(Y), X*sin(Y), col = co[k], lwd = lw[k], lty = ly[k])
556 unrooted.xy <- function(Ntip, Nnode, edge, edge.length, nb.sp, rotate.tree)
558 foo <- function(node, ANGLE, AXIS) {
559 ind <- which(edge[, 1] == node)
561 start <- AXIS - ANGLE/2
562 for (i in 1:length(sons)) {
563 h <- edge.length[ind[i]]
564 angle[sons[i]] <<- alpha <- ANGLE*nb.sp[sons[i]]/nb.sp[node]
565 axis[sons[i]] <<- beta <- start + alpha/2
566 start <- start + alpha
567 xx[sons[i]] <<- h*cos(beta) + xx[node]
568 yy[sons[i]] <<- h*sin(beta) + yy[node]
571 if (i > Ntip) foo(i, angle[i], axis[i])
573 Nedge <- dim(edge)[1]
574 yy <- xx <- numeric(Ntip + Nnode)
575 ## `angle': the angle allocated to each node wrt their nb of tips
576 ## `axis': the axis of each branch
577 axis <- angle <- numeric(Ntip + Nnode)
578 ## start with the root...
579 foo(Ntip + 1L, 2*pi, 0 + rotate.tree)
582 axe <- axis[1:Ntip] # the axis of the terminal branches (for export)
584 ## insures that returned angles are in [-PI, +PI]:
585 axe[axeGTpi] <- axe[axeGTpi] - 2*pi
586 list(M = M, axe = axe)
589 node.depth <- function(phy)
591 n <- length(phy$tip.label)
593 N <- dim(phy$edge)[1]
594 phy <- reorder(phy, order = "pruningwise")
595 .C("node_depth", as.integer(n), as.integer(m),
596 as.integer(phy$edge[, 1]), as.integer(phy$edge[, 2]),
597 as.integer(N), double(n + m), DUP = FALSE, PACKAGE = "ape")[[6]]
600 node.depth.edgelength <- function(phy)
602 n <- length(phy$tip.label)
604 N <- dim(phy$edge)[1]
605 phy <- reorder(phy, order = "pruningwise")
606 .C("node_depth_edgelength", as.integer(n), as.integer(n),
607 as.integer(phy$edge[, 1]), as.integer(phy$edge[, 2]),
608 as.integer(N), as.double(phy$edge.length), double(n + m),
609 DUP = FALSE, PACKAGE = "ape")[[7]]
612 node.height <- function(phy)
614 n <- length(phy$tip.label)
616 N <- dim(phy$edge)[1]
617 phy <- reorder(phy, order = "pruningwise")
626 .C("node_height", as.integer(n), as.integer(m),
627 as.integer(e1), as.integer(e2), as.integer(N),
628 as.double(yy), DUP = FALSE, PACKAGE = "ape")[[6]]
631 node.height.clado <- function(phy)
633 n <- length(phy$tip.label)
635 N <- dim(phy$edge)[1]
636 phy <- reorder(phy, order = "pruningwise")
645 .C("node_height_clado", as.integer(n), as.integer(m),
646 as.integer(e1), as.integer(e2), as.integer(N),
647 double(n + m), as.double(yy), DUP = FALSE,
648 PACKAGE = "ape")[[7]]
651 plot.multiPhylo <- function(x, layout = 1, ...)
653 layout(matrix(1:layout, ceiling(sqrt(layout)), byrow = TRUE))
654 if (!devAskNewPage() && interactive()) {
656 on.exit(devAskNewPage(FALSE))
658 for (i in 1:length(x)) plot(x[[i]], ...)
661 trex <- function(phy, title = TRUE, subbg = "lightyellow3",
662 return.tree = FALSE, ...)
664 lastPP <- get("last_plot.phylo", envir = .PlotPhyloEnv)
665 devmain <- dev.cur() # where the main tree is plotted
667 restore <- function() {
669 assign("last_plot.phylo", lastPP, envir = .PlotPhyloEnv)
674 cat("Click close to a node. Right-click to exit.\n")
676 x <- identify.phylo(phy, quiet = TRUE)
677 if (is.null(x)) return(invisible(NULL)) else {
679 if (is.null(x)) cat("Try again!\n") else {
685 } else dev.set(devsub)
687 tr <- extract.clade(phy, x)
689 if (is.character(title)) title(title)
692 if (is.null(phy$node.label))
693 paste("From node #", x, sep = "")
694 else paste("From", phy$node.label[x - Ntip(phy)])
697 if (return.tree) return(tr)
704 kronoviz <- function(x, layout = length(x), horiz = TRUE, ...)
706 par(mar = rep(0.5, 4), oma = rep(2, 4))
707 rts <- sapply(x, function(x) branching.times(x)[1])
709 lim <- cbind(rts - maxrts, rts)
711 Ntips <- sapply(x, Ntip)
721 layout(matrix(1:layout, nrow), widths = w, heights = h)
722 if (layout < Ntree && !devAskNewPage() && interactive()) {
724 on.exit(devAskNewPage(FALSE))
728 plot(x[[i]], x.lim = lim[i, ], ...)
731 plot(x[[i]], y.lim = lim[i, ], direction = "u", ...)
733 axisPhylo(if (horiz) 1 else 4) # better if the deepest tree is last ;)