From 57c9d4fd10a4449c27c65873eaf9315b5b7ed874 Mon Sep 17 00:00:00 2001 From: paradis Date: Tue, 14 Dec 2010 15:51:09 +0000 Subject: [PATCH] several fixes git-svn-id: https://svn.mpl.ird.fr/ape/dev/ape@142 6e262413-ae40-0410-9e79-b911bd7a66b7 --- ChangeLog | 14 +- DESCRIPTION | 4 +- R/ace.R | 5 +- R/as.phylo.R | 3 +- R/compar.gee.R | 2 +- R/plot.phylo.R | 483 ++++++++++++++++++++++++++++--------------------- Thanks | 11 +- 7 files changed, 303 insertions(+), 219 deletions(-) diff --git a/ChangeLog b/ChangeLog index 057cf58..53f4a62 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,12 @@ + CHANGES IN APE VERSION 2.6-3 + + +BUG FIXES + + o as.hclust.phylo() now returns an error with unrooted trees. + + + CHANGES IN APE VERSION 2.6-2 @@ -32,7 +41,10 @@ BUG FIXES cutree() or rect.hclust(). o write.tree() did not output correctly node labels (thanks to Naim - Matasci for the fix). + Matasci and Jeremy Beaulieu for the fix). + + o ace(type = "discrete") has been improved thanks to Naim Marasci and + Jeremy Beaulieu. diff --git a/DESCRIPTION b/DESCRIPTION index 8fd651a..a9d4511 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: ape -Version: 2.6-2 -Date: 2010-12-07 +Version: 2.6-3 +Date: 2010-12-14 Title: Analyses of Phylogenetics and Evolution Author: Emmanuel Paradis, Ben Bolker, Julien Claude, Hoa Sien Cuong, Richard Desper, Benoit Durand, Julien Dutheil, Olivier Gascuel, Christoph Heibl, Daniel Lawson, Vincent Lefort, Pierre Legendre, Jim Lemon, Yvonnick Noel, Johan Nylander, Rainer Opgen-Rhein, Klaus Schliep, Korbinian Strimmer, Damien de Vienne Maintainer: Emmanuel Paradis diff --git a/R/ace.R b/R/ace.R index 33c72de..b67549f 100644 --- a/R/ace.R +++ b/R/ace.R @@ -1,4 +1,4 @@ -## ace.R (2010-05-12) +## ace.R (2010-12-08) ## Ancestral Character Estimation @@ -174,7 +174,8 @@ as the number of categories in `x'") liks[anc, ] <- v/comp[anc] } if (output.liks) return(liks[-TIPS, ]) - -2 * sum(log(comp[-TIPS])) + dev <- -2 * sum(log(comp[-TIPS])) + if (is.na(dev)) Inf else dev } out <- nlminb(rep(ip, length.out = np), function(p) dev(p), lower = rep(0, np), upper = rep(1e50, np)) diff --git a/R/as.phylo.R b/R/as.phylo.R index c0304b5..f5a2569 100644 --- a/R/as.phylo.R +++ b/R/as.phylo.R @@ -1,4 +1,4 @@ -## as.phylo.R (2010-11-30) +## as.phylo.R (2010-12-14) ## Conversion Among Tree Objects @@ -88,6 +88,7 @@ as.hclust.phylo <- function(x, ...) { if (!is.ultrametric(x)) stop("the tree is not ultrametric") if (!is.binary.tree(x)) stop("the tree is not binary") + if (!is.rooted(x)) stop("the tree is not rooted") n <- length(x$tip.label) bt <- sort(branching.times(x)) inode <- as.numeric(names(bt)) diff --git a/R/compar.gee.R b/R/compar.gee.R index b2a9c1e..8c9bf86 100644 --- a/R/compar.gee.R +++ b/R/compar.gee.R @@ -45,7 +45,7 @@ compar.gee <- } id <- rep(1, dim(R)[1]) - geemod <<- geemod <- do.call("gee", list(formula, id, data = data, family = family, R = R, + geemod <- do.call("gee", list(formula, id, data = data, family = family, R = R, corstr = "fixed", scale.fix = scale.fix, scale.value = scale.value)) W <- geemod$naive.variance diff --git a/R/plot.phylo.R b/R/plot.phylo.R index ca3ff1c..df246bd 100644 --- a/R/plot.phylo.R +++ b/R/plot.phylo.R @@ -1,8 +1,8 @@ -## plot.phylo.R (2008-05-08) +## plot.phylo.R (2010-08-12) ## Plot Phylogenies -## Copyright 2002-2008 Emmanuel Paradis +## Copyright 2002-2010 Emmanuel Paradis ## This file is part of the R-package `ape'. ## See the file ../COPYING for licensing issues. @@ -10,17 +10,40 @@ plot.phylo <- function(x, type = "phylogram", use.edge.length = TRUE, node.pos = NULL, show.tip.label = TRUE, show.node.label = FALSE, edge.color = "black", - edge.width = 1, font = 3, cex = par("cex"), + edge.width = 1, edge.lty = 1, font = 3, cex = par("cex"), adj = NULL, srt = 0, no.margin = FALSE, root.edge = FALSE, label.offset = 0, underscore = FALSE, x.lim = NULL, y.lim = NULL, direction = "rightwards", lab4ut = "horizontal", tip.color = "black", ...) { Ntip <- length(x$tip.label) - if (Ntip == 1) stop("found only one tip in the tree!") - Nedge <- dim(x$edge)[1] + if (Ntip == 1) { + warning("found only one tip in the tree") + return(NULL) + } if (any(tabulate(x$edge[, 1]) == 1)) - stop("there are single (non-splitting) nodes in your tree; you may need to use collapse.singles().") + stop("there are single (non-splitting) nodes in your tree; you may need to use collapse.singles()") + + .nodeHeight <- function(Ntip, Nnode, edge, Nedge, yy) + .C("node_height", as.integer(Ntip), as.integer(Nnode), + as.integer(edge[, 1]), as.integer(edge[, 2]), + as.integer(Nedge), as.double(yy), + DUP = FALSE, PACKAGE = "ape")[[6]] + + .nodeDepth <- function(Ntip, Nnode, edge, Nedge) + .C("node_depth", as.integer(Ntip), as.integer(Nnode), + as.integer(edge[, 1]), as.integer(edge[, 2]), + as.integer(Nedge), double(Ntip + Nnode), + DUP = FALSE, PACKAGE = "ape")[[6]] + + .nodeDepthEdgelength <- function(Ntip, Nnode, edge, Nedge, edge.length) + .C("node_depth_edgelength", as.integer(Ntip), + as.integer(Nnode), as.integer(edge[, 1]), + as.integer(edge[, 2]), as.integer(Nedge), + as.double(edge.length), double(Ntip + Nnode), + DUP = FALSE, PACKAGE = "ape")[[7]] + + Nedge <- dim(x$edge)[1] Nnode <- x$Nnode ROOT <- Ntip + 1 type <- match.arg(type, c("phylogram", "cladogram", "fan", @@ -28,171 +51,170 @@ plot.phylo <- function(x, type = "phylogram", use.edge.length = TRUE, direction <- match.arg(direction, c("rightwards", "leftwards", "upwards", "downwards")) if (is.null(x$edge.length)) use.edge.length <- FALSE - if (type == "unrooted" || !use.edge.length) root.edge <- FALSE + + ## the order of the last two conditions is important: + if (type %in% c("unrooted", "radial") || !use.edge.length || + is.null(x$root.edge) || !x$root.edge) root.edge <- FALSE + if (type == "fan" && root.edge) { + warning("drawing root edge with type = 'fan' is not yet supported") + root.edge <- FALSE + } + phyloORclado <- type %in% c("phylogram", "cladogram") horizontal <- direction %in% c("rightwards", "leftwards") + xe <- x$edge # to save if (phyloORclado) { ## we first compute the y-coordinates of the tips. - ## Fix from Klaus Schliep (2007-06-16): - if (!is.null(attr(x, "order"))) - if (attr(x, "order") == "pruningwise") - x <- reorder(x) - ## End of fix + phyOrder <- attr(x, "order") + ## make sure the tree is in cladewise order: + if (is.null(phyOrder) || phyOrder != "cladewise") { + x <- reorder(x) # fix from Klaus Schliep (2007-06-16) + if (!identical(x$edge, xe)) { + ## modified from Li-San Wang's fix (2007-01-23): + ereorder <- match(x$edge[, 2], xe[, 2]) + if (length(edge.color) > 1) { + edge.color <- rep(edge.color, length.out = Nedge) + edge.color <- edge.color[ereorder] + } + if (length(edge.width) > 1) { + edge.width <- rep(edge.width, length.out = Nedge) + edge.width <- edge.width[ereorder] + } + if (length(edge.lty) > 1) { + edge.lty <- rep(edge.lty, length.out = Nedge) + edge.lty <- edge.lty[ereorder] + } + } + } +### By contrats to ape (< 2.4), the arguments edge.color, etc., are +### not elongated before being passed to segments(), except if needed +### to be reordered yy <- numeric(Ntip + Nnode) TIPS <- x$edge[x$edge[, 2] <= Ntip, 2] yy[TIPS] <- 1:Ntip } - edge.color <- rep(edge.color, length.out = Nedge) - edge.width <- rep(edge.width, length.out = Nedge) - ## fix from Li-San Wang (2007-01-23): - xe <- x$edge - x <- reorder(x, order = "pruningwise") - ereorder <- match(x$edge[, 2], xe[, 2]) - edge.color <- edge.color[ereorder] - edge.width <- edge.width[ereorder] - ## End of fix + ## 'z' is the tree in pruningwise order used in calls to .C + z <- reorder(x, order = "pruningwise") +### edge.color <- rep(edge.color, length.out = Nedge) +### edge.width <- rep(edge.width, length.out = Nedge) +### edge.lty <- rep(edge.lty, length.out = Nedge) +### ## fix from Li-San Wang (2007-01-23): +### xe <- x$edge +### x <- reorder(x, order = "pruningwise") +### ereorder <- match(x$edge[, 2], xe[, 2]) +### edge.color <- edge.color[ereorder] +### edge.width <- edge.width[ereorder] +### edge.lty <- edge.lty[ereorder] +### ## end of fix if (phyloORclado) { if (is.null(node.pos)) { node.pos <- 1 if (type == "cladogram" && !use.edge.length) node.pos <- 2 } if (node.pos == 1) - yy <- .C("node_height", as.integer(Ntip), as.integer(Nnode), - as.integer(x$edge[, 1]), as.integer(x$edge[, 2]), - as.integer(Nedge), as.double(yy), - DUP = FALSE, PACKAGE = "ape")[[6]] + yy <- .nodeHeight(Ntip, Nnode, z$edge, Nedge, yy) else { ## node_height_clado requires the number of descendants ## for each node, so we compute `xx' at the same time ans <- .C("node_height_clado", as.integer(Ntip), - as.integer(Nnode), as.integer(x$edge[, 1]), - as.integer(x$edge[, 2]), as.integer(Nedge), + as.integer(Nnode), as.integer(z$edge[, 1]), + as.integer(z$edge[, 2]), as.integer(Nedge), double(Ntip + Nnode), as.double(yy), DUP = FALSE, PACKAGE = "ape") xx <- ans[[6]] - 1 yy <- ans[[7]] } if (!use.edge.length) { - if(node.pos != 2) - xx <- .C("node_depth", as.integer(Ntip), as.integer(Nnode), - as.integer(x$edge[, 1]), as.integer(x$edge[, 2]), - as.integer(Nedge), double(Ntip + Nnode), - DUP = FALSE, PACKAGE = "ape")[[6]] - 1 + if (node.pos != 2) xx <- .nodeDepth(Ntip, Nnode, z$edge, Nedge) - 1 xx <- max(xx) - xx } else { - xx <- .C("node_depth_edgelength", as.integer(Ntip), - as.integer(Nnode), as.integer(x$edge[, 1]), - as.integer(x$edge[, 2]), as.integer(Nedge), - as.double(x$edge.length), double(Ntip + Nnode), - DUP = FALSE, PACKAGE = "ape")[[7]] + xx <- .nodeDepthEdgelength(Ntip, Nnode, z$edge, Nedge, z$edge.length) } - } - if (type == "fan") { + } else switch(type, "fan" = { ## if the tips are not in the same order in tip.label ## and in edge[, 2], we must reorder the angles: we ## use `xx' to store temporarily the angles - TIPS <- xe[which(xe[, 2] <= Ntip), 2] + TIPS <- x$edge[which(x$edge[, 2] <= Ntip), 2] xx <- seq(0, 2*pi*(1 - 1/Ntip), 2*pi/Ntip) theta <- double(Ntip) theta[TIPS] <- xx theta <- c(theta, numeric(Nnode)) - theta <- .C("node_height", as.integer(Ntip), as.integer(Nnode), - as.integer(x$edge[, 1]), as.integer(x$edge[, 2]), - as.integer(Nedge), theta, DUP = FALSE, - PACKAGE = "ape")[[6]] + theta <- .nodeHeight(Ntip, Nnode, z$edge, Nedge, theta) if (use.edge.length) { - r <- .C("node_depth_edgelength", as.integer(Ntip), - as.integer(Nnode), as.integer(x$edge[, 1]), - as.integer(x$edge[, 2]), as.integer(Nedge), - as.double(x$edge.length), double(Ntip + Nnode), - DUP = FALSE, PACKAGE = "ape")[[7]] + r <- .nodeDepthEdgelength(Ntip, Nnode, z$edge, Nedge, z$edge.length) } else { - r <- .C("node_depth", as.integer(Ntip), as.integer(Nnode), - as.integer(x$edge[, 1]), as.integer(x$edge[, 2]), - as.integer(Nedge), double(Ntip + Nnode), - DUP = FALSE, PACKAGE = "ape")[[6]] + r <- .nodeDepth(Ntip, Nnode, z$edge, Nedge) r <- 1/r } xx <- r*cos(theta) yy <- r*sin(theta) - - } - if (type == "unrooted") { + }, "unrooted" = { + nb.sp <- .nodeDepth(Ntip, Nnode, z$edge, Nedge) XY <- if (use.edge.length) - unrooted.xy(Ntip, Nnode, x$edge, x$edge.length) + unrooted.xy(Ntip, Nnode, z$edge, z$edge.length, nb.sp) else - unrooted.xy(Ntip, Nnode, x$edge, rep(1, Nedge)) + unrooted.xy(Ntip, Nnode, z$edge, rep(1, Nedge), nb.sp) ## rescale so that we have only positive values xx <- XY$M[, 1] - min(XY$M[, 1]) yy <- XY$M[, 2] - min(XY$M[, 2]) - } - if (type == "radial") { - X <- .C("node_depth", as.integer(Ntip), as.integer(Nnode), - as.integer(x$edge[, 1]), as.integer(x$edge[, 2]), - as.integer(Nedge), double(Ntip + Nnode), - DUP = FALSE, PACKAGE = "ape")[[6]] + }, "radial" = { + X <- .nodeDepth(Ntip, Nnode, z$edge, Nedge) X[X == 1] <- 0 ## radius: X <- 1 - X/Ntip ## angle (1st compute the angles for the tips): yy <- c((1:Ntip)*2*pi/Ntip, rep(0, Nnode)) - Y <- .C("node_height", as.integer(Ntip), as.integer(Nnode), - as.integer(x$edge[, 1]), as.integer(x$edge[, 2]), - as.integer(Nedge), as.double(yy), - DUP = FALSE, PACKAGE = "ape")[[6]] + Y <- .nodeHeight(Ntip, Nnode, z$edge, Nedge, yy) xx <- X * cos(Y) yy <- X * sin(Y) - } - if (phyloORclado && direction != "rightwards") { - if (direction == "leftwards") { - xx <- -xx - xx <- xx - min(xx) - } + }) + if (phyloORclado) { if (!horizontal) { tmp <- yy yy <- xx xx <- tmp - min(tmp) + 1 - if (direction == "downwards") { - yy <- -yy - yy <- yy - min(yy) - } } - } - if (phyloORclado && root.edge) { - if (direction == "rightwards") xx <- xx + x$root.edge - if (direction == "upwards") yy <- yy + x$root.edge + if (root.edge) { + if (direction == "rightwards") xx <- xx + x$root.edge + if (direction == "upwards") yy <- yy + x$root.edge + } } if (no.margin) par(mai = rep(0, 4)) if (is.null(x.lim)) { if (phyloORclado) { if (horizontal) { x.lim <- c(0, NA) - tmp <- - if (show.tip.label) nchar(x$tip.label) * 0.018 * max(xx) * cex - else 0 - x.lim[2] <- - if (direction == "leftwards") max(xx[ROOT] + tmp) - else max(xx[1:Ntip] + tmp) + pin1 <- par("pin")[1] # width of the device in inches + strWi <- strwidth(x$tip.label, "inches") # id. for the tip labels + ## 1.04 comes from that we are using a regular axis system + ## with 4% on both sides of the range of x: + xx.tips <- xx[1:Ntip] * 1.04 + ## 'alp' is the conversion coefficient from + ## user coordinates to inches: + alp <- try(uniroot(function(a) max(a*xx.tips + strWi) - pin1, + c(0, 1e6))$root, silent = TRUE) + ## if the above fails, give 1/3 of the device for the tip labels: + if (is.character(alp)) tmp <- max(xx.tips)*1.5 else { + tmp <- if (show.tip.label) max(xx.tips + strWi/alp) else max(xx.tips) + } + x.lim[2] <- tmp } else x.lim <- c(1, Ntip) - } - if (type == "fan") { + } else switch(type, "fan" = { if (show.tip.label) { offset <- max(nchar(x$tip.label) * 0.018 * max(yy) * cex) x.lim <- c(min(xx) - offset, max(xx) + offset) } else x.lim <- c(min(xx), max(xx)) - } - if (type == "unrooted") { + }, "unrooted" = { if (show.tip.label) { offset <- max(nchar(x$tip.label) * 0.018 * max(yy) * cex) x.lim <- c(0 - offset, max(xx) + offset) } else x.lim <- c(0, max(xx)) - } - if (type == "radial") { + }, "radial" = { if (show.tip.label) { offset <- max(nchar(x$tip.label) * 0.03 * cex) x.lim <- c(-1 - offset, 1 + offset) } else x.lim <- c(-1, 1) - } + }) } else if (length(x.lim) == 1) { x.lim <- c(0, x.lim) if (phyloORclado && !horizontal) x.lim[1] <- 1 @@ -203,36 +225,43 @@ plot.phylo <- function(x, type = "phylogram", use.edge.length = TRUE, if (show.tip.label) -1 - max(nchar(x$tip.label) * 0.03 * cex) else -1 } + ## mirror the xx: + if (phyloORclado && direction == "leftwards") xx <- x.lim[2] - xx if (is.null(y.lim)) { if (phyloORclado) { if (horizontal) y.lim <- c(1, Ntip) else { y.lim <- c(0, NA) - tmp <- - if (show.tip.label) nchar(x$tip.label) * 0.018 * max(yy) * cex - else 0 - y.lim[2] <- - if (direction == "downwards") max(yy[ROOT] + tmp) - else max(yy[1:Ntip] + tmp) + pin2 <- par("pin")[2] # height of the device in inches + strWi <- strwidth(x$tip.label, "inches") + ## 1.04 comes from that we are using a regular axis system + ## with 4% on both sides of the range of x: + yy.tips <- yy[1:Ntip] * 1.04 + ## 'alp' is the conversion coefficient from + ## user coordinates to inches: + alp <- try(uniroot(function(a) max(a*yy.tips + strWi) - pin2, + c(0, 1e6))$root, silent = TRUE) + ## if the above fails, give 1/3 of the device for the tip labels: + if (is.character(alp)) tmp <- max(yy.tips)*1.5 else { + tmp <- if (show.tip.label) max(yy.tips + strWi/alp) else max(yy.tips) + } + y.lim[2] <- tmp } - } - if (type == "fan") { + } else switch(type, "fan" = { if (show.tip.label) { offset <- max(nchar(x$tip.label) * 0.018 * max(yy) * cex) y.lim <- c(min(yy) - offset, max(yy) + offset) } else y.lim <- c(min(yy), max(yy)) - } - if (type == "unrooted") { + }, "unrooted" = { if (show.tip.label) { offset <- max(nchar(x$tip.label) * 0.018 * max(yy) * cex) y.lim <- c(0 - offset, max(yy) + offset) } else y.lim <- c(0, max(yy)) - } - if (type == "radial") { + }, "radial" = { if (show.tip.label) { offset <- max(nchar(x$tip.label) * 0.03 * cex) y.lim <- c(-1 - offset, 1 + offset) } else y.lim <- c(-1, 1) - } + }) } else if (length(y.lim) == 1) { y.lim <- c(0, y.lim) if (phyloORclado && horizontal) y.lim[1] <- 1 @@ -241,49 +270,61 @@ plot.phylo <- function(x, type = "phylogram", use.edge.length = TRUE, if (type == "radial") y.lim[1] <- if (show.tip.label) -1 - max(nchar(x$tip.label) * 0.018 * max(yy) * cex) else -1 } + ## mirror the yy: + if (phyloORclado && direction == "downwards") yy <- y.lim[2] - yy if (phyloORclado && root.edge) { if (direction == "leftwards") x.lim[2] <- x.lim[2] + x$root.edge if (direction == "downwards") y.lim[2] <- y.lim[2] + x$root.edge } - ## fix by Klaus Schliep (2008-03-28): - asp <- if (type %in% c("fan", "radial")) 1 else NA - plot(0, type = "n", xlim = x.lim, ylim = y.lim, xlab = "", - ylab = "", xaxt = "n", yaxt = "n", bty = "n", asp = asp, ...) + asp <- if (type %in% c("fan", "radial", "unrooted")) 1 else NA # fixes by Klaus Schliep (2008-03-28 and 2010-08-12) + plot(0, type = "n", xlim = x.lim, ylim = y.lim, ann = FALSE, axes = FALSE, asp = asp, ...) if (is.null(adj)) - adj <- if (phyloORclado && direction == "leftwards") 1 else 0 - if (phyloORclado) { + adj <- if (phyloORclado && direction == "leftwards") 1 else 0 + if (phyloORclado && show.tip.label) { MAXSTRING <- max(strwidth(x$tip.label, cex = cex)) + loy <- 0 if (direction == "rightwards") { lox <- label.offset + MAXSTRING * 1.05 * adj - loy <- 0 } if (direction == "leftwards") { lox <- -label.offset - MAXSTRING * 1.05 * (1 - adj) - loy <- 0 - xx <- xx + MAXSTRING + #xx <- xx + MAXSTRING } if (!horizontal) { psr <- par("usr") - MAXSTRING <- MAXSTRING * 1.09 * (psr[4] - psr[3]) / (psr[2] - psr[1]) + MAXSTRING <- MAXSTRING * 1.09 * (psr[4] - psr[3])/(psr[2] - psr[1]) loy <- label.offset + MAXSTRING * 1.05 * adj lox <- 0 srt <- 90 + srt if (direction == "downwards") { loy <- -loy - yy <- yy + MAXSTRING + ##yy <- yy + MAXSTRING srt <- 180 + srt } } } if (type == "phylogram") { phylogram.plot(x$edge, Ntip, Nnode, xx, yy, - horizontal, edge.color, edge.width) + horizontal, edge.color, edge.width, edge.lty) } else { - if (type == "fan") - circular.plot(x$edge, Ntip, Nnode, xx, yy, theta, - r, edge.color, edge.width) - else - cladogram.plot(x$edge, xx, yy, edge.color, edge.width) + if (type == "fan") { + ereorder <- match(z$edge[, 2], x$edge[, 2]) + if (length(edge.color) > 1) { + edge.color <- rep(edge.color, length.out = Nedge) + edge.color <- edge.color[ereorder] + } + if (length(edge.width) > 1) { + edge.width <- rep(edge.width, length.out = Nedge) + edge.width <- edge.width[ereorder] + } + if (length(edge.lty) > 1) { + edge.lty <- rep(edge.lty, length.out = Nedge) + edge.lty <- edge.lty[ereorder] + } + circular.plot(z$edge, Ntip, Nnode, xx, yy, theta, + r, edge.color, edge.width, edge.lty) + } else + cladogram.plot(x$edge, xx, yy, edge.color, edge.width, edge.lty) } if (root.edge) switch(direction, @@ -292,11 +333,13 @@ plot.phylo <- function(x, type = "phylogram", use.edge.length = TRUE, "upwards" = segments(xx[ROOT], 0, xx[ROOT], x$root.edge), "downwards" = segments(xx[ROOT], yy[ROOT], xx[ROOT], yy[ROOT] + x$root.edge)) if (show.tip.label) { + if (is.expression(x$tip.label)) underscore <- TRUE if (!underscore) x$tip.label <- gsub("_", " ", x$tip.label) - if (phyloORclado) { + + if (phyloORclado) text(xx[1:Ntip] + lox, yy[1:Ntip] + loy, x$tip.label, adj = adj, font = font, srt = srt, cex = cex, col = tip.color) - } + if (type == "unrooted") { if (lab4ut == "horizontal") { y.adj <- x.adj <- numeric(Ntip) @@ -315,36 +358,24 @@ plot.phylo <- function(x, type = "phylogram", use.edge.length = TRUE, adj <- as.numeric(abs(XY$axe) > pi/2) srt <- 180*XY$axe/pi srt[as.logical(adj)] <- srt[as.logical(adj)] - 180 - ## temporary check of the values of `srt': - ## set to 0 if "-0.000001 < srt < 0" - sel <- srt > -1e-6 & srt < 0 - if (any(sel)) srt[sel] <- 0 - ## - ## `srt' takes only a single value, so we cannot vectorize this: + ## `srt' takes only a single value, so can't vectorize this: for (i in 1:Ntip) text(xx[i], yy[i], cex = cex, x$tip.label[i], adj = adj[i], font = font, srt = srt[i], col = tip.color[i]) } } if (type %in% c("fan", "radial")) { - xx.scaled <- xx[1:Ntip] - if (type == "fan") { # no need if type == "radial" - maxx <- max(abs(xx.scaled)) - if (maxx > 1) xx.scaled <- xx.scaled/maxx - } - angle <- acos(xx.scaled)*180/pi - s1 <- angle > 90 & yy[1:Ntip] > 0 - s2 <- angle < 90 & yy[1:Ntip] < 0 - s3 <- angle > 90 & yy[1:Ntip] < 0 - angle[s1] <- angle[s1] + 180 - angle[s2] <- -angle[s2] - angle[s3] <- 180 - angle[s3] + xx.tips <- xx[1:Ntip] + ## using atan2 considerably facilitates things compared to acos... + angle <- atan2(yy[1:Ntip], xx.tips)*180/pi + s <- xx.tips < 0 + angle[s] <- angle[s] + 180 adj <- numeric(Ntip) - adj[xx[1:Ntip] < 0] <- 1 - ## `srt' takes only a single value, so we cannot vectorize this: + adj[xx.tips < 0] <- 1 + ## `srt' takes only a single value, so can't vectorize this: for (i in 1:Ntip) - text(xx[i], yy[i], x$tip.label[i], font = font, cex = cex, - srt = angle[i], adj = adj[i], col = tip.color[i]) + text(xx[i], yy[i], x$tip.label[i], font = font, cex = cex, + srt = angle[i], adj = adj[i], col = tip.color[i]) } } if (show.node.label) @@ -362,8 +393,8 @@ plot.phylo <- function(x, type = "phylogram", use.edge.length = TRUE, invisible(L) } -phylogram.plot <- function(edge, Ntip, Nnode, xx, yy, - horizontal, edge.color, edge.width) +phylogram.plot <- function(edge, Ntip, Nnode, xx, yy, horizontal, + edge.color, edge.width, edge.lty) { nodes <- (Ntip + 1):(Ntip + Nnode) if (!horizontal) { @@ -374,90 +405,134 @@ phylogram.plot <- function(edge, Ntip, Nnode, xx, yy, ## un trait vertical à chaque noeud... x0v <- xx[nodes] y0v <- y1v <- numeric(Nnode) + ## store the index of each node in the 1st column of edge: + NodeInEdge1 <- vector("list", Nnode) for (i in nodes) { - j <- edge[which(edge[, 1] == i), 2] - y0v[i - Ntip] <- min(yy[j]) - y1v[i - Ntip] <- max(yy[j]) + ii <- i - Ntip + j <- NodeInEdge1[[ii]] <- which(edge[, 1] == i) + tmp <- range(yy[edge[j, 2]]) + y0v[ii] <- tmp[1] + y1v[ii] <- tmp[2] } ## ... et un trait horizontal partant de chaque tip et chaque noeud ## vers la racine - sq <- if (Nnode == 1) 1:Ntip else c(1:Ntip, nodes[-1]) - y0h <- yy[sq] - x1h <- xx[sq] - ## match() is very useful here becoz each element in edge[, 2] is - ## unique (not sure this is so useful in edge[, 1]; needs to be checked) - ## `pos' gives for each element in `sq' its index in edge[, 2] - pos <- match(sq, edge[, 2]) - x0h <- xx[edge[pos, 1]] + x0h <- xx[edge[, 1]] + x1h <- xx[edge[, 2]] + y0h <- yy[edge[, 2]] - e.w <- unique(edge.width) - if (length(e.w) == 1) width.v <- rep(e.w, Nnode) - else { - width.v <- rep(1, Nnode) - for (i in 1:Nnode) { - br <- edge[which(edge[, 1] == i + Ntip), 2] - width <- unique(edge.width[br]) - if (length(width) == 1) width.v[i] <- width - } - } - e.c <- unique(edge.color) - if (length(e.c) == 1) color.v <- rep(e.c, Nnode) - else { + nc <- length(edge.color) + nw <- length(edge.width) + nl <- length(edge.lty) + + if (nc + nw + nl == 3) { + color.v <- edge.color + width.v <- edge.width + lty.v <- edge.lty + } else { + Nedge <- dim(edge)[1] + edge.color <- rep(edge.color, length.out = Nedge) + edge.width <- rep(edge.width, length.out = Nedge) + edge.lty <- rep(edge.lty, length.out = Nedge) + DF <- data.frame(edge.color, edge.width, edge.lty, stringsAsFactors = FALSE) color.v <- rep("black", Nnode) + width.v <- rep(1, Nnode) + lty.v <- rep(1, Nnode) for (i in 1:Nnode) { - br <- which(edge[, 1] == i + Ntip) - #br <- edge[which(edge[, 1] == i + Ntip), 2] - color <- unique(edge.color[br]) - if (length(color) == 1) color.v[i] <- color + br <- NodeInEdge1[[i]] + if (length(br) > 2) { + x <- unique(DF[br, 1]) + if (length(x) == 1) color.v[i] <- x + x <- unique(DF[br, 2]) + if (length(x) == 1) width.v[i] <- x + x <- unique(DF[br, 3]) + if (length(x) == 1) lty.v[i] <- x + } else { + A <- br[1] + B <- br[2] + if (any(DF[A, ] != DF[B, ])) { + color.v[i] <- edge.color[B] + width.v[i] <- edge.width[B] + lty.v[i] <- edge.lty[B] + ## add a new line: + y0v <- c(y0v, y0v[i]) + y1v <- c(y1v, yy[i + Ntip]) + x0v <- c(x0v, x0v[i]) + color.v <- c(color.v, edge.color[A]) + width.v <- c(width.v, edge.width[A]) + lty.v <- c(lty.v, edge.lty[A]) + ## shorten the line: + y0v[i] <- yy[i + Ntip] + } else { + color.v[i] <- edge.color[A] + width.v[i] <- edge.width[A] + lty.v[i] <- edge.lty[A] + } + } } } - ## we need to reorder `edge.color' and `edge.width': - edge.width <- edge.width[pos] - edge.color <- edge.color[pos] if (horizontal) { - segments(x0v, y0v, x0v, y1v, col = color.v, lwd = width.v) # draws vertical lines - segments(x0h, y0h, x1h, y0h, col = edge.color, lwd = edge.width) # draws horizontal lines + segments(x0h, y0h, x1h, y0h, col = edge.color, lwd = edge.width, lty = edge.lty) # draws horizontal lines + segments(x0v, y0v, x0v, y1v, col = color.v, lwd = width.v, lty = lty.v) # draws vertical lines } else { - segments(y0v, x0v, y1v, x0v, col = color.v, lwd = width.v) # draws horizontal lines - segments(y0h, x0h, y0h, x1h, col = edge.color, lwd = edge.width) # draws vertical lines + segments(y0h, x0h, y0h, x1h, col = edge.color, lwd = edge.width, lty = edge.lty) # draws vertical lines + segments(y0v, x0v, y1v, x0v, col = color.v, lwd = width.v, lty = lty.v) # draws horizontal lines } } -cladogram.plot <- function(edge, xx, yy, edge.color, edge.width) - segments(xx[edge[, 1]], yy[edge[, 1]], xx[edge[, 2]], yy[edge[, 2]], - col = edge.color, lwd = edge.width) +cladogram.plot <- function(edge, xx, yy, edge.color, edge.width, edge.lty) + segments(xx[edge[, 1]], yy[edge[, 1]], xx[edge[, 2]], yy[edge[, 2]], + col = edge.color, lwd = edge.width, lty = edge.lty) circular.plot <- function(edge, Ntip, Nnode, xx, yy, theta, - r, edge.color, edge.width) + r, edge.color, edge.width, edge.lty) +### 'edge' must be in pruningwise order { r0 <- r[edge[, 1]] r1 <- r[edge[, 2]] theta0 <- theta[edge[, 2]] + costheta0 <- cos(theta0) + sintheta0 <- sin(theta0) - x0 <- r0*cos(theta0) - y0 <- r0*sin(theta0) - x1 <- r1*cos(theta0) - y1 <- r1*sin(theta0) + x0 <- r0 * costheta0 + y0 <- r0 * sintheta0 + x1 <- r1 * costheta0 + y1 <- r1 * sintheta0 - segments(x0, y0, x1, y1, col = edge.color, lwd = edge.width) + segments(x0, y0, x1, y1, col = edge.color, lwd = edge.width, lty = edge.lty) tmp <- which(diff(edge[, 1]) != 0) start <- c(1, tmp + 1) - end <- c(tmp, dim(edge)[1]) + Nedge <- dim(edge)[1] + end <- c(tmp, Nedge) + + ## function dispatching the features to the arcs + foo <- function(edge.feat, default) { + if (length(edge.feat) == 1) return(rep(edge.feat, Nnode)) + else { + edge.feat <- rep(edge.feat, length.out = Nedge) + feat.arc <- rep(default, Nnode) + for (k in 1:Nnode) { + tmp <- edge.feat[start[k]] + if (tmp == edge.feat[end[k]]) feat.arc[k] <- tmp + } + } + feat.arc + } + co <- foo(edge.color, "black") + lw <- foo(edge.width, 1) + ly <- foo(edge.lty, 1) for (k in 1:Nnode) { i <- start[k] j <- end[k] X <- rep(r[edge[i, 1]], 100) Y <- seq(theta[edge[i, 2]], theta[edge[j, 2]], length.out = 100) - co <- if (edge.color[i] == edge.color[j]) edge.color[i] else "black" - lw <- if (edge.width[i] == edge.width[j]) edge.width[i] else 1 - lines(X*cos(Y), X*sin(Y), col = co, lwd = lw) + lines(X*cos(Y), X*sin(Y), col = co[k], lwd = lw[k], lty = ly[k]) } } -unrooted.xy <- function(Ntip, Nnode, edge, edge.length) +unrooted.xy <- function(Ntip, Nnode, edge, edge.length, nb.sp) { foo <- function(node, ANGLE, AXIS) { ind <- which(edge[, 1] == node) @@ -474,19 +549,13 @@ unrooted.xy <- function(Ntip, Nnode, edge, edge.length) for (i in sons) if (i > Ntip) foo(i, angle[i], axis[i]) } - root <- Ntip + 1 Nedge <- dim(edge)[1] yy <- xx <- numeric(Ntip + Nnode) - nb.sp <- .C("node_depth", as.integer(Ntip), as.integer(Nnode), - as.integer(edge[, 1]), as.integer(edge[, 2]), - as.integer(Nedge), double(Ntip + Nnode), - DUP = FALSE, PACKAGE = "ape")[[6]] ## `angle': the angle allocated to each node wrt their nb of tips ## `axis': the axis of each branch axis <- angle <- numeric(Ntip + Nnode) ## start with the root... - ## xx[root] <- yy[root] <- 0 # already set! - foo(root, 2*pi, 0) + foo(Ntip + 1L, 2*pi, 0) M <- cbind(xx, yy) axe <- axis[1:Ntip] # the axis of the terminal branches (for export) diff --git a/Thanks b/Thanks index aa4d19a..402c81a 100644 --- a/Thanks +++ b/Thanks @@ -4,11 +4,12 @@ APE at one stage or another. Many users gave important feed-back with their encouragements, comments, or bug reports: thanks to all of you! -Significant bug fixes were provided by Cécile Ané, James Bullard, Otto -Cordero, Éric Durand, Olivier François, Rich FitzJohn, , Jos Käfer, -Bret Larget, Naim Matasci, Nick Matzke, Michael Phelan, Elizabeth -Purdom, Dan Rabosky, Filipe Vieira, Tim Wallstrom, Li-San Wang, Yan -Wong, Peter Wragg, and Janet Young. Contact me if I forgot someone. +Significant bug fixes were provided by Cécile Ané, Jeremy Beaulieu, +James Bullard, Otto Cordero, Éric Durand, Olivier François, Rich +FitzJohn, , Jos Käfer, Bret Larget, Naim Matasci, Nick Matzke, Michael +Phelan, Elizabeth Purdom, Dan Rabosky, Filipe Vieira, Tim Wallstrom, +Li-San Wang, Yan Wong, Peter Wragg, and Janet Young. Contact me if I +forgot someone. Kurt Hornik, of the R Core Team, helped in several occasions to fix some problems and bugs. -- 2.39.2