From b9e04a6e6af3beda74b916eda00b42ac38875563 Mon Sep 17 00:00:00 2001 From: paradis Date: Tue, 2 Oct 2012 04:00:29 +0000 Subject: [PATCH] some updates for ape 3.0-6 git-svn-id: https://svn.mpl.ird.fr/ape/dev/ape@196 6e262413-ae40-0410-9e79-b911bd7a66b7 --- NEWS | 17 +++++++++++++-- R/CDF.birth.death.R | 11 ++++++---- R/DNA.R | 16 +++++++------- R/dbd.R | 6 +---- R/evonet.R | 6 ++--- R/ladderize.R | 4 ++-- R/me.R | 22 +++++++++++-------- R/pic.R | 6 ++--- R/plot.phylo.R | 53 ++++++++++++++++++++++++++++++++++++++++++++- R/read.nexus.R | 4 +++- R/read.tree.R | 5 +++-- R/rtree.R | 3 ++- man/ape-internal.Rd | 2 -- man/dbd.Rd | 9 ++++---- man/node.depth.Rd | 22 ++++++++++++++----- man/rTraitCont.Rd | 2 +- src/plot_phylo.c | 21 +++--------------- 17 files changed, 138 insertions(+), 71 deletions(-) diff --git a/NEWS b/NEWS index 46af03f..729e3ed 100644 --- a/NEWS +++ b/NEWS @@ -5,13 +5,18 @@ NEW FEATURES o reorder.phylo() has a new order, "postorder", and a new option index.only = TRUE to return only the vector of indices (the tree - is unmodified; see ?reorder.phylo for details). + is unmodified, see ?reorder.phylo for details). + + o The three new functions node.depth.length, node.height, and + node.height.clado make some internal code available from R. See + ?node.depth (which was already available and documented) for + details. BUG FIXES o reorder(, "pruningwise") made R crash if the rows of the edge - matrix are in random order. + matrix are in random order: this is now fixed. OTHER CHANGES @@ -22,6 +27,14 @@ OTHER CHANGES visible for small trees (n < 1000) but this can be more than 1000 faster for big trees (n >= 1e4). + o The attribute "order" of the objects of class "phylo" is now + strongly recommended, though not mandatory. Most functions in + ape should return a tree with this attribute correctly set. + + o dbd() is now vectorized on both arguments 'x' (number of species + in clade) and 't' (clade age) to make likelihood calculations + easier and faster. + CHANGES IN APE VERSION 3.0-5 diff --git a/R/CDF.birth.death.R b/R/CDF.birth.death.R index 847c37d..7ba5f5e 100644 --- a/R/CDF.birth.death.R +++ b/R/CDF.birth.death.R @@ -1,9 +1,9 @@ -## CDF.birth.death.R (2010-09-27) +## CDF.birth.death.R (2012-09-14) -## Functions to simulate and fit +## Functions to Simulate and Fit ## Time-Dependent Birth-Death Models -## Copyright 2010 Emmanuel Paradis +## Copyright 2010-2012 Emmanuel Paradis ## This file is part of the R-package `ape'. ## See the file ../COPYING for licensing issues. @@ -158,7 +158,9 @@ CDF.birth.death <- if (!case %in% c(1, 3, 6)) Pi <- Vectorize(Pi) - denom <- if (fast) integrateTrapeze(Pi, 0, Tmax) else integrate(Pi, 0, Tmax)$value + denom <- + if (fast) integrateTrapeze(Pi, 0, Tmax) + else integrate(Pi, 0, Tmax)$value n <- length(x) p <- numeric(n) if (fast) { @@ -178,6 +180,7 @@ CDF.birth.death <- phy <- list(edge = edge, edge.length = edge.length, tip.label = paste("t", 1:(i + 1), sep = ""), Nnode = i) class(phy) <- "phylo" + attr(phy, "order") <- "cladewise" phy } diff --git a/R/DNA.R b/R/DNA.R index 2bd7520..b426dbd 100644 --- a/R/DNA.R +++ b/R/DNA.R @@ -1,4 +1,4 @@ -## DNA.R (2012-06-19) +## DNA.R (2012-09-13) ## Manipulations and Comparisons of DNA Sequences @@ -295,7 +295,7 @@ base.freq <- function(x, freq = FALSE, all = FALSE) { if (is.list(x)) x <- unlist(x) n <- length(x) - BF <-.C("BaseProportion", x, n, double(17), + BF <-.C("BaseProportion", x, as.integer(n), double(17), DUP = FALSE, NAOK = TRUE, PACKAGE = "ape")[[3]] names(BF) <- c("a", "c", "g", "t", "r", "m", "w", "s", "k", "y", "v", "h", "d", "b", "n", "-", "?") @@ -351,8 +351,8 @@ seg.sites <- function(x) n <- n[1] } if (n == 1) return(integer(0)) - ans <- .C("SegSites", x, n, s, integer(s), - DUP = FALSE, NAOK = TRUE, PACKAGE = "ape") + ans <- .C("SegSites", x, as.integer(n), as.integer(s), + integer(s), DUP = FALSE, NAOK = TRUE, PACKAGE = "ape") which(as.logical(ans[[4]])) } @@ -397,10 +397,10 @@ dist.dna <- function(x, model = "K80", variance = FALSE, gamma = FALSE, var <- if (variance) double(Ndist) else 0 if (!gamma) gamma <- alpha <- 0 else alpha <- gamma <- 1 - d <- .C("dist_dna", x, n, s, imod, double(Ndist), BF, - as.integer(pairwise.deletion), as.integer(variance), - var, as.integer(gamma), alpha, DUP = FALSE, NAOK = TRUE, - PACKAGE = "ape") + d <- .C("dist_dna", x, as.integer(n), as.integer(s), imod, + double(Ndist), BF, as.integer(pairwise.deletion), + as.integer(variance), var, as.integer(gamma), + alpha, DUP = FALSE, NAOK = TRUE, PACKAGE = "ape") if (variance) var <- d[[9]] d <- d[[5]] if (imod == 11) { diff --git a/R/dbd.R b/R/dbd.R index 3a325dc..b0beb7a 100644 --- a/R/dbd.R +++ b/R/dbd.R @@ -1,4 +1,4 @@ -## dbd.R (2012-03-19) +## dbd.R (2012-09-14) ## Probability Density Under Birth--Death Models @@ -27,10 +27,6 @@ dbd <- function(x, lambda, mu, t, conditional = FALSE, log = FALSE) mu <- mu[1] warning("only the first value of 'mu' was considered") } - if (length(t) > 1) { - t <- t[1] - warning("only the first value of 't' was considered") - } if (mu == 0) return(dyule(x, lambda, t, log)) diff --git a/R/evonet.R b/R/evonet.R index 9dbaaa3..7f6b45f 100644 --- a/R/evonet.R +++ b/R/evonet.R @@ -1,8 +1,8 @@ -## evonet.R (2011-06-09) +## evonet.R (2012-09-14) ## Evolutionary Networks -## Copyright 2011 Emmanuel Paradis +## Copyright 2011-2012 Emmanuel Paradis ## This file is part of the R-package `ape'. ## See the file ../COPYING for licensing issues. @@ -67,7 +67,7 @@ as.networx.evonet <- function(x, weight = NA, ...) { if (any(x$reticulation <= Ntip(x))) stop("some tips are involved in reticulations: cannot convert to \"networx\"") - x <- reorder(x, "p") + x <- reorder(x, "pruningwise") ned <- Nedge(x) nrt <- nrow(x$reticulation) x$edge <- rbind(x$edge, x$reticulation) diff --git a/R/ladderize.R b/R/ladderize.R index 8d2e94a..596a8fd 100644 --- a/R/ladderize.R +++ b/R/ladderize.R @@ -23,7 +23,7 @@ ladderize <- function(phy, right = TRUE) start <- start[o] neworder[newpos] <<- start for (i in 1:Nclade) - if (desc[i] > nb.tip) foo(desc[i], end[i], newpos[i] + 1) + if (desc[i] > nb.tip) foo(desc[i], end[i], newpos[i] + 1) } nb.tip <- length(phy$tip.label) nb.node <- phy$Nnode @@ -37,6 +37,6 @@ ladderize <- function(phy, right = TRUE) foo(nb.tip + 1, nb.edge, 1) phy$edge <- phy$edge[neworder, ] if (!is.null(phy$edge.length)) - phy$edge.length <- phy$edge.length[neworder] + phy$edge.length <- phy$edge.length[neworder] phy } diff --git a/R/me.R b/R/me.R index 1006c64..533c40e 100644 --- a/R/me.R +++ b/R/me.R @@ -1,4 +1,4 @@ -## me.R (2012-04-30) +## me.R (2012-09-14) ## Tree Estimation Based on Minimum Evolution Algorithm @@ -20,10 +20,12 @@ fastme.bal <- function(X, nni = TRUE, spr = TRUE, tbr = TRUE) labels <- attr(X, "Labels") if (is.null(labels)) labels <- as.character(1:N) labels <- labels[ans[[3]]] - structure(list(edge = cbind(ans[[7]], ans[[8]]), - edge.length = ans[[9]], - tip.label = labels, Nnode = N - 2L), - class = "phylo") + obj <- list(edge = cbind(ans[[7]], ans[[8]]), + edge.length = ans[[9]], + tip.label = labels, Nnode = N - 2L) + class(obj) <- "phylo" + attr(obj, "order") <- "cladewise" + obj } fastme.ols <- function(X, nni = TRUE) @@ -37,10 +39,12 @@ fastme.ols <- function(X, nni = TRUE) labels <- attr(X, "Labels") if (is.null(labels)) labels <- as.character(1:N) labels <- labels[ans[[3]]] - structure(list(edge = cbind(ans[[5]], ans[[6]]), - edge.length = ans[[7]], - tip.label = labels, Nnode = N - 2L), - class = "phylo") + obj <- list(edge = cbind(ans[[5]], ans[[6]]), + edge.length = ans[[7]], + tip.label = labels, Nnode = N - 2L) + class(obj) <- "phylo" + attr(obj, "order") <- "cladewise" + obj } bionj <- function(X) diff --git a/R/pic.R b/R/pic.R index 816a09b..b1cb20b 100644 --- a/R/pic.R +++ b/R/pic.R @@ -1,8 +1,8 @@ -## pic.R (2011-03-01) +## pic.R (2012-09-11) ## Phylogenetically Independent Contrasts -## Copyright 2002-2011 Emmanuel Paradis +## Copyright 2002-2012 Emmanuel Paradis ## This file is part of the R-package `ape'. ## See the file ../COPYING for licensing issues. @@ -22,7 +22,7 @@ pic <- function(x, phy, scaled = TRUE, var.contrasts = FALSE) if (any(is.na(x))) stop("missing data in 'x': you may consider removing the species with missing data from your tree with the function 'drop.tip'.") - phy <- reorder(phy, "pruningwise") + phy <- reorder(phy, "postorder") phenotype <- numeric(nb.tip + nb.node) if (is.null(names(x))) { diff --git a/R/plot.phylo.R b/R/plot.phylo.R index 3fd0637..862291d 100644 --- a/R/plot.phylo.R +++ b/R/plot.phylo.R @@ -1,4 +1,4 @@ -## plot.phylo.R (2012-03-22) +## plot.phylo.R (2012-10-02) ## Plot Phylogenies @@ -595,6 +595,57 @@ node.depth <- function(phy) as.integer(N), double(n + m), DUP = FALSE, PACKAGE = "ape")[[6]] } +node.depth.edgelength <- function(phy) +{ + n <- length(phy$tip.label) + m <- phy$Nnode + N <- dim(phy$edge)[1] + phy <- reorder(phy, order = "pruningwise") + .C("node_depth_edgelength", as.integer(n), as.integer(n), + as.integer(phy$edge[, 1]), as.integer(phy$edge[, 2]), + as.integer(N), as.double(phy$edge.length), double(n + m), + DUP = FALSE, PACKAGE = "ape")[[7]] +} + +node.height <- function(phy) +{ + n <- length(phy$tip.label) + m <- phy$Nnode + N <- dim(phy$edge)[1] + phy <- reorder(phy, order = "pruningwise") + + e1 <- phy$edge[, 1] + e2 <- phy$edge[, 2] + + yy <- numeric(n + m) + TIPS <- e2[e2 <= n] + yy[TIPS] <- 1:n + + .C("node_height", as.integer(n), as.integer(m), + as.integer(e1), as.integer(e2), as.integer(N), + as.double(yy), DUP = FALSE, PACKAGE = "ape")[[6]] +} + +node.height.clado <- function(phy) +{ + n <- length(phy$tip.label) + m <- phy$Nnode + N <- dim(phy$edge)[1] + phy <- reorder(phy, order = "pruningwise") + + e1 <- phy$edge[, 1] + e2 <- phy$edge[, 2] + + yy <- numeric(n + m) + TIPS <- e2[e2 <= n] + yy[TIPS] <- 1:n + + .C("node_height_clado", as.integer(n), as.integer(m), + as.integer(e1), as.integer(e2), as.integer(N), + double(n + m), as.double(yy), DUP = FALSE, + PACKAGE = "ape")[[7]] +} + plot.multiPhylo <- function(x, layout = 1, ...) { if (layout > 1) diff --git a/R/read.nexus.R b/R/read.nexus.R index 0268e68..2f4177e 100644 --- a/R/read.nexus.R +++ b/R/read.nexus.R @@ -1,4 +1,4 @@ -## read.nexus.R (2012-02-09) +## read.nexus.R (2012-09-28) ## Read Tree File in Nexus Format @@ -16,6 +16,7 @@ names(phy) <- nms if (all(phy$node.label == "")) phy$node.label <- NULL class(phy) <- "phylo" + attr(phy, "order") <- "cladewise" phy } @@ -97,6 +98,7 @@ clado.build <- function(tp) if (all(obj$node.label == "NA")) NULL else gsub("^NA", "", obj$node.label) class(obj) <- "phylo" + attr(obj, "order") <- "cladewise" obj } diff --git a/R/read.tree.R b/R/read.tree.R index cba51d1..47a432c 100644 --- a/R/read.tree.R +++ b/R/read.tree.R @@ -1,8 +1,8 @@ -## read.tree.R (2010-09-27) +## read.tree.R (2012-09-14) ## Read Tree Files in Parenthetic Format -## Copyright 2002-2010 Emmanuel Paradis, Daniel Lawson and Klaus Schliep +## Copyright 2002-2012 Emmanuel Paradis, Daniel Lawson and Klaus Schliep ## This file is part of the R-package `ape'. ## See the file ../COPYING for licensing issues. @@ -95,6 +95,7 @@ tree.build <- function(tp) if (any(nzchar(node.label))) obj$node.label <- node.label if (!is.na(root.edge)) obj$root.edge <- root.edge class(obj) <- "phylo" + attr(obj, "order") <- "cladewise" obj } diff --git a/R/rtree.R b/R/rtree.R index 0df3596..eb00b60 100644 --- a/R/rtree.R +++ b/R/rtree.R @@ -1,4 +1,4 @@ -## rtree.R (2012-02-09) +## rtree.R (2012-09-14) ## Generates Trees @@ -100,6 +100,7 @@ rtree <- function(n, rooted = TRUE, tip.label = NULL, br = runif, ...) } phy$Nnode <- n - 2L + rooted class(phy) <- "phylo" + attr(phy, "order") <- "cladewise" phy } diff --git a/man/ape-internal.Rd b/man/ape-internal.Rd index 838a3a3..c923093 100644 --- a/man/ape-internal.Rd +++ b/man/ape-internal.Rd @@ -13,8 +13,6 @@ \alias{circular.plot} \alias{unrooted.plot} \alias{unrooted.xy} -\alias{node.height} -\alias{node.height.clado} \alias{birth.step} \alias{death.step} \alias{ht.move} diff --git a/man/dbd.Rd b/man/dbd.Rd index ad2e620..4f5ba15 100644 --- a/man/dbd.Rd +++ b/man/dbd.Rd @@ -43,11 +43,12 @@ dbdTime(x, birth, death, t, conditional = FALSE, \code{dbdTime} is for time-varying \code{lambda} and \code{mu} specified as \R functions. - Only \code{dyule} is vectorized simultaneously on its three arguments + \code{dyule} is vectorized simultaneously on its three arguments \code{x}, \code{lambda}, and \code{t}, according to \R's rules of - recycling arguments. The two others are vectorized only on \code{x}; - the other arguments are eventually shortened with a warning if - necessary. + recycling arguments. \code{dbd} is vectorized simultaneously \code{x} + and \code{t} (to make likelihood calculations easy), and + \code{dbdTime} is vectorized only on \code{x}; the other arguments are + eventually shortened with a warning if necessary. The returned value is, logically, zero for values of \code{x} out of range, i.e., negative or zero for \code{dyule} or if \code{conditional diff --git a/man/node.depth.Rd b/man/node.depth.Rd index aced56a..e4c7bf9 100644 --- a/man/node.depth.Rd +++ b/man/node.depth.Rd @@ -1,19 +1,31 @@ \name{node.depth} \alias{node.depth} -\title{Depth of Nodes and Tips} +\alias{node.depth.length} +\alias{node.height} +\alias{node.height.clado} +\title{Depth and Heights of Nodes and Tips} \description{ - This function returns the depth of nodes and tips given by the number - of descendants (1 is returned for tips). + These functions return the depth or height of nodes and tips. } \usage{ node.depth(phy) +node.depth.length(phy) +node.height(phy) +node.height.clado(phy) } \arguments{ \item{phy}{an object of class "phylo".} } \details{ - The depth of a node is computed as the number of tips which are its - descendants. The value of 1 is given to the tips. + \code{node.depth} computes the depth of a node as the number of tips + which are its descendants. The value of 1 is given to the tips. + + \code{node.depth.length} does the same but using branch lengths. + + \code{node.height} computes the heights of nodes and tips as plotted + by a phylogram. + + \code{node.height.clado} does the same but for a cladogram. } \value{ A numeric vector indexed with the node numbers of the matrix `edge' of diff --git a/man/rTraitCont.Rd b/man/rTraitCont.Rd index 9e96f9e..2251c85 100644 --- a/man/rTraitCont.Rd +++ b/man/rTraitCont.Rd @@ -73,7 +73,7 @@ rTraitCont(phy, model = "BM", sigma = 0.1, alpha = 1, theta = 0, data(bird.orders) rTraitCont(bird.orders) # BM with sigma = 0.1 ### OU model with two optima: -tr <- reorder(bird.orders, "p") +tr <- reorder(bird.orders, "postorder") plot(tr) edgelabels() theta <- rep(0, Nedge(tr)) diff --git a/src/plot_phylo.c b/src/plot_phylo.c index d8666a2..9b38d2a 100644 --- a/src/plot_phylo.c +++ b/src/plot_phylo.c @@ -1,6 +1,6 @@ -/* plot_phylo.c (2011-06-23) */ +/* plot_phylo.c (2012-10-01) */ -/* Copyright 2004-2011 Emmanuel Paradis +/* Copyright 2004-2012 Emmanuel Paradis /* This file is part of the R-package `ape'. */ /* See the file ../COPYING for licensing issues. */ @@ -35,7 +35,7 @@ void node_depth(int *ntip, int *nnode, int *edge1, int *edge2, } void node_height(int *ntip, int *nnode, int *edge1, int *edge2, - int *nedge, double *yy) + int *nedge, double *yy) { int i, n; double S; @@ -77,18 +77,3 @@ void node_height_clado(int *ntip, int *nnode, int *edge1, int *edge2, } } } - -void get_single_index_integer(int *x, int *val, int *index) -{ - while (x[*index] != *val) (*index)++; - *index += 1; -} - -void get_two_index_integer(int *x, int *val, int *index) -{ - while (x[index[0]] != *val) index[0]++; - index[1] = index[0] + 1; - while (x[index[1]] != *val) index[1]++; - index[0] += 1; - index[1] += 1; -} -- 2.39.2