version 0.4

DESCRIPTION

@@ -1,19 +1,19 @@
 Package: coalescentMCMC
-Version: 0.3
-Date: 2013-10-27
+Version: 0.4
+Date: 2013-12-04
 Title: MCMC Algorithms for the Coalescent
 Authors@R: c(person("Emmanuel", "Paradis", role = c("aut", "cre", "cph"), email = "Emmanuel.Paradis@ird.fr"))
 Depends: ape, coda
-Imports: phangorn, stats
+Imports: Matrix, phangorn, stats
 ZipData: no
 Description: coalescentMCMC provides a flexible framework for
   coalescent analyses in R. It includes a main function running the
   MCMC algorithm, auxiliary functions for tree rearrangement, and some
   functions to compute population genetic parameters.
 License: GPL (>= 2)
-Packaged: 2013-10-27 04:34:21 UTC; paradis
+Packaged: 2013-12-04 08:23:47 UTC; paradis
 Author: Emmanuel Paradis [aut, cre, cph]
 Maintainer: Emmanuel Paradis <Emmanuel.Paradis@ird.fr>
 NeedsCompilation: yes
 Repository: CRAN
-Date/Publication: 2013-10-27 07:08:27
+Date/Publication: 2013-12-04 09:41:46

MD5

@@ -1,24 +1,24 @@
 0632cc7589c46bbf4981f0a778f1ef2a *COPYING
-055c18b2e6e0a17502d01cd1f894166b *DESCRIPTION
-4c0e9e2412b04458ac4fe1f24368b785 *NAMESPACE
-d9e6cdd8fe023cb932c9202efbc3de92 *NEWS
-0f81220bc1e0b8d81bec4efdc50214a6 *R/coalescentMCMC.R
-a52942fca58804291fa4cca2a5e9c3dd *R/dcoal.R
-5855c71436ebb147dfc8d5b038d500ef *R/treeOperators.R
-d3bc00b2363b4155824eff32034801ad *R/zzz.R
+8b2d664f670dceb8eb4b7939fe295c8e *DESCRIPTION
+1c9c9c190266b2db4ab3f7585abb4413 *NAMESPACE
+7baac206f2f3c02bafedd8c234570dbd *NEWS
+15406815fb4cbdc755f970201072fac0 *R/coalescentMCMC.R
+5f2eb5752f0e5d2e8109756e9489eb85 *R/dcoal.R
+c10bd794fbc949a2537ddd990f902ced *R/treeOperators.R
+d85e17220990885e7d898858088064ea *R/zzz.R
 01319ab1a68cc304cbb0593461e10541 *build/vignette.rds
 dfa3ae7f0e0bc0918691ff5894094c03 *inst/doc/CoalescentModels.R
-38240cef46b41eeb490ec9514930a48c *inst/doc/CoalescentModels.Rnw
-8da39d8199433c96cf9ee1c1fa9ff5e1 *inst/doc/CoalescentModels.pdf
-0d5a83d0784207bb0911e0352170ab6b *inst/doc/Running_coalescentMCMC.R
-e2e9e5b527e95804201ce3d789cc2eed *inst/doc/Running_coalescentMCMC.Rnw
-90a4f0bf20a91e8348586a180be3c104 *inst/doc/Running_coalescentMCMC.pdf
-5776e0661b44e5213b0689cc176aa3d9 *man/coalesceMCMC.Rd
+3c1e57130cbd9694fa95ddcefd4846be *inst/doc/CoalescentModels.Rnw
+95e7d0650c8a8983fd0a380d5e32f687 *inst/doc/CoalescentModels.pdf
+bd3f5e6b53673e4c4fef7ddc8f210fa6 *inst/doc/Running_coalescentMCMC.R
+bdb91a1646d610bed17b9724df2f4fb6 *inst/doc/Running_coalescentMCMC.Rnw
+70ce8a00421ca3363249921b3598e408 *inst/doc/Running_coalescentMCMC.pdf
+5067e244a087c5ede08c7b231f096460 *man/coalesceMCMC.Rd
 9aa980aea36886b18ac6a09b823765da *man/coalescentMCMC-internal.Rd
 f768c542f46f407ec6cce19410f749cc *man/dcoal.Rd
 90b1945d199985a0b7162999b6f140d0 *man/treeOperators.Rd
 2a6f9e9e044a78154d3cfda5936d6f48 *src/Makevars
 c3e22b2a119684f9a878926188ddd39e *src/coalescentMCMC.c
-38240cef46b41eeb490ec9514930a48c *vignettes/CoalescentModels.Rnw
-e2e9e5b527e95804201ce3d789cc2eed *vignettes/Running_coalescentMCMC.Rnw
-8bc7717e10b96e671d31408560cf5981 *vignettes/coalescentMCMC.bib
+3c1e57130cbd9694fa95ddcefd4846be *vignettes/CoalescentModels.Rnw
+bdb91a1646d610bed17b9724df2f4fb6 *vignettes/Running_coalescentMCMC.Rnw
+926f24b46ead272ee5d5f742d386b807 *vignettes/coalescentMCMC.bib

NAMESPACE

@@ -1,10 +1,13 @@
 useDynLib(coalescentMCMC, .registration = TRUE)
 
-export(.coalescentMCMCenv, coalescentMCMC, dcoal, dcoal.linear,
-       dcoal.step, dcoal.time, dcoal.time2, EdgeLengthJittering,
-       getMCMCtrees, NeighborhoodRearrangement, TipInterchange)
+export(.coalescentMCMCenv, cleanMCMCtrees, coalescentMCMC, dcoal,
+       dcoal.linear, dcoal.step, dcoal.time, dcoal.time2,
+       EdgeLengthJittering, getLastTree, getMCMCstats, getMCMCtrees,
+       NeighborhoodRearrangement, saveMCMCtrees, TipInterchange)
 
-importFrom(ape, as.phylo, branching.times, dist.dna, Ntip, reorder.phylo)
+importFrom(ape, as.phylo, branching.times, dist.dna, Ntip,
+           reorder.phylo, write.nexus, write.tree)
 importFrom(coda, mcmc)
-importFrom(phangorn, phyDat, pml)
+importFrom(phangorn, edQt, lli, phyDat, pml, pml.fit, pml.free, pml.init)
+importFrom(Matrix, Matrix)
 importFrom(stats, hclust, rbinom, reorder, runif)

NEWS

@@ -1,3 +1,24 @@
+		CHANGES IN coalescentMCMC VERSION 0.4
+
+
+NEW FEATURES
+
+    o Three new functions, getLastTree, saveMCMCtrees, and
+      cleanMCMCtrees, provide tools to manage lists of trees.
+
+    o Two new models are available in coalescentMCMC(): THETA varying
+      linearly with time, or in a stepwise way (with a single
+      breakpoint in time). All available models are described in a
+      vignette.
+
+    o It is now possible to interrupt coalescentMCMC(): the results at
+      this stage are output normally.
+
+    o getMCMCtrees() has a new option (chain) to extract directly a
+      specific list of trees bypassing the interactive menu.
+
+
+
 		CHANGES IN coalescentMCMC VERSION 0.3
 
 

R/coalescentMCMC.R

@@ -1,4 +1,4 @@
-## coalescentMCMC.R (2013-10-18)
+## coalescentMCMC.R (2013-12-04)
 
 ##   Run MCMC for Coalescent Trees
 
@@ -11,6 +11,39 @@ coalescentMCMC <-
     function(x, ntrees = 3000, burnin = 1000, frequency = 1,
              tree0 = NULL, model = NULL, quiet = FALSE)
 {
+    on.exit({
+        if (k < nOut) {
+            if (!k) stop("burn-in period not yet finished")
+            TREES <- TREES[seq_len(k)]
+            LL <- LL[seq_len(i)]
+            params <- params[seq_len(i), ]
+            warning(paste("MCMC interrupted after", i, "generations"))
+        }
+        ## compress the list of trees:
+        attr(TREES, "TipLabel") <- TREES[[1L]]$tip.label
+        lapply(TREES, function(x) x$tip.label <- NULL)
+        class(TREES) <- "multiPhylo"
+
+        suffix <- 1
+        list.trees <- .get.list.trees()
+        if (l <- length(list.trees))
+            suffix <- 1 + as.numeric(sub("TREES_", "", list.trees[l]))
+        suffix <- sprintf("%03d", suffix)
+        assign(paste("TREES", suffix, sep = "_"), TREES,
+               envir = .coalescentMCMCenv)
+
+        i <- i - 1L
+
+        MCMCstats <- get("MCMCstats", envir = .coalescentMCMCenv)
+        MCMCstats[[suffix]] <- c(k, burnin, frequency, i, j)
+        assign("MCMCstats", MCMCstats, envir = .coalescentMCMCenv)
+
+        LL <- cbind(LL, params)
+        colnames(LL) <- c("logLik", para.nms)
+        LL <- mcmc(LL, start = 1, end = i)
+        return(LL)
+    })
+
     if (is.null(tree0)) {
         d <- dist.dna(x, "JC69")
         tree0 <- as.phylo(hclust(d, "average"))
@@ -24,6 +57,21 @@ coalescentMCMC <-
 
     getlogLik <- function(phy, X) pml(phy, X)$logLik
 
+    if (packageVersion("phangorn") >= "1.99.5") {
+        INV <- Matrix(lli(X, tree0), sparse = TRUE)
+        ll.0 <- numeric(attr(X, "nr"))
+        bf <- rep(0.25, 4)
+        eig <- edQt()
+        getlogLik <- function(phy, X) {
+            phy <- reorder(phy, "postorder")
+            pml.init(X)
+            loglik <- pml.fit(phy, X, bf = bf, eig = eig,
+                              INV = INV, ll.0 = ll.0)
+            pml.free()
+            loglik
+        }
+    }
+
     TREES <- vector("list", nOut)
     LL <- numeric(nOut2)
     TREES[[1L]] <- tree0
@@ -36,26 +84,55 @@ coalescentMCMC <-
         ## quantities to calculate THETA:
         two2n <- 2:n
         K4theta <- length(two2n)
-        tmp <- choose(two2n, 2)
+        tmp <- two2n * (two2n - 1) # == 2 * choose(two2n, 2)
         getparams <- function(phy, bt) {
             x4theta <- rev(diff(c(0, sort(bt))))
             sum(x4theta * tmp)/K4theta
         }
         f.theta <- function(t, p) p
-    } else { # only "time" (ie, exponential model) for the moment
-        np <- 2L
-        para.nms <- c("theta0", "rho")
-        getparams <- function(phy, bt) { # 'bt' is not used but is needed to have the same arguments than above
-            out <- nlminb(c(0.02, 0),
-                          function(p) -dcoal.time(phy, p[1], p[2], log = TRUE))
-            out$par
-        }
-        f.theta <- function(t, p) p[1] * exp(p[2] * t)
+    } else {
+        switch(model, time = {
+            np <- 2L
+            para.nms <- c("theta0", "rho")
+            getparams <- function(phy, bt) { # 'bt' is not used but is needed to have the same arguments than above
+                halfdev <- function(p) {
+                    if (any(p <= 0) || any(is.nan(p))) return(1e100)
+                    -dcoal.time(phy, p[1], p[2], log = TRUE)
+                }
+                out <- nlminb(c(0.02, 0), halfdev)
+                out$par
+            }
+            f.theta <- function(t, p) p[1] * exp(p[2] * t)
+        }, step = {
+            np <- 3L
+            para.nms <- c("theta0", "theta1", "tau")
+            getparams <- function(phy, bt) {
+                halfdev <- function(p) {
+                    if (any(p <= 0) || any(is.nan(p))) return(1e100)
+                    -dcoal.step(phy, p[1], p[2], p[3], log = TRUE)
+                }
+                out <- nlminb(c(0.02, 0.02, bt[1]/2), halfdev)
+                out$par
+            }
+            f.theta <- function(t, p) ifelse(t <= p[3], p[1], p[2])
+        }, linear = {
+            np <- 3L
+            para.nms <- c("theta0", "thetaT", "TMRCA")
+            getparams <- function(phy, bt) {
+                halfdev <- function(p) {
+                    if (any(p <= 0) || any(is.nan(p))) return(1e100)
+                    -dcoal.linear(phy, p[1], p[2], p[3], log = TRUE)
+                }
+                out <- nlminb(c(0.02, 0.02, bt[1]), halfdev)
+                out$par
+            }
+            f.theta <- function(t, p) p[1] + t * (p[2] - p[1])/p[3]
+        })
     }
     params <- matrix(0, nOut2, np)
 
     i <- 2L
-    j <- 0L # number of accepted trees
+    j <- 0L # number of accepted moves
     k <- 0L # number of sampled trees
 
     if (!quiet) {
@@ -64,10 +141,9 @@ coalescentMCMC <-
         cat("  Burn-in period:", burnin, "\n")
         cat("  Sampling frequency:", frequency, "\n")
         cat("  Number of generations to run:", ntrees * frequency + burnin, "\n")
-        cat("Generation    Nb of accepted trees\n")
+        cat("Generation    Nb of accepted moves\n")
     }
 
-    ##
     bt0 <- branching.times(tree0)
     params[1L, ] <- para0 <- getparams(tree0, bt0)
 
@@ -84,8 +160,9 @@ coalescentMCMC <-
         tr.b <- NeighborhoodRearrangement(tree0, n, nodeMax, target, THETA, bt0)
         ## do TipInterchange() every 10 steps:
         ## tr.b <-
-        ##     if (!i %% 10) TipInterchange(tree0, n)
+        ##     if (! i %% 10) TipInterchange(tree0, n)
         ##     else NeighborhoodRearrangement(tree0, n, nodeMax, target, THETA, bt0)
+
         if (!(i %% frequency) && i > burnin) {
             k <- k + 1L
             TREES[[k]] <- tr.b
@@ -108,36 +185,70 @@ coalescentMCMC <-
             bt0 <- bt
         }
     }
-
-    #dim(LL) <- c(i - 1, 1)
-    LL <- cbind(LL, params)
-    colnames(LL) <- c("logLik", para.nms)
-    LL <- mcmc(LL, start = 1, end = i - 1)
-
-    ## compress the list of trees:
-    attr(TREES, "TipLabel") <- TREES[[1L]]$tip.label
-    for (i in seq_len(nOut)) TREES[[i]]$tip.label <- NULL
-    class(TREES) <- "multiPhylo"
-
-    j <- 1
-    list.trees <- ls(envir = .coalescentMCMCenv)
-    if (l <- length(list.trees))
-        j <- 1 + as.numeric(sub("TREES_", "", list.trees[l]))
-    assign(paste("TREES", j, sep = "_"), TREES, envir = .coalescentMCMCenv)
     if (!quiet) cat("\nDone.\n")
-    LL
 }
 
-getMCMCtrees <- function() {
-    list.trees <- ls(envir = .coalescentMCMCenv)
+.get.list.trees <- function()
+    ls(envir = .coalescentMCMCenv, pattern = "^TREES_")
+
+getMCMCtrees <- function(chain = NULL)
+{
+    list.trees <- .get.list.trees()
+    l <- length(list.trees)
+    if (is.null(chain)) {
+        if (!l) return(NULL)
+        if (l == 1)
+            return(get(list.trees, envir = .coalescentMCMCenv))
+        ## l > 1:
+        cat("Several lists of MCMC trees are stored:\n\n")
+        for (i in 1:l) cat(i, ":", list.trees[i], "\n")
+        cat("\nReturn which number? ")
+        chain <- as.numeric(readLines(n = 1))
+    } else {
+        if (!l) {
+            warning("no list of MCMC trees stored")
+            return(NULL)
+        }
+        if (l < chain) {
+            warning("no enough lists of MCMC trees stored")
+            return(NULL)
+        }
+    }
+    get(paste("TREES", sprintf("%03d", chain), sep = "_"),
+        envir = .coalescentMCMCenv)
+}
+
+saveMCMCtrees <- function(destdir = ".", format = "RDS", ...)
+{
+    format <- match.arg(toupper(format), c("RDS", "NEWICK", "NEXUS"))
+    switch(format, RDS = {
+        FUN <- saveRDS
+        suffix <- ".rds"
+    }, NEWICK = {
+        FUN <- write.tree
+        suffix <- ".tre"
+    }, NEXUS = {
+        FUN <- write.nexus
+        suffix <- ".nex"
+    })
+    list.trees <- .get.list.trees()
     l <- length(list.trees)
-    if (!l) return(NULL)
-    if (l == 1)
-        return(get(list.trees, envir = .coalescentMCMCenv))
-    ## l > 1:
-    cat("Several lists of MCMC trees are stored:\n\n")
-    for (i in 1:l) cat(i, ":", list.trees[i], "\n")
-    cat("\nReturn which number? ")
-    i <- as.numeric(readLines(n = 1))
-    get(paste("TREES", i, sep = "_"), envir = .coalescentMCMCenv)
+    if (!l) warning("no list of trees to save") else {
+        for (i in 1:l) {
+            f <- list.trees[i]
+            outfile <- paste(destdir, "/", f, suffix, sep = "")
+            FUN(get(f, envir = .coalescentMCMCenv), outfile, ...)
+        }
+    }
+}
+
+cleanMCMCtrees <- function()
+    rm(list = .get.list.trees(), envir = .coalescentMCMCenv)
+
+getLastTree <- function(X) X[[length(X)]]
+
+getMCMCstats <- function()
+{
+    cat("MCMC chain summaries (chains as columns):\n\n")
+    get("MCMCstats", envir = .coalescentMCMCenv)
 }