X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=R%2Frtree.R;h=b36f9d307a40c655cbaf85b2395fada04839d64b;hb=b0d1251527d8dd48ca1703b1cfaf217f413eda0e;hp=2e1a8a4778ba11d3b35859da41658e48f360b7a3;hpb=c827059eeafc8cbe41c812b26979543ab287803e;p=ape.git diff --git a/R/rtree.R b/R/rtree.R index 2e1a8a4..b36f9d3 100644 --- a/R/rtree.R +++ b/R/rtree.R @@ -1,8 +1,8 @@ -## rtree.R (2007-12-22) +## rtree.R (2013-04-02) -## Generates Random Trees +## Generates Trees -## Copyright 2004-2007 Emmanuel Paradis +## Copyright 2004-2013 Emmanuel Paradis ## This file is part of the R-package `ape'. ## See the file ../COPYING for licensing issues. @@ -10,43 +10,43 @@ 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 - nod <<- nod + 1 + nod <<- nod + 1L if (n1 > 2) { edge[pos, 2] <<- nod foo(n1, pos + 1) } else if (n1 == 2) { edge[c(pos + 1, pos + 2), 1] <<- edge[pos, 2] <<- nod - nod <<- nod + 1 + nod <<- nod + 1L } if (n2 > 2) { edge[po2, 2] <<- nod foo(n2, po2 + 1) } else if (n2 == 2) { edge[c(po2 + 1, po2 + 2), 1] <<- edge[po2, 2] <<- nod - nod <<- nod + 1 + nod <<- nod + 1L } } if (n < 2) stop("a tree must have at least 2 tips.") - nbr <- 2 * n - 2 - if (!rooted) nbr <- nbr - 1 + nbr <- 2 * n - 3 + rooted edge <- matrix(NA, nbr, 2) + n <- as.integer(n) if (n == 2) { - if (rooted) edge[] <- c(3, 3, 1, 2) + if (rooted) edge[] <- c(3L, 3L, 1L, 2L) else stop("an unrooted tree must have at least 3 tips.") } else if (n == 3) { edge[] <- - if (rooted) c(4, 5, 5, 4, 5, 1:3) - else c(4, 4, 4, 1:3) + if (rooted) c(4L, 5L, 5L, 4L, 5L, 1:3) + else c(4L, 4L, 4L, 1:3) } else if (n == 4 && !rooted) { - edge[] <- c(5, 6, 6, 5, 5, 6, 1:4) + edge[] <- c(5L, 6L, 6L, 5L, 5L, 6L, 1:4) } else { - nod <- n + 1 + nod <- n + 1L if (rooted) { # n > 3 foo(n, 1) ## The following is slightly more efficient than affecting the @@ -54,39 +54,39 @@ 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)) - if (n1 == n - 2) { - n2 <- n3 <- 1 + n1 <- sample.int(n - 2L, 1L) + if (n1 == n - 2L) { + n2 <- n3 <- 1L } else { - n2 <- .Internal(sample(n - n1 - 1, 1, FALSE, NULL)) + n2 <- sample.int(n - n1 - 1L, 1L) n3 <- n - n1 - n2 } - po2 <- 2*n1 - po3 <- 2*(n1 + n2) - 1 - edge[c(1, po2, po3), 1] <- nod - nod <- nod + 1 - if (n1 > 2) { - edge[1, 2] <- nod - foo(n1, 2) - } else if (n1 == 2) { - edge[2:3, 1] <- edge[1, 2] <- nod - nod <- nod + 1 + po2 <- 2L * n1 + po3 <- 2L * (n1 + n2) - 1L + edge[c(1, po2, po3), 1L] <- nod + nod <- nod + 1L + if (n1 > 2L) { + edge[1L, 2L] <- nod + foo(n1, 2L) + } else if (n1 == 2L) { + edge[2:3, 1L] <- edge[1L, 2L] <- nod + nod <- nod + 1L } - if (n2 > 2) { - edge[po2, 2] <- nod - foo(n2, po2 + 1) - } else if (n2 == 2) { - edge[c(po2 + 1, po2 + 2), 1] <- edge[po2, 2] <- nod - nod <- nod + 1 + if (n2 > 2L) { + edge[po2, 2L] <- nod + foo(n2, po2 + 1L) + } else if (n2 == 2L) { + edge[c(po2 + 1L, po2 + 2), 1L] <- edge[po2, 2L] <- nod + nod <- nod + 1L } if (n3 > 2) { - edge[po3, 2] <- nod - foo(n3, po3 + 1) - } else if (n3 == 2) { - edge[c(po3 + 1, po3 + 2), 1] <- edge[po3, 2] <- nod - ## nod <- nod + 1 + edge[po3, 2L] <- nod + foo(n3, po3 + 1L) + } else if (n3 == 2L) { + edge[c(po3 + 1L, po3 + 2), 1L] <- edge[po3, 2L] <- nod + ## nod <- nod + 1L } - i <- which(is.na(edge[, 2])) + i <- which(is.na(edge[, 2L])) edge[i, 2] <- 1:n } } @@ -94,29 +94,36 @@ rtree <- function(n, rooted = TRUE, tip.label = NULL, br = runif, ...) phy$tip.label <- if (is.null(tip.label)) paste("t", sample(n), sep = "") else sample(tip.label) - if (is.function(br)) phy$edge.length <- br(nbr, ...) - phy$Nnode <- if (rooted) n - 1 else n - 2 + if (!is.null(br)) { + phy$edge.length <- + if (is.function(br)) br(nbr, ...) else rep(br, length.out = nbr) + } + phy$Nnode <- n - 2L + rooted class(phy) <- "phylo" + attr(phy, "order") <- "cladewise" phy } -rcoal <- function(n, tip.label = NULL, br = rexp, ...) +rcoal <- function(n, tip.label = NULL, br = "coalescent", ...) { + n <- as.integer(n) nbr <- 2*n - 2 edge <- matrix(NA, nbr, 2) - x <- br(n - 1, ...) # coalescence times + ## coalescence times by default: + x <- if (is.character(br)) 2*rexp(n - 1)/(as.double(n:2) * as.double((n - 1):1)) + else if (is.numeric(br)) rep(br, length.out = n - 1) else br(n - 1, ...) if (n == 2) { - edge[] <- c(3, 3, 1:2) + edge[] <- c(3L, 3L, 1:2) edge.length <- rep(x, 2) } else if (n == 3) { - edge[] <- c(4, 5, 5, 4, 5, 1:3) - edge.length <- c(x[2], x[1], x[1], sum(x)) + edge[] <- c(4L, 5L, 5L, 4L, 5L, 1:3) + edge.length <- c(x[c(2, 1, 1)], sum(x)) } else { edge.length <- numeric(nbr) h <- numeric(2*n - 1) # initialized with 0's node.height <- cumsum(x) pool <- 1:n - nextnode <- 2*n - 1 + nextnode <- 2L*n - 1L for (i in 1:(n - 1)) { y <- sample(pool, size = 2) ind <- (i - 1)*2 + 1:2 @@ -125,16 +132,82 @@ rcoal <- function(n, tip.label = NULL, br = rexp, ...) edge.length[ind] <- node.height[i] - h[y] h[nextnode] <- node.height[i] pool <- c(pool[! pool %in% y], nextnode) - nextnode <- nextnode - 1 + nextnode <- nextnode - 1L } } phy <- list(edge = edge, edge.length = edge.length) - phy$tip.label <- - if (is.null(tip.label)) paste("t", 1:n, sep = "") - else tip.label - phy$Nnode <- n - 1 + if (is.null(tip.label)) + tip.label <- paste("t", 1:n, sep = "") + phy$tip.label <- sample(tip.label) + phy$Nnode <- n - 1L class(phy) <- "phylo" - ##reorder(phy) + phy <- reorder(phy) ## to avoid crossings when converting with as.hclust: - read.tree(text = write.tree(phy)) + phy$edge[phy$edge[, 2] <= n, 2] <- 1:n + phy +} + +rmtree <- function(N, n, rooted = TRUE, tip.label = NULL, br = runif, ...) +{ + a <- replicate(N, rtree(n, rooted = rooted, tip.label = tip.label, + br = br, ...), simplify = FALSE) + class(a) <- "multiPhylo" + a +} + +stree <- function(n, type = "star", tip.label = NULL) +{ + type <- match.arg(type, c("star", "balanced", "left", "right")) + n <- as.integer(n) + if (type == "star") { + N <- n + m <- 1L + } else { + m <- n - 1L + N <- n + m - 1L + } + edge <- matrix(0L, N, 2) + + switch(type, "star" = { + edge[, 1] <- n + 1L + edge[, 2] <- 1:n + }, "balanced" = { + if (log2(n) %% 1) + stop("'n' is not a power of 2: cannot make a balanced tree") + foo <- function(node, size) { + if (size == 2) { + edge[c(i, i + 1L), 1L] <<- node + edge[c(i, i + 1L), 2L] <<- c(nexttip, nexttip + 1L) + nexttip <<- nexttip + 2L + i <<- i + 2L + } else { + for (k in 1:2) { # do the 2 subclades + edge[i, ] <<- c(node, nextnode) + nextnode <<- nextnode + 1L + i <<- i + 1L + foo(nextnode - 1L, size/2) + } + } + } + i <- 1L + nexttip <- 1L + nextnode <- n + 2L + foo(n + 1L, n) + }, "left" = { + edge[c(seq.int(from = 1, to = N - 1, by = 2), N), 2L] <- 1:n + nodes <- (n + 1L):(n + m) + edge[seq.int(from = 2, to = N - 1, by = 2), 2L] <- nodes[-1] + edge[, 1L] <- rep(nodes, each = 2) + }, "right" = { + nodes <- (n + 1L):(n + m) + edge[, 1L] <- c(nodes, rev(nodes)) + edge[, 2L] <- c(nodes[-1], 1:n) + }) + + if (is.null(tip.label)) + tip.label <- paste("t", 1:n, sep = "") + phy <- list(edge = edge, tip.label = tip.label, Nnode = m) + class(phy) <- "phylo" + attr(phy, "order") <- "cladewise" + phy }