From: paradis Date: Wed, 24 Nov 2010 09:26:27 +0000 (+0000) Subject: various bug fixes X-Git-Url: https://git.donarmstrong.com/?p=ape.git;a=commitdiff_plain;h=52008daf7f708f3fcdc735f22af308dd1a461670 various bug fixes git-svn-id: https://svn.mpl.ird.fr/ape/dev/ape@139 6e262413-ae40-0410-9e79-b911bd7a66b7 --- diff --git a/ChangeLog b/ChangeLog index 2858cd4..a021f9c 100644 --- a/ChangeLog +++ b/ChangeLog @@ -21,6 +21,13 @@ BUG FIXES o mcmc.popsize() terminated unexpectedly if the progress bar was turned off. + o prop.part(x) made R frozen if 'x' is of class "multiPhylo". + + o Compilation under Mandriva failed (thanks to Jos Käfer for the fix). + + o drop.tip() shuffled tip labels with subtree = TRUE or trim.internal + = FALSE. + CHANGES IN APE VERSION 2.6-1 diff --git a/DESCRIPTION b/DESCRIPTION index b9d0974..04f43c6 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: ape Version: 2.6-2 -Date: 2010-11-18 +Date: 2010-11-24 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/dist.topo.R b/R/dist.topo.R index 86b396b..fa9e869 100644 --- a/R/dist.topo.R +++ b/R/dist.topo.R @@ -1,4 +1,4 @@ -## dist.topo.R (2010-05-25) +## dist.topo.R (2010-11-18) ## Topological Distances, Tree Bipartitions, ## Consensus Trees, and Bootstrapping Phylogenies @@ -100,17 +100,16 @@ prop.part <- function(..., check.labels = TRUE) if (length(obj) == 1 && class(obj[[1]]) != "phylo") obj <- obj[[1]] ## - ## class(obj) <- NULL # needed? + ## class(obj) <- NULL # needed? apparently not, see below (2010-11-18) ## ntree <- length(obj) if (ntree == 1) check.labels <- FALSE - if (check.labels) obj <- .compressTipLabel(obj) + if (check.labels) .compressTipLabel(obj) # no need to store cause uncompress later for (i in 1:ntree) storage.mode(obj[[i]]$Nnode) <- "integer" ## ## The 1st must have tip labels ## Maybe simply pass the number of tips to the C code?? - if (!is.null(attr(obj, "TipLabel"))) - for (i in 1:ntree) obj[[i]]$tip.label <- attr(obj, "TipLabel") + obj <- .uncompressTipLabel(obj) # fix a bug (2010-11-18) ## clades <- .Call("prop_part", obj, ntree, TRUE, PACKAGE = "ape") attr(clades, "number") <- attr(clades, "number")[1:length(clades)] diff --git a/R/drop.tip.R b/R/drop.tip.R index c5eb07c..948b226 100644 --- a/R/drop.tip.R +++ b/R/drop.tip.R @@ -1,4 +1,4 @@ -## drop.tip.R (2010-07-21) +## drop.tip.R (2010-11-24) ## Remove Tips in a Phylogenetic Tree @@ -96,6 +96,8 @@ drop.tip <- if (is.character(tip)) tip <- which(phy$tip.label %in% tip) } + if (any(tip > Ntip)) + warning("some tip numbers were higher than the number of tips") if (!rooted && subtree) { phy <- root(phy, (1:Ntip)[-tip][1]) @@ -119,15 +121,6 @@ drop.tip <- edge2 <- phy$edge[, 2] # keep <- !logical(Nedge) - ## find the tips to drop: - if (is.character(tip)) - tip <- which(phy$tip.label %in% tip) - - if (!rooted && subtree) { - phy <- root(phy, (1:Ntip)[-tip][1]) - root.edge <- 0 - } - ## delete the terminal edges given by `tip': keep[match(tip, edge2)] <- FALSE @@ -179,21 +172,27 @@ drop.tip <- ## get the old No. of the nodes and tips that become tips: oldNo.ofNewTips <- phy$edge[TERMS, 2] + ## in case some tips are dropped but kept because of 'subtree = TRUE': + if (subtree) { + i <- which(tip %in% oldNo.ofNewTips) + if (length(i)) { + phy$tip.label[tip[i]] <- "[1_tip]" + tip <- tip[-i] + } + } + n <- length(oldNo.ofNewTips) # the new number of tips in the tree - ## the tips may not be sorted in increasing order of their - ## in the 2nd col of edge, so no need to reorder $tip.label + ## the tips may not be sorted in increasing order in the + ## 2nd col of edge, so no need to reorder $tip.label phy$edge[TERMS, 2] <- rank(phy$edge[TERMS, 2]) + phy$tip.label <- phy$tip.label[-tip] ## make new tip labels if necessary: if (subtree || !trim.internal) { - ## get the logical indices of the tips that are kept within 'oldNo.ofNewTips': - tips.kept <- oldNo.ofNewTips <= Ntip & !(oldNo.ofNewTips %in% tip) - new.tip.label <- character(n) - new.tip.label[tips.kept] <- phy$tip.label[-tip] ## get the numbers of the nodes that become tips: - node2tip <- oldNo.ofNewTips[!tips.kept] - new.tip.label[!tips.kept] <- if (subtree) { + node2tip <- oldNo.ofNewTips[oldNo.ofNewTips > Ntip] + new.tip.label <- if (subtree) { paste("[", N[node2tip], "_tips]", sep = "") } else { if (is.null(phy$node.label)) rep("NA", length(node2tip)) @@ -201,8 +200,8 @@ drop.tip <- } if (!is.null(phy$node.label)) phy$node.label <- phy$node.label[-(node2tip - Ntip)] - phy$tip.label <- new.tip.label - } else phy$tip.label <- phy$tip.label[-tip] + phy$tip.label <- c(phy$tip.label, new.tip.label) + } ## update node.label if needed: if (!is.null(phy$node.label)) @@ -222,3 +221,4 @@ drop.tip <- storage.mode(phy$edge) <- "integer" collapse.singles(phy) } + diff --git a/R/plot.phylo.R b/R/plot.phylo.R index df246bd..ca3ff1c 100644 --- a/R/plot.phylo.R +++ b/R/plot.phylo.R @@ -1,8 +1,8 @@ -## plot.phylo.R (2010-08-12) +## plot.phylo.R (2008-05-08) ## Plot Phylogenies -## Copyright 2002-2010 Emmanuel Paradis +## Copyright 2002-2008 Emmanuel Paradis ## This file is part of the R-package `ape'. ## See the file ../COPYING for licensing issues. @@ -10,40 +10,17 @@ 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, edge.lty = 1, font = 3, cex = par("cex"), + edge.width = 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) { - 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()") - - .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]] - + if (Ntip == 1) stop("found only one tip in the tree!") Nedge <- dim(x$edge)[1] + if (any(tabulate(x$edge[, 1]) == 1)) + stop("there are single (non-splitting) nodes in your tree; you may need to use collapse.singles().") Nnode <- x$Nnode ROOT <- Ntip + 1 type <- match.arg(type, c("phylogram", "cladogram", "fan", @@ -51,170 +28,171 @@ 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 - - ## 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 - } - + if (type == "unrooted" || !use.edge.length) 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. - 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 + ## Fix from Klaus Schliep (2007-06-16): + if (!is.null(attr(x, "order"))) + if (attr(x, "order") == "pruningwise") + x <- reorder(x) + ## End of fix yy <- numeric(Ntip + Nnode) TIPS <- x$edge[x$edge[, 2] <= Ntip, 2] yy[TIPS] <- 1:Ntip } - ## '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 + 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 if (phyloORclado) { if (is.null(node.pos)) { node.pos <- 1 if (type == "cladogram" && !use.edge.length) node.pos <- 2 } if (node.pos == 1) - yy <- .nodeHeight(Ntip, Nnode, z$edge, Nedge, yy) + 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]] 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(z$edge[, 1]), - as.integer(z$edge[, 2]), as.integer(Nedge), + as.integer(Nnode), as.integer(x$edge[, 1]), + as.integer(x$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 <- .nodeDepth(Ntip, Nnode, z$edge, Nedge) - 1 + 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 xx <- max(xx) - xx } else { - xx <- .nodeDepthEdgelength(Ntip, Nnode, z$edge, Nedge, z$edge.length) + 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]] } - } else switch(type, "fan" = { + } + if (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 <- x$edge[which(x$edge[, 2] <= Ntip), 2] + TIPS <- xe[which(xe[, 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 <- .nodeHeight(Ntip, Nnode, z$edge, Nedge, theta) + 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]] if (use.edge.length) { - r <- .nodeDepthEdgelength(Ntip, Nnode, z$edge, Nedge, z$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]] } else { - r <- .nodeDepth(Ntip, Nnode, z$edge, Nedge) + 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 <- 1/r } xx <- r*cos(theta) yy <- r*sin(theta) - }, "unrooted" = { - nb.sp <- .nodeDepth(Ntip, Nnode, z$edge, Nedge) + + } + if (type == "unrooted") { XY <- if (use.edge.length) - unrooted.xy(Ntip, Nnode, z$edge, z$edge.length, nb.sp) + unrooted.xy(Ntip, Nnode, x$edge, x$edge.length) else - unrooted.xy(Ntip, Nnode, z$edge, rep(1, Nedge), nb.sp) + unrooted.xy(Ntip, Nnode, x$edge, rep(1, Nedge)) ## rescale so that we have only positive values xx <- XY$M[, 1] - min(XY$M[, 1]) yy <- XY$M[, 2] - min(XY$M[, 2]) - }, "radial" = { - X <- .nodeDepth(Ntip, Nnode, z$edge, Nedge) + } + 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]] 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 <- .nodeHeight(Ntip, Nnode, z$edge, Nedge, yy) + 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]] xx <- X * cos(Y) yy <- X * sin(Y) - }) - if (phyloORclado) { + } + if (phyloORclado && direction != "rightwards") { + if (direction == "leftwards") { + xx <- -xx + xx <- xx - min(xx) + } if (!horizontal) { tmp <- yy yy <- xx xx <- tmp - min(tmp) + 1 - } - if (root.edge) { - if (direction == "rightwards") xx <- xx + x$root.edge - if (direction == "upwards") yy <- yy + x$root.edge + 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 (no.margin) par(mai = rep(0, 4)) if (is.null(x.lim)) { if (phyloORclado) { if (horizontal) { x.lim <- c(0, NA) - 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 + 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) } else x.lim <- c(1, Ntip) - } else switch(type, "fan" = { + } + if (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)) - }, "unrooted" = { + } + if (type == "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)) - }, "radial" = { + } + if (type == "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 @@ -225,43 +203,36 @@ 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) - 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 + 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) } - } else switch(type, "fan" = { + } + if (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)) - }, "unrooted" = { + } + if (type == "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)) - }, "radial" = { + } + if (type == "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 @@ -270,61 +241,49 @@ 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 } - 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, ...) + ## 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, ...) if (is.null(adj)) - adj <- if (phyloORclado && direction == "leftwards") 1 else 0 - if (phyloORclado && show.tip.label) { + adj <- if (phyloORclado && direction == "leftwards") 1 else 0 + if (phyloORclado) { 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) - #xx <- xx + MAXSTRING + loy <- 0 + 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, edge.lty) + horizontal, edge.color, edge.width) } else { - 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 (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 (root.edge) switch(direction, @@ -333,13 +292,11 @@ 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) @@ -358,24 +315,36 @@ 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 - ## `srt' takes only a single value, so can't vectorize this: + ## 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: 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.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 + 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] adj <- numeric(Ntip) - adj[xx.tips < 0] <- 1 - ## `srt' takes only a single value, so can't vectorize this: + adj[xx[1:Ntip] < 0] <- 1 + ## `srt' takes only a single value, so we cannot 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) @@ -393,8 +362,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, edge.lty) +phylogram.plot <- function(edge, Ntip, Nnode, xx, yy, + horizontal, edge.color, edge.width) { nodes <- (Ntip + 1):(Ntip + Nnode) if (!horizontal) { @@ -405,134 +374,90 @@ phylogram.plot <- function(edge, Ntip, Nnode, xx, yy, horizontal, ## 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) { - ii <- i - Ntip - j <- NodeInEdge1[[ii]] <- which(edge[, 1] == i) - tmp <- range(yy[edge[j, 2]]) - y0v[ii] <- tmp[1] - y1v[ii] <- tmp[2] + j <- edge[which(edge[, 1] == i), 2] + y0v[i - Ntip] <- min(yy[j]) + y1v[i - Ntip] <- max(yy[j]) } ## ... et un trait horizontal partant de chaque tip et chaque noeud ## vers la racine - x0h <- xx[edge[, 1]] - x1h <- xx[edge[, 2]] - y0h <- yy[edge[, 2]] - - nc <- length(edge.color) - nw <- length(edge.width) - nl <- length(edge.lty) + 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]] - 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) + e.w <- unique(edge.width) + if (length(e.w) == 1) width.v <- rep(e.w, Nnode) + else { width.v <- rep(1, Nnode) - lty.v <- rep(1, Nnode) for (i in 1:Nnode) { - 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] - } - } + 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 { + color.v <- rep("black", 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 } } + ## we need to reorder `edge.color' and `edge.width': + edge.width <- edge.width[pos] + edge.color <- edge.color[pos] if (horizontal) { - 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 + 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 } else { - 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 + 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 } } -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) +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) circular.plot <- function(edge, Ntip, Nnode, xx, yy, theta, - r, edge.color, edge.width, edge.lty) -### 'edge' must be in pruningwise order + r, edge.color, edge.width) { r0 <- r[edge[, 1]] r1 <- r[edge[, 2]] theta0 <- theta[edge[, 2]] - costheta0 <- cos(theta0) - sintheta0 <- sin(theta0) - x0 <- r0 * costheta0 - y0 <- r0 * sintheta0 - x1 <- r1 * costheta0 - y1 <- r1 * sintheta0 + x0 <- r0*cos(theta0) + y0 <- r0*sin(theta0) + x1 <- r1*cos(theta0) + y1 <- r1*sin(theta0) - segments(x0, y0, x1, y1, col = edge.color, lwd = edge.width, lty = edge.lty) + segments(x0, y0, x1, y1, col = edge.color, lwd = edge.width) tmp <- which(diff(edge[, 1]) != 0) start <- c(1, tmp + 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) + end <- c(tmp, dim(edge)[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) - lines(X*cos(Y), X*sin(Y), col = co[k], lwd = lw[k], lty = ly[k]) + 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) } } -unrooted.xy <- function(Ntip, Nnode, edge, edge.length, nb.sp) +unrooted.xy <- function(Ntip, Nnode, edge, edge.length) { foo <- function(node, ANGLE, AXIS) { ind <- which(edge[, 1] == node) @@ -549,13 +474,19 @@ unrooted.xy <- function(Ntip, Nnode, edge, edge.length, nb.sp) 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... - foo(Ntip + 1L, 2*pi, 0) + ## xx[root] <- yy[root] <- 0 # already set! + foo(root, 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 606d7a8..eb5ecd1 100644 --- a/Thanks +++ b/Thanks @@ -5,10 +5,10 @@ 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, Bret -Larget, 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. +Otto Cordero, Éric Durand, Olivier François, Rich FitzJohn, , Jos +Käfer, Bret Larget, 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. diff --git a/man/vcv.phylo.Rd b/man/vcv.phylo.Rd index 98def31..702ca63 100644 --- a/man/vcv.phylo.Rd +++ b/man/vcv.phylo.Rd @@ -1,4 +1,4 @@ -\name{vcv.phylo} +\name{vcv} \alias{vcv} \alias{vcv.phylo} \alias{vcv.corPhyl} diff --git a/src/newick.c b/src/newick.c index d0a4f15..e94ded5 100644 --- a/src/newick.c +++ b/src/newick.c @@ -1,4 +1,4 @@ -/* newick.c 2008-01-14 */ +/* newick.c 2010-11-23 */ /* Copyright 2007-2008 Vincent Lefort */ @@ -230,7 +230,7 @@ tree *readNewickString (char *str, int numLeaves) } } centerNode = decodeNewickSubtree (str, T, &uCount); - snprintf (centerNode->label, MAX_LABEL_LENGTH, rootLabel); + snprintf (centerNode->label, MAX_LABEL_LENGTH, "%s", rootLabel); /* added "%s" following Jos Kafer's suggestion (2010-11-23) */ T->root = centerNode; return(T); }