X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=R%2Frtree.R;h=46f9abd3c3cc435de4adc92fb0a32d1b73bbb1e7;hb=90f18c75d642f56b020bc6e0cdd0c5949c1d9a1d;hp=6d0172c1ef8b5029739dbd91eb7d027064f1e37f;hpb=2e3b8be33765d5311e071efb39871249d16f7e30;p=ape.git diff --git a/R/rtree.R b/R/rtree.R index 6d0172c..46f9abd 100644 --- a/R/rtree.R +++ b/R/rtree.R @@ -1,8 +1,8 @@ -## rtree.R (2008-01-08) +## rtree.R (2010-03-09) -## Generates Random Trees +## Generates Trees -## Copyright 2004-2008 Emmanuel Paradis +## Copyright 2004-2010 Emmanuel Paradis ## This file is part of the R-package `ape'. ## See the file ../COPYING for licensing issues. @@ -14,39 +14,39 @@ rtree <- function(n, rooted = TRUE, tip.label = NULL, br = runif, ...) 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 @@ -70,21 +70,21 @@ rtree <- function(n, rooted = TRUE, tip.label = NULL, br = runif, ...) foo(n1, 2) } else if (n1 == 2) { edge[2:3, 1] <- edge[1, 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 (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 + ## nod <- nod + 1L } i <- which(is.na(edge[, 2])) edge[i, 2] <- 1:n @@ -94,29 +94,35 @@ 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" 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,18 +131,19 @@ 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, ...) @@ -146,3 +153,60 @@ rmtree <- function(N, n, rooted = TRUE, tip.label = NULL, br = runif, ...) 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 +}