From 48681c35f7904c17b064de9c6bf96d84581e9e52 Mon Sep 17 00:00:00 2001 From: paradis Date: Fri, 10 Feb 2012 06:53:03 +0000 Subject: [PATCH] final commit for ape 3.0 git-svn-id: https://svn.mpl.ird.fr/ape/dev/ape@177 6e262413-ae40-0410-9e79-b911bd7a66b7 --- DESCRIPTION | 2 +- NEWS | 24 ++++- R/PGLS.R | 251 ++++++++++++++++++++++----------------------- R/ace.R | 2 +- R/chronopl.R | 104 +++++++++++-------- R/is.ultrametric.R | 32 +++--- R/nodelabels.R | 11 +- R/rTrait.R | 6 +- R/read.nexus.R | 8 +- R/rtree.R | 12 +-- R/vcv.phylo.R | 68 ++++++------ R/write.nexus.R | 54 +++------- man/chronopl.Rd | 21 ++-- man/nodelabels.Rd | 4 +- man/pcoa.Rd | 2 + man/rTraitDisc.Rd | 9 ++ man/read.nexus.Rd | 6 -- man/write.nexus.Rd | 10 +- src/BIONJ.c | 14 +-- src/me.c | 14 +-- src/me_ols.c | 4 +- src/newick.c | 14 +-- 22 files changed, 341 insertions(+), 331 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index b16b849..9f97add 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: ape Version: 3.0 -Date: 2012-02-03 +Date: 2012-02-10 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, Andrei-Alin Popescu, Klaus Schliep, Korbinian Strimmer, Damien de Vienne Maintainer: Emmanuel Paradis diff --git a/NEWS b/NEWS index 0ca6a14..0a4293e 100644 --- a/NEWS +++ b/NEWS @@ -15,10 +15,18 @@ NEW FEATURES o boot.phylo() and prop.clades() have a new option rooted = FALSE. Note that the behaviour of prop.part() is unchanged. + o edgelabels() has a new option 'date' to place labels on edges of + chronograms using the time scale (suggestion by Rob Lanfear). + BUG FIXES - o mantel.test() could return P-values > 1 with the default + o In chronopl(), the code setting the initial dates failed in + complicated settings (several dates known within intervals). + This has been generally improved and should result in faster + and more efficient convergence even in simple settings. + + o mantel.test() sometimes returned P-values > 1 with the default two-tailed test. o extract.clade() sometimes shuffled some tip labels (thanks to @@ -38,6 +46,16 @@ OTHER CHANGES o The C code of base.freq() has been improved and is now nearly 8 times faster. + o The option 'original.data' of write.nexus() is now deprecated and + will be removed in a future release. + + o The code of is.ultrametric() has been improved and is now 3 times + faster. + + o The code of vcv.phylo() has been improved and is now 10 or 30 + times faster for 100 or 1000 tips, respectively. Consequently, + fitting models with PGLS will be faster overall. + CHANGES IN APE VERSION 2.8 @@ -191,8 +209,8 @@ BUG FIXES OTHER CHANGES - o identify.phylo() now returns NULL if the user right-(instead of - left-)clicks (an error was returned previously). + o identify.phylo() now returns NULL if the user right- (instead of + left-) clicks (an error was returned previously). o read.nexus() should be robust to commands inserted in the TREES block. diff --git a/R/PGLS.R b/R/PGLS.R index 6b6f4af..5378ed1 100644 --- a/R/PGLS.R +++ b/R/PGLS.R @@ -1,184 +1,173 @@ -## PGLS.R (2008-05-08) +## PGLS.R (2012-02-09) ## Phylogenetic Generalized Least Squares -## Copyright 2004 Julien Dutheil, and 2006-2008 Emmanuel Paradis +## Copyright 2004 Julien Dutheil, and 2006-2012 Emmanuel Paradis ## This file is part of the R-package `ape'. ## See the file ../COPYING for licensing issues. -corBrownian <- function(value = 1, phy, form=~1) +corBrownian <- function(value = 1, phy, form = ~1) { - if (!("phylo" %in% class(phy))) - stop("object \"phy\" is not of class \"phylo\"") - attr(value, "formula") <- form - attr(value, "fixed") <- TRUE - attr(value, "tree") <- phy - class(value) <- c("corBrownian", "corPhyl", "corStruct") - return(value) + if (!inherits(phy, "phylo")) + stop('object "phy" is not of class "phylo"') + attr(value, "formula") <- form + attr(value, "fixed") <- TRUE + attr(value, "tree") <- phy + class(value) <- c("corBrownian", "corPhyl", "corStruct") + value } -corMartins <- function(value, phy, form=~1, fixed=FALSE) +corMartins <- function(value, phy, form = ~1, fixed = FALSE) { - if(length(value) > 1) - stop("only one parameter is allowed in corPGLS structure.") - if(value < 0) stop("parameter alpha must be positive.") - if (!("phylo" %in% class(phy))) - stop("object \"phy\" is not of class \"phylo\"") - attr(value, "formula") <- form - attr(value, "fixed") <- fixed - attr(value, "tree") <- phy - class(value) <- c("corMartins", "corPhyl", "corStruct") - return(value) + if (length(value) > 1) + stop("only one parameter is allowed") + if (value < 0) stop("the parameter alpha must be positive") + if (!inherits(phy, "phylo")) + stop('object "phy" is not of class "phylo"') + attr(value, "formula") <- form + attr(value, "fixed") <- fixed + attr(value, "tree") <- phy + class(value) <- c("corMartins", "corPhyl", "corStruct") + value } -corGrafen <- function(value, phy, form=~1, fixed=FALSE) +corGrafen <- function(value, phy, form = ~1, fixed = FALSE) { - if(length(value) > 1) - stop("only one parameter is allowed in corGrafen structure.") - if(value < 0) stop("parameter rho must be positive.") - value <- log(value) # Optimization under constraint, use exponential transform. - if (!("phylo" %in% class(phy))) - stop("object \"phy\" is not of class \"phylo\"") - attr(value, "formula") <- form - attr(value, "fixed") <- fixed - attr(value, "tree") <- phy - class(value) <- c("corGrafen", "corPhyl", "corStruct") - return(value) + if (length(value) > 1) + stop("only one parameter is allowed") + if (value < 0) stop("parameter rho must be positive") + value <- log(value) # Optimization under constraint, use exponential transform. + if (!inherits(phy, "phylo")) + stop('object "phy" is not of class "phylo"') + attr(value, "formula") <- form + attr(value, "fixed") <- fixed + attr(value, "tree") <- phy + class(value) <- c("corGrafen", "corPhyl", "corStruct") + value } Initialize.corPhyl <- function(object, data, ...) { - # The same as in Initialize corStruct: - form <- formula(object) - ## Obtaining the group information, if any - if(!is.null(getGroupsFormula(form))) { - attr(object, "groups") <- getGroups(object, form, data = data) - attr(object, "Dim") <- Dim(object, attr(object, "groups")) - } else { # no groups - attr(object, "Dim") <- Dim(object, as.factor(rep(1, nrow(data)))) - } - ## Obtaining the covariate(s) - attr(object, "covariate") <- getCovariate(object, data = data) + ## The same as in Initialize corStruct: + form <- formula(object) + ## Obtaining the group information, if any + if (!is.null(getGroupsFormula(form))) { + attr(object, "groups") <- getGroups(object, form, data = data) + attr(object, "Dim") <- Dim(object, attr(object, "groups")) + } else { # no groups + attr(object, "Dim") <- Dim(object, as.factor(rep(1, nrow(data)))) + } + ## Obtaining the covariate(s) + attr(object, "covariate") <- getCovariate(object, data = data) - # Specific to corPhyl: - phy <- attr(object, "tree") - if (is.null(data)) - data <- parent.frame() - ## Added by EP 29 May 2006: - if (nrow(data) != length(phy$tip.label)) - stop("number of observations and number of tips in the tree are not equal.") - ## END - if(is.null(rownames(data))) { - warning("No row names supplied in dataframe, data taken to be in the same order as in tree.") - attr(object, "index") <- 1:dim(data)[1] - } else { - index <- match(rownames(data), phy$tip.label) - if(any(is.na(index))) { - warning("Row names in dataframe do not match tree tip names. data taken to be in the same order as in tree.") - attr(object, "index") <- 1:dim(data)[1] + ## Specific to corPhyl: + phy <- attr(object, "tree") + if (is.null(data)) data <- parent.frame() + ## Added by EP 29 May 2006: + if (nrow(data) != length(phy$tip.label)) + stop("number of observations and number of tips in the tree are not equal.") + ## END + if (is.null(rownames(data))) { + warning("No rownames supplied in data frame, data taken to be in the same order than in tree") + attr(object, "index") <- 1:dim(data)[1] } else { - attr(object, "index") <- index + index <- match(rownames(data), phy$tip.label) + if (any(is.na(index))) { + warning("Rownames in data frame do not match tree tip names; data taken to be in the same order as in tree") + attr(object, "index") <- 1:dim(data)[1] + } else { + attr(object, "index") <- index + } } - } - return(object) + object } -corMatrix.corBrownian <- function(object, covariate = getCovariate(object), corr = TRUE, ...) +corMatrix.corBrownian <- + function(object, covariate = getCovariate(object), corr = TRUE, ...) { - if (!("corBrownian" %in% class(object))) - stop("object is not of class \"corBrownian\".") - if(!any(attr(object, "index"))) - stop("object have not been initialized.") - tree <- attr(object, "tree") - mat <- vcv.phylo(tree, cor = corr) - n <- dim(mat)[1] - # reorder matrix: - matr <- matrix(nrow=n, ncol=n) - index <- attr(object, "index") - for(i in 1:n) - for(j in i:n) - matr[i,j] <- matr[j,i] <- mat[index[i], index[j]] - return(matr) + if (!("corBrownian" %in% class(object))) + stop('object is not of class "corBrownian"') + if (!any(attr(object, "index"))) + stop("object has not been initialized.") + tree <- attr(object, "tree") + mat <- vcv.phylo(tree, corr = corr) + n <- dim(mat)[1] + ## reorder matrix: + index <- attr(object, "index") + mat[index, index] } corMatrix.corMartins <- function(object, covariate = getCovariate(object), corr = TRUE, ...) { if (!("corMartins" %in% class(object))) - stop("object is not of class \"corMartins\".") - if(!any(attr(object, "index"))) - stop("object have not been initialized.") + stop('object is not of class "corMartins"') + if (!any(attr(object, "index"))) + stop("object has not been initialized.") tree <- attr(object, "tree") dist <- cophenetic.phylo(tree) mat <- exp(-object[1] * dist) if (corr) mat <- cov2cor(mat) n <- dim(mat)[1] ## reorder matrix: - matr <- matrix(nrow=n, ncol=n) index <- attr(object, "index") - for(i in 1:n) - for(j in i:n) - matr[i,j] <- matr[j,i] <- mat[index[i], index[j]] - matr + mat[index, index] } -corMatrix.corGrafen <- function(object, covariate = getCovariate(object), corr = TRUE, ...) +corMatrix.corGrafen <- + function(object, covariate = getCovariate(object), corr = TRUE, ...) { - if (!("corGrafen" %in% class(object))) - stop("object is not of class \"corGrafen\".") - if (!any(attr(object, "index"))) - stop("object have not been initialized.") - tree <- compute.brlen(attr(object, "tree"), - method = "Grafen", power = exp(object[1])) - mat <- vcv.phylo(tree, cor = corr) - n <- dim(mat)[1] - # reorder matrix: - matr <- matrix(nrow=n, ncol=n) - index <- attr(object, "index") - for(i in 1:n) - for(j in i:n) - matr[i,j] <- matr[j,i] <- mat[index[i], index[j]] - return(matr) + if (!("corGrafen" %in% class(object))) + stop('object is not of class "corGrafen"') + if (!any(attr(object, "index"))) + stop("object has not been initialized.") + tree <- compute.brlen(attr(object, "tree"), + method = "Grafen", power = exp(object[1])) + mat <- vcv.phylo(tree, corr = corr) + n <- dim(mat)[1] + ## reorder matrix: + index <- attr(object, "index") + mat[index, index] } coef.corBrownian <- function(object, unconstrained = TRUE, ...) { - if (!("corBrownian" %in% class(object))) - stop("object is not of class \"corBrownian\".") - return(numeric(0)) + if (!("corBrownian" %in% class(object))) + stop('object is not of class "corBrownian"') + numeric(0) } coef.corMartins <- function(object, unconstrained = TRUE, ...) { - if (!("corMartins" %in% class(object))) - stop("object is not of class \"corMartins\".") - if(unconstrained) { - if(attr(object, "fixed")) { - return(numeric(0)) - } else { - return(as.vector(object)) + if (!("corMartins" %in% class(object))) + stop('object is not of class "corMartins"') + if (unconstrained) { + if (attr(object, "fixed")) { + return(numeric(0)) + } else { + return(as.vector(object)) + } } - } - aux <- as.vector(object) - names(aux) <- "alpha" - return(aux) + aux <- as.vector(object) + names(aux) <- "alpha" + aux } coef.corGrafen <- function(object, unconstrained = TRUE, ...) { - if (!("corGrafen" %in% class(object))) - stop("object is not of class \"corGrafen\".") - if(unconstrained) { - if(attr(object, "fixed")) { - return(numeric(0)) - } else { - return(as.vector(object)) + if (!("corGrafen" %in% class(object))) + stop('object is not of class "corGrafen"') + if (unconstrained) { + if (attr(object, "fixed")) { + return(numeric(0)) + } else { + return(as.vector(object)) + } } - } - aux <- exp(as.vector(object)) - names(aux) <- "rho" - return(aux) + aux <- exp(as.vector(object)) + names(aux) <- "rho" + aux } ### removed node.sons() and node.leafnumber() (2006-10-12) @@ -219,7 +208,7 @@ corPagel <- function(value, phy, form = ~1, fixed = FALSE) { if (value < 0 || value > 1) stop("the value of lambda must be between 0 and 1.") - if (!("phylo" %in% class(phy))) + if (!inherits(phy, "phylo")) stop('object "phy" is not of class "phylo"') attr(value, "formula") <- form attr(value, "fixed") <- fixed @@ -233,7 +222,7 @@ corMatrix.corPagel <- { if (!any(attr(object, "index"))) stop("object has not been initialized") - mat <- vcv.phylo(attr(object, "tree"), cor = corr) + mat <- vcv.phylo(attr(object, "tree"), corr = corr) index <- attr(object, "index") mat <- mat[index, index] tmp <- diag(mat) @@ -257,7 +246,7 @@ corBlomberg <- function(value, phy, form = ~1, fixed = FALSE) { if (value <= 0) stop("the value of g must be greater than 0.") - if (!("phylo" %in% class(phy))) + if (!inherits(phy, "phylo")) stop('object "phy" is not of class "phylo"') attr(value, "formula") <- form attr(value, "fixed") <- fixed @@ -278,7 +267,7 @@ probably need to set 'fixed = TRUE' in corBlomberg().") phy <- attr(object, "tree") d <- (dist.nodes(phy)[length(phy$tip.label) + 1, ])^(1/object[1]) phy$edge.length <- d[phy$edge[, 2]] - d[phy$edge[, 1]] - mat <- vcv.phylo(phy, cor = corr) + mat <- vcv.phylo(phy, corr = corr) mat[index, index] } diff --git a/R/ace.R b/R/ace.R index b51b58e..ed86fc6 100644 --- a/R/ace.R +++ b/R/ace.R @@ -12,7 +12,7 @@ ace <- function(x, phy, type = "continuous", method = "ML", CI = TRUE, scaled = TRUE, kappa = 1, corStruct = NULL, ip = 0.1) { if (!inherits(phy, "phylo")) - stop('object "phy" is not of class "phylo".') + stop('object "phy" is not of class "phylo"') if (is.null(phy$edge.length)) stop("tree has no branch lengths") type <- match.arg(type, c("continuous", "discrete")) diff --git a/R/chronopl.R b/R/chronopl.R index 0e72a3f..98f27fb 100644 --- a/R/chronopl.R +++ b/R/chronopl.R @@ -1,8 +1,8 @@ -## chronopl.R (2011-07-04) +## chronopl.R (2012-02-09) ## Molecular Dating With Penalized Likelihood -## Copyright 2005-2011 Emmanuel Paradis +## Copyright 2005-2012 Emmanuel Paradis ## This file is part of the R-package `ape'. ## See the file ../COPYING for licensing issues. @@ -13,8 +13,8 @@ chronopl <- CV = FALSE, eval.max = 500, iter.max = 500, ...) { n <- length(phy$tip.label) - ROOT <- n + 1 - if (length(node) == 1 && node == "root") node <- ROOT + ROOT <- n + 1L + if (identical(node, "root")) node <- ROOT if (any(node <= n)) stop("node numbers should be greater than the number of tips") zerobl <- which(phy$edge.length <= 0) @@ -42,8 +42,9 @@ chronopl <- } m <- phy$Nnode el <- phy$edge.length - e <- phy$edge - N <- dim(e)[1] + e1 <- phy$edge[, 1L] + e2 <- phy$edge[, 2L] + N <- length(e1) TIPS <- 1:n EDGES <- 1:N @@ -52,7 +53,7 @@ chronopl <- ## `basal' contains the indices of the basal edges ## (ie, linked to the root): - basal <- which(e[, 1] == ROOT) + basal <- which(e1 == ROOT) Nbasal <- length(basal) ## `ind' contains in its 1st column the index of all nonbasal @@ -68,43 +69,54 @@ chronopl <- ind <- matrix(0L, N - Nbasal, 2) ind[, 1] <- EDGES[-basal] - ind[, 2] <- match(e[EDGES[-basal], 1], e[, 2]) + ind[, 2] <- match(e1[EDGES[-basal]], e2) age <- numeric(n + m) - ##ini.time <- node.depth(phy)[-TIPS] - 1 - ini.time <- node.depth(phy) - 1 +############################################################################# +### This bit sets 'ini.time' and should result in no negative branch lengths + + seq.nod <- .Call("seq_root2tip", phy$edge, n, phy$Nnode, PACKAGE = "ape") + + ini.time <- age + ini.time[ROOT:(n + m)] <- NA + ini.time[node] <- if (is.null(age.max)) age.min else (age.min + age.max) / 2 + + ## if no age given for the root, find one approximately: + if (is.na(ini.time[ROOT])) + ini.time[ROOT] <- if (is.null(age.max)) 3 * max(age.min) else 3 * max(age.max) + + ISnotNA.ALL <- unlist(lapply(seq.nod, function(x) sum(!is.na(ini.time[x])))) + o <- order(ISnotNA.ALL, decreasing = TRUE) + + for (y in seq.nod[o]) { + ISNA <- is.na(ini.time[y]) + if (any(ISNA)) { + i <- 2L # we know the 1st value is not NA, so we start at the 2nd one + while (i <= length(y)) { + if (ISNA[i]) { # we stop at the next NA + j <- i + 1L + while (ISNA[j]) j <- j + 1L # look for the next non-NA + nb.val <- j - i + by <- (ini.time[y[i - 1L]] - ini.time[y[j]]) / (nb.val + 1) + ini.time[y[i:(j - 1L)]] <- ini.time[y[i - 1L]] - by * seq_len(nb.val) + i <- j + 1L + } else i <- i + 1L + } + } + } - ## first, rescale all times with respect to - ## the age of the 1st element of `node': - ratio <- age.min[1]/ini.time[node[1]] - ini.time <- ini.time*ratio + real.edge.length <- ini.time[e1] - ini.time[e2] - ## because if (!is.null(age.max)), 'node' is modified, so we copy it in case CV = TRUE: + if (any(real.edge.length <= 0)) + stop("some initial branch lengths are zero or negative; + maybe you need to adjust the given dates -- see '?chronopl' for details") +############################################################################# + + ## because if (!is.null(age.max)), 'node' is modified, + ## so we copy it in case CV = TRUE: node.bak <- node - if (length(node) > 1) { - ini.time[node] <- age.min - real.edge.length <- ini.time[e[, 1]] - ini.time[e[, 2]] - while (any(real.edge.length <= 0)) { - for (i in EDGES) { - if (real.edge.length[i] > 0) next - if (e[i, 1] %in% node) { - ini.time[e[i, 2]] <- ini.time[e[1, 2]] - 2 * real.edge.length[i] - next - } - if (e[i, 2] %in% node) { - ini.time[e[i, 1]] <- ini.time[e[1, 1]] + 2 * real.edge.length[i] - next - } - ##browser() - ini.time[e[i, 2]] <- ini.time[e[1, 2]] - real.edge.length[i] - ini.time[e[i, 1]] <- ini.time[e[1, 1]] + real.edge.length[i] - } - real.edge.length <- ini.time[e[, 1]] - ini.time[e[, 2]] - ##print(min(real.edge.length)) - } - } ## `unknown.ages' will contain the index of the nodes of unknown age: unknown.ages <- n + 1:m @@ -115,13 +127,15 @@ chronopl <- if (!is.null(age.max)) { # are some nodes known within some intervals? lower[node - n] <- age.min upper[node - n] <- age.max + ## find nodes known within an interval: interv <- which(age.min != age.max) + ## drop them from the 'node' since they will be estimated: node <- node[-interv] - if (length(node)) age[node] <- age.min[-interv] + if (length(node)) age[node] <- age.min[-interv] # update 'age' } else age[node] <- age.min if (length(node)) { - unknown.ages <- unknown.ages[n - node] + unknown.ages <- unknown.ages[n - node] # 'n - node' is simplification for '-(node - n)' lower <- lower[n - node] upper <- upper[n - node] } @@ -137,7 +151,7 @@ chronopl <- minusploglik.gr <- function(rate, node.time) { grad <- numeric(N + length(unknown.ages)) age[unknown.ages] <- node.time - real.edge.length <- age[e[, 1]] - age[e[, 2]] + real.edge.length <- age[e1] - age[e2] if (any(real.edge.length < 0)) { grad[] <- 0 return(grad) @@ -158,7 +172,7 @@ chronopl <- } for (i in EDGES) { - ii <- c(which(e[, 2] == e[i, 1]), which(e[, 1] == e[i, 2])) + ii <- c(which(e2 == e1[i]), which(e1 == e2[i])) if (!length(ii)) next grad[i] <- grad[i] + lambda*(2*length(ii)*rate[i] - 2*sum(rate[ii])) } @@ -166,11 +180,11 @@ chronopl <- ## gradient for the 'node times' for (i in 1:length(unknown.ages)) { nd <- unknown.ages[i] - ii <- which(e[, 1] == nd) + ii <- which(e1 == nd) grad[i + N] <- sum(rate[ii] - el[ii]/real.edge.length[ii])#, na.rm = TRUE) if (nd != ROOT) { - ii <- which(e[, 2] == nd) + ii <- which(e2 == nd) grad[i + N] <- grad[i + N] - rate[ii] + el[ii]/real.edge.length[ii] } @@ -180,7 +194,7 @@ chronopl <- minusploglik <- function(rate, node.time) { age[unknown.ages] <- node.time - real.edge.length <- age[e[, 1]] - age[e[, 2]] + real.edge.length <- age[e1] - age[e2] if (any(real.edge.length < 0)) return(1e50) B <- rate*real.edge.length loglik <- sum(-B + el*log(B) - lfactorial(el)) @@ -199,7 +213,7 @@ chronopl <- attr(phy, "message") <- out$message age[unknown.ages] <- out$par[-EDGES] if (CV) ophy <- phy - phy$edge.length <- age[e[, 1]] - age[e[, 2]] + phy$edge.length <- age[e1] - age[e2] if (CV) attr(phy, "D2") <- chronopl.cv(ophy, lambda, age.min, age.max, node.bak, n, S, tol, eval.max, iter.max, ...) diff --git a/R/is.ultrametric.R b/R/is.ultrametric.R index ae4c439..792684f 100644 --- a/R/is.ultrametric.R +++ b/R/is.ultrametric.R @@ -1,8 +1,8 @@ -## is.ultrametric.R (2009-05-10) +## is.ultrametric.R (2012-02-09) ## Test if a Tree is Ultrametric -## Copyright 2003-2009 Emmanuel Paradis +## Copyright 2003-2012 Emmanuel Paradis ## This file is part of the R-package `ape'. ## See the file ../COPYING for licensing issues. @@ -10,21 +10,27 @@ is.ultrametric <- function(phy, tol = .Machine$double.eps^0.5) { if (!inherits(phy, "phylo")) - stop('object "phy" is not of class "phylo".') + stop('object "phy" is not of class "phylo"') if (is.null(phy$edge.length)) - stop("the tree has no branch lengths.") + stop("the tree has no branch lengths") + phy <- reorder(phy) n <- length(phy$tip.label) - n.node <- phy$Nnode + e1 <- phy$edge[, 1] + e2 <- phy$edge[, 2] + EL <- phy$edge.length + + ## xx: vecteur donnant la distance d'un noeud + ## ou d'un tip à partir de la racine + xx <- numeric(n + phy$Nnode) - ## xx: vecteur donnant la distance d'un - ## noeud ou tip à partir de la racine - xx <- numeric(n + n.node) + ## the following must start at the root and follow the + ## edges contiguously; so the tree must be either in cladewise + ## order (or in pruningwise but the for loop must start from + ## the bottom of the edge matrix) - for (i in 1:dim(phy$edge)[1]) - xx[phy$edge[i, 2]] <- xx[phy$edge[i, 1]] + phy$edge.length[i] + for (i in seq_len(length(e1))) + xx[e2[i]] <- xx[e1[i]] + EL[i] - if (identical(all.equal.numeric(var(xx[1:n]), - 0, tolerance = tol), TRUE)) TRUE - else FALSE + isTRUE(all.equal.numeric(var(xx[1:n]), 0, tolerance = tol)) } diff --git a/R/nodelabels.R b/R/nodelabels.R index 40bfab1..152761d 100644 --- a/R/nodelabels.R +++ b/R/nodelabels.R @@ -1,8 +1,8 @@ -## nodelabels.R (2010-07-17) +## nodelabels.R (2012-02-10) ## Labelling Trees -## Copyright 2004-2010 Emmanuel Paradis, 2006 Ben Bolker, and 2006 Jim Lemon +## Copyright 2004-2012 Emmanuel Paradis, 2006 Ben Bolker, and 2006 Jim Lemon ## This file is part of the R-package `ape'. ## See the file ../COPYING for licensing issues. @@ -178,7 +178,7 @@ edgelabels <- function(text, edge, adj = c(0.5, 0.5), frame = "rect", pch = NULL, thermo = NULL, pie = NULL, piecol = NULL, col = "black", bg = "lightgreen", horiz = FALSE, - width = NULL, height = NULL, ...) + width = NULL, height = NULL, date = NULL, ...) { lastPP <- get("last_plot.phylo", envir = .PlotPhyloEnv) if (missing(edge)) { @@ -200,6 +200,11 @@ edgelabels <- XX <- (lastPP$xx[subedge[, 1]] + lastPP$xx[subedge[, 2]]) / 2 YY <- (lastPP$yy[subedge[, 1]] + lastPP$yy[subedge[, 2]]) / 2 } + + ## suggestion by Rob Lanfear: + if (!is.null(date)) + XX[] <- max(lastPP$xx) - date + BOTHlabels(text, sel, XX, YY, adj, frame, pch, thermo, pie, piecol, col, bg, horiz, width, height, ...) } diff --git a/R/rTrait.R b/R/rTrait.R index 12dc5bf..1b5852d 100644 --- a/R/rTrait.R +++ b/R/rTrait.R @@ -1,8 +1,8 @@ -## rTrait.R (2011-06-15) +## rTrait.R (2012-02-09) ## Trait Evolution -## Copyright 2010-2011 Emmanuel Paradis +## Copyright 2010-2012 Emmanuel Paradis ## This file is part of the R-package `ape'. ## See the file ../COPYING for licensing issues. @@ -64,7 +64,7 @@ rTraitDisc <- diag(Q) <- -rowSums(Q) for (i in N:1) { p <- matexpo(Q * el[i])[x[anc[i]], ] - x[des[i]] <- .Internal(sample(k, size = 1, FALSE, prob = p)) + x[des[i]] <- sample.int(k, size = 1, FALSE, prob = p) } } diff --git a/R/read.nexus.R b/R/read.nexus.R index 642bee7..0268e68 100644 --- a/R/read.nexus.R +++ b/R/read.nexus.R @@ -1,8 +1,8 @@ -## read.nexus.R (2011-03-26) +## read.nexus.R (2012-02-09) ## Read Tree File in Nexus Format -## Copyright 2003-2011 Emmanuel Paradis and 2010 Klaus Schliep +## Copyright 2003-2012 Emmanuel Paradis and 2010 Klaus Schliep ## This file is part of the R-package `ape'. ## See the file ../COPYING for licensing issues. @@ -247,9 +247,5 @@ read.nexus <- function(file, tree.names = NULL) class(trees) <- "multiPhylo" if (!all(nms.trees == "")) names(trees) <- nms.trees } - if (length(grep("[\\/]", file)) == 1) - if (!file.exists(file)) # suggestion by Francois Michonneau - file <- paste(getwd(), file, sep = "/") - attr(trees, "origin") <- file trees } diff --git a/R/rtree.R b/R/rtree.R index 46f9abd..0df3596 100644 --- a/R/rtree.R +++ b/R/rtree.R @@ -1,8 +1,8 @@ -## rtree.R (2010-03-09) +## rtree.R (2012-02-09) ## Generates Trees -## Copyright 2004-2010 Emmanuel Paradis +## Copyright 2004-2012 Emmanuel Paradis ## This file is part of the R-package `ape'. ## See the file ../COPYING for licensing issues. @@ -10,7 +10,7 @@ rtree <- function(n, rooted = TRUE, tip.label = NULL, br = runif, ...) { foo <- function(n, pos) { - n1 <- .Internal(sample(n - 1, 1, FALSE, NULL)) + n1 <- sample.int(n - 1, 1, FALSE, NULL) n2 <- n - n1 po2 <- pos + 2*n1 - 1 edge[c(pos, po2), 1] <<- nod @@ -54,11 +54,11 @@ rtree <- function(n, rooted = TRUE, tip.label = NULL, br = runif, ...) i <- which(is.na(edge[, 2])) edge[i, 2] <- 1:n } else { # n > 4 - n1 <- .Internal(sample(n - 2, 1, FALSE, NULL)) + n1 <- sample.int(n - 2, 1, FALSE, NULL) if (n1 == n - 2) { n2 <- n3 <- 1 } else { - n2 <- .Internal(sample(n - n1 - 1, 1, FALSE, NULL)) + n2 <- sample.int(n - n1 - 1, 1, FALSE, NULL) n3 <- n - n1 - n2 } po2 <- 2*n1 @@ -207,6 +207,6 @@ stree <- function(n, type = "star", tip.label = NULL) tip.label <- paste("t", 1:n, sep = "") phy <- list(edge = edge, tip.label = tip.label, Nnode = m) class(phy) <- "phylo" - attr(phy, "order" <- "cladewise") + attr(phy, "order") <- "cladewise" phy } diff --git a/R/vcv.phylo.R b/R/vcv.phylo.R index 2473280..294b767 100644 --- a/R/vcv.phylo.R +++ b/R/vcv.phylo.R @@ -1,8 +1,8 @@ -## vcv.phylo.R (2009-11-19) +## vcv.phylo.R (2012-02-09) ## Phylogenetic Variance-Covariance or Correlation Matrix -## Copyright 2002-2009 Emmanuel Paradis +## Copyright 2002-2012 Emmanuel Paradis ## This file is part of the R-package `ape'. ## See the file ../COPYING for licensing issues. @@ -12,39 +12,41 @@ vcv <- function(phy, ...) UseMethod("vcv") vcv.phylo <- function(phy, model = "Brownian", corr = FALSE, ...) { if (is.null(phy$edge.length)) - stop("the tree has no branch lengths") - - foo <- function(node, var, endofclade) { - ## First, get the extent of clade descending - ## from `node' in the matrix `edge': - from <- which(phy$edge[, 1] == node) - to <- c(from[-1] - 1, endofclade) - ## Get the #'s of the descendants of `node': - desc <- phy$edge[from, 2] - ## The variance of each of these is easy: - vcv[desc, desc] <<- var + phy$edge.length[from] - nd <- length(desc) - ## The variance of `node' is equal to the covariance of - ## each possible pair among its descendant clades. - for (i in 1:(nd - 1)) - for (j in (i + 1):nd) - for (k in phy$edge[from[i]:to[i], 2]) - for (l in phy$edge[from[j]:to[j], 2]) - vcv[k, l] <<- vcv[l, k] <<- var - for (i in 1:nd) { - if (desc[i] <= n) next - foo(desc[i], vcv[desc[i], desc[i]], to[i]) + stop("the tree has no branch lengths") + + pp <- prop.part(phy) + + phy <- reorder(phy, "pruningwise") + + n <- length(phy$tip.label) + e1 <- phy$edge[, 1] + e2 <- phy$edge[, 2] + EL <- phy$edge.length + + ## xx: vecteur donnant la distance d'un noeud + ## ou d'un tip à partir de la racine + ## (same than in is.ultrametric) + xx <- numeric(n + phy$Nnode) + + vcv <- matrix(0, n, n) + + ## the loop below starts from the bottom of the edge matrix, so + ## from the root + + for (i in length(e1):1) { + var.cur.node <- xx[e1[i]] + xx[e2[i]] <- var.cur.node + EL[i] # sets the variance + j <- i - 1L + while (e1[j] == e1[i] && j > 0) { + left <- if (e2[j] > n) pp[[e2[j] - n]] else e2[j] + right <- if (e2[i] > n) pp[[e2[i] - n]] else e2[i] + vcv[left, right] <- vcv[right, left] <- var.cur.node # sets the covariance + j <- j - 1L } } - phy <- reorder(phy) - n <- length(phy$tip.label) - n.node <- phy$Nnode - N <- n.node + n - vcv <- matrix(0, N, N) - foo(n + 1, 0, dim(phy$edge)[1]) + diag(vcv) <- xx[1:n] - vcv <- vcv[1:n, 1:n] if (corr) { ## This is inspired from the code of `cov2cor' (2005-09-08): M <- vcv @@ -52,7 +54,9 @@ vcv.phylo <- function(phy, model = "Brownian", corr = FALSE, ...) vcv[] <- Is * M * rep(Is, each = n) vcv[1 + 0:(n - 1)*(n + 1)] <- 1 } - rownames(vcv) <- colnames(vcv) <- phy$tip.label + + dimnames(vcv)[1:2] <- list(phy$tip.label) + vcv } diff --git a/R/write.nexus.R b/R/write.nexus.R index 1db2252..fd44078 100644 --- a/R/write.nexus.R +++ b/R/write.nexus.R @@ -1,13 +1,13 @@ -## write.nexus.R (2011-03-26) +## write.nexus.R (2012-02-09) ## Write Tree File in Nexus Format -## Copyright 2003-2011 Emmanuel Paradis +## Copyright 2003-2012 Emmanuel Paradis ## This file is part of the R-package `ape'. ## See the file ../COPYING for licensing issues. -write.nexus <- function(..., file = "", translate = TRUE, original.data = TRUE) +write.nexus <- function(..., file = "", translate = TRUE, original.data = NULL) { obj <- list(...) ## We insure that all trees are in a list, even if there is a single one: @@ -21,40 +21,21 @@ write.nexus <- function(..., file = "", translate = TRUE, original.data = TRUE) cat("#NEXUS\n", file = file) cat(paste("[R-package APE, ", date(), "]\n\n", sep = ""), file = file, append = TRUE) - if (original.data) { - if (!is.null(attr(obj[[1]], "origin"))) { - if (!file.exists(attr(obj[[1]], "origin"))) { - warning(paste("the file", attr(obj[[1]], "origin"), - "cannot be found, -the original data won't be written with the tree.")) - original.data <- FALSE - } - else { - ORI <- scan(file = attr(obj[[1]], "origin"), what = character(), - sep = "\n", skip = 1) - start <- grep("BEGIN TAXA;", ORI) - ORI <- ORI[-(1:(start - 1))] - ORI <- gsub("ENDBLOCK;", "END;", ORI) - endblock <- grep("END;", ORI) - start <- grep("BEGIN TREES;", ORI) - end <- endblock[endblock > start][1] - cat(ORI[1:(start - 1)], file = file, append = TRUE, sep = "\n") - ORI <- ORI[-(1:end)] - } - } - else original.data <- FALSE - } + + if (!is.null(original.data)) + warning("the option 'original.data' is deprecated and will be removed soon. Please update your code.") + N <- length(obj[[1]]$tip.label) - if (!original.data) { - cat("BEGIN TAXA;\n", file = file, append = TRUE) - cat(paste("\tDIMENSIONS NTAX = ", N, ";\n", sep = ""), - file = file, append = TRUE) - cat("\tTAXLABELS\n", file = file, append = TRUE) - cat(paste("\t\t", obj[[1]]$tip.label, sep = ""), - sep = "\n", file = file, append = TRUE) - cat("\t;\n", file = file, append = TRUE) - cat("END;\n", file = file, append = TRUE) - } + + cat("BEGIN TAXA;\n", file = file, append = TRUE) + cat(paste("\tDIMENSIONS NTAX = ", N, ";\n", sep = ""), + file = file, append = TRUE) + cat("\tTAXLABELS\n", file = file, append = TRUE) + cat(paste("\t\t", obj[[1]]$tip.label, sep = ""), + sep = "\n", file = file, append = TRUE) + cat("\t;\n", file = file, append = TRUE) + cat("END;\n", file = file, append = TRUE) + cat("BEGIN TREES;\n", file = file, append = TRUE) if (translate) { cat("\tTRANSLATE\n", file = file, append = TRUE) @@ -93,5 +74,4 @@ the original data won't be written with the tree.")) "\n", sep = "", file = file, append = TRUE) } cat("END;\n", file = file, append = TRUE) - if(original.data) cat(ORI, file = file, append = TRUE, sep = "\n") } diff --git a/man/chronopl.Rd b/man/chronopl.Rd index 6a5366c..13f183f 100644 --- a/man/chronopl.Rd +++ b/man/chronopl.Rd @@ -1,6 +1,12 @@ \name{chronopl} \alias{chronopl} \title{Molecular Dating With Penalized Likelihood} +\description{ + This function estimates the node ages of a tree using a + semi-parametric method based on penalized likelihood (Sanderson + 2002). The branch lengths of the input tree are interpreted as mean + numbers of substitutions (i.e., per site). +} \usage{ chronopl(phy, lambda, age.min = 1, age.max = NULL, node = "root", S = 1, tol = 1e-8, @@ -27,12 +33,6 @@ chronopl(phy, lambda, age.min = 1, age.max = NULL, algorithm.} \item{\dots}{further arguments passed to control \code{nlminb}.} } -\description{ - This function estimates the node ages of a tree using a - semi-parametric method based on penalized likelihood (Sanderson - 2002). The branch lengths of the input tree are interpreted as mean - numbers of substitutions (i.e., per site). -} \details{ The idea of this method is to use a trade-off between a parametric formulation where each branch has its own rate, and a nonparametric @@ -58,6 +58,15 @@ chronopl(phy, lambda, age.min = 1, age.max = NULL, c(10, 6), node = c(15, 18)} means that the age of node 15 is 10 units of time, and the age of node 18 is between 5 and 6). + If two nodes are linked (i.e., one is the ancestor of the other) and + have the same values of \code{age.min} and \code{age.max} (say, 10 and + 15) this will result in an error because the medians of these values + are used as initial times (here 12.5) giving initial branch length(s) + equal to zero. The easiest way to solve this is to change slightly the + given values, for instance use \code{age.max = 14.9} for the youngest + node, or \code{age.max = 15.1} for the oldest one (or similarly for + \code{age.min}). + The input tree may have multichotomies. If some internal branches are of zero-length, they are collapsed (with a warning), and the returned tree will have less nodes than the input one. The presence of diff --git a/man/nodelabels.Rd b/man/nodelabels.Rd index 9a7f51a..c499694 100644 --- a/man/nodelabels.Rd +++ b/man/nodelabels.Rd @@ -20,7 +20,7 @@ tiplabels(text, tip, adj = c(0.5, 0.5), frame = "rect", edgelabels(text, edge, adj = c(0.5, 0.5), frame = "rect", pch = NULL, thermo = NULL, pie = NULL, piecol = NULL, col = "black", bg = "lightgreen", horiz = FALSE, - width = NULL, height = NULL, ...) + width = NULL, height = NULL, date = NULL, ...) } \arguments{ @@ -63,6 +63,8 @@ edgelabels(text, edge, adj = c(0.5, 0.5), frame = "rect", \item{horiz, width, height}{parameters controlling the aspect of thermometers; by default, their width and height are determined automatically.} + \item{date}{specifies the positions of labels on edges of chronograms + with respect to the time scale.} } \details{ These three functions have the same optional arguments and the same diff --git a/man/pcoa.Rd b/man/pcoa.Rd index 04933bf..8a9baac 100644 --- a/man/pcoa.Rd +++ b/man/pcoa.Rd @@ -88,6 +88,7 @@ Lingoes, J. C. (1971) Some boundary conditions for a monotone analysis of symmet \examples{ # Oribatid mite data from Borcard and Legendre (1994) +\dontrun{ if (require(vegan)) { data(mite) # Community composition data, 70 peat cores, 35 species @@ -113,5 +114,6 @@ biplot(res, mite.log, dir.axis1=-1, dir.axis2=-1) biplot(res, mite.log.st, dir.axis1=-1, dir.axis2=-1) } } +} \keyword{ multivariate } diff --git a/man/rTraitDisc.Rd b/man/rTraitDisc.Rd index e2c43d9..a800c46 100644 --- a/man/rTraitDisc.Rd +++ b/man/rTraitDisc.Rd @@ -72,8 +72,17 @@ data(bird.orders) ### the two followings are the same: rTraitDisc(bird.orders) rTraitDisc(bird.orders, model = matrix(c(0, 0.1, 0.1, 0), 2)) + ### two-state model with irreversibility: rTraitDisc(bird.orders, model = matrix(c(0, 0, 0.1, 0), 2)) + +### simple two-state model: +tr <- rcoal(n <- 40, br = runif) +x <- rTraitDisc(tr, ancestor = TRUE) +plot(tr, show.tip.label = FALSE) +nodelabels(pch = 19, col = x[-(1:n)]) +tiplabels(pch = 19, col = x[1:n]) + ### an imaginary model with stasis 0.5 time unit after a node, then ### random evolution: foo <- function(x, l) { diff --git a/man/read.nexus.Rd b/man/read.nexus.Rd index 3c8d9c1..8a668d1 100644 --- a/man/read.nexus.Rd +++ b/man/read.nexus.Rd @@ -60,12 +60,6 @@ read.nexus(file, tree.names = NULL) If several trees are read in the file, the returned object is of class \code{"multiPhylo"}, and is a list of objects of class \code{"phylo"}. - - An attribute \code{"origin"} is further given to the returned object - which gives the name of the source file (with its path). This is used - to write a tree in a NEXUS file where all the original data must be - written (not only the tree) in accordance to the specifications of - Maddison et al. (1997). } \references{ Maddison, D. R., Swofford, D. L. and Maddison, W. P. (1997) NEXUS: an diff --git a/man/write.nexus.Rd b/man/write.nexus.Rd index a3e66aa..faf2ef6 100644 --- a/man/write.nexus.Rd +++ b/man/write.nexus.Rd @@ -2,7 +2,7 @@ \alias{write.nexus} \title{Write Tree File in Nexus Format} \usage{ -write.nexus(..., file = "", translate = TRUE, original.data = TRUE) +write.nexus(..., file = "", translate = TRUE, original.data = NULL) } \arguments{ \item{\dots}{either (i) a single object of class \code{"phylo"}, (ii) a @@ -14,18 +14,12 @@ write.nexus(..., file = "", translate = TRUE, original.data = TRUE) \item{translate}{a logical, if \code{TRUE} (the default) a translation of the tip labels is done which are replaced in the parenthetic representation with tokens.} - \item{original.data}{a logical, if \code{TRUE} (the default) the - data in the original NEXUS file are eventually written in - \code{"file"} (see below).} + \item{original.data}{deprecated; will be removed soon.} } \description{ This function writes trees in a file with the NEXUS format. } \details{ - If \code{original.data = TRUE}, the file as specified by the attribute - \code{"origin"} of the first tree is read and its data (except the - trees) are written in \code{file}. - If several trees are given, they must have all the same tip labels. If among the objects given some are not trees of class \code{"phylo"}, diff --git a/src/BIONJ.c b/src/BIONJ.c index f8bf683..d08d69d 100644 --- a/src/BIONJ.c +++ b/src/BIONJ.c @@ -1,4 +1,4 @@ -/* BIONJ.c 2008-01-14 */ +/* BIONJ.c 2012-02-09 */ /* Copyright 2007-2008 Olivier Gascuel, Hoa Sien Cuong, R port by Vincent Lefort, bionj() below modified by Emmanuel Paradis */ @@ -117,8 +117,7 @@ void Initialize(float **delta, double *X, char **labels, int n, POINTERS *trees) name=(WORD *)calloc(1,sizeof(WORD)); /* taxon name is */ if(name == NULL) /* put in trees */ { - Rprintf("Out of memories !!"); - exit(0); + error("out of memory"); } else { @@ -206,15 +205,13 @@ void bionj(double *X, int *N, char **labels, delta[i]=(float *)calloc(n+1, sizeof(float)); if(delta[i] == NULL) { - Rprintf("Out of memories!!"); - exit(0); + error("out of memory"); } } trees=(POINTERS *)calloc(n+1,sizeof(POINTERS)); if(trees == NULL) { - Rprintf("Out of memories!!"); - exit(0); + error("out of memory"); } /* initialise and symmetrize the running delta matrix */ r=n; @@ -372,8 +369,7 @@ void Concatenate(char chain1[MAX_LABEL_LENGTH], int ind, POINTERS *trees, int po bran=(WORD *)calloc(1,sizeof(WORD)); if(bran == NULL) { - Rprintf("Out of memories"); - exit(0); + error("out of memory"); } else { diff --git a/src/me.c b/src/me.c index 21a2396..6841802 100644 --- a/src/me.c +++ b/src/me.c @@ -1,4 +1,4 @@ -/* me.c 2008-01-14 */ +/* me.c 2012-02-09 */ /* Copyright 2007-2008 Olivier Gascuel, Rick Desper, R port by Vincent Lefort, me_*() below modified by Emmanuel Paradis */ @@ -274,8 +274,7 @@ tree *detrifurcate(tree *T) return(T); if (NULL != v->parentEdge) { - Rprintf ("Error: root %s is poorly rooted.\n",v->label); - exit(0); + error("root %s is poorly rooted.", v->label); } for(e = v->middleEdge, v->middleEdge = NULL; NULL != e; e = f ) { @@ -326,8 +325,7 @@ void compareSets(tree *T, set *S) } if (-1 == v->index2) { - Rprintf("Error leaf %s in tree not in distance matrix.\n",v->label); - exit(0); + error("leaf %s in tree not in distance matrix.", v->label); } e = depthFirstTraverse(T,NULL); while (NULL != e) @@ -335,16 +333,14 @@ void compareSets(tree *T, set *S) v = e->head; if ((leaf(v)) && (-1 == v->index2)) { - Rprintf("Error leaf %s in tree not in distance matrix.\n",v->label); - exit(0); + error("leaf %s in tree not in distance matrix.", v->label); } e = depthFirstTraverse(T,e); } for(X = S; NULL != X; X = X->secondNode) if (X->firstNode->index2 > -1) { - Rprintf("Error node %s in matrix but not a leaf in tree.\n",X->firstNode->label); - exit(0); + error("node %s in matrix but not a leaf in tree.", X->firstNode->label); } return; } diff --git a/src/me_ols.c b/src/me_ols.c index 3c002c9..d9ce297 100644 --- a/src/me_ols.c +++ b/src/me_ols.c @@ -234,8 +234,8 @@ void printDoubleTable(double **A, int d) for(i=0;imiddleEdge = thisEdge; else { - Rprintf("Error: node %s has too many (>3) children.\n",centerNode->label); - exit(0); + error("node %s has too many (>3) children.", centerNode->label); } //sprintf(thisEdge->label,"E%d",edgeCount++); //snprintf(thisEdge->label,MAX_LABEL_LENGTH,"E%d",edgeCount++); @@ -197,8 +195,7 @@ tree *readNewickString (char *str, int numLeaves) if ('(' != str[0]) { - Rprintf("Error reading generated tree - does not start with '('.\n"); - exit(0); + error("generated tree does not start with '('"); } inputLength = strlen (str)+1; for(i = 0; i < inputLength; i++) @@ -225,8 +222,7 @@ tree *readNewickString (char *str, int numLeaves) } else if (parCount < 0) { - Rprintf("Error reading generated tree. Too many right parentheses.\n"); - exit(0); + error("generated tree has too many right parentheses"); } } centerNode = decodeNewickSubtree (str, T, &uCount); -- 2.39.2