Updates for simulation example

NAMESPACE

@@ -14,10 +14,12 @@ export(calcNdBlnc)
 export(calcNdsBlnc)
 export(calcOvrlp)
 export(calcPhyDv)
+export(calcPrtFrPrp)
 export(checkTreeMan)
 export(checkTreeMen)
 export(fastCheckTreeMan)
 export(getAge)
+export(getCnnctdNds)
 export(getDcsd)
 export(getLvng)
 export(getNdAge)
@@ -45,6 +47,7 @@ export(getSpnAge)
 export(getSpnsAge)
 export(getSubtree)
 export(getTxnyms)
+export(getUnqNds)
 export(loadTreeMan)
 export(pinTips)
 export(pstMnp)

R/calc-methods.R

@@ -253,6 +253,7 @@ calcDstRF <- function(tree_1, tree_2, nrmlsd=FALSE,
 #' @details Faith's phylogenetic diversity is calculated as the sum of all connected
 #' branches for specified tips in a tree. It can be used to investigate how biodviersity
 #' as measured by the phylogeny changes. Parallelizable.
+#' The function uses \code{getCnntdNds()}.
 #' @param tree \code{TreeMan} object
 #' @param tids tip ids
 #' @param parallel logical, make parallel?
@@ -261,7 +262,7 @@ calcDstRF <- function(tree_1, tree_2, nrmlsd=FALSE,
 #' Faith, D. (1992). Conservation evaluation and phylogenetic diversity.
 #'  Biological Conservation, 61, 1-10.
 #' @seealso
-#' \code{\link{calcFrPrp}}, \code{\link{calcOvrlp}}
+#' \code{\link{calcFrPrp}}, \code{\link{calcOvrlp}}, \code{\link{getCnntdNds}}
 #' \url{https://github.com/DomBennett/treeman/wiki/calc-methods}
 #' @export
 #' @examples
@@ -270,10 +271,7 @@ calcDstRF <- function(tree_1, tree_2, nrmlsd=FALSE,
 #' calcPhyDv(tree, tree['tips'])
 calcPhyDv <- function(tree, tids,
                       parallel=FALSE, progress="none") {
-  prids <- c(unique(unlist(getNdsPrids(tree, tids))),
-             tids)
-  counts <- table(prids)
-  prids <- names(counts)[counts < length(tids)]
+  prids <- getCnnctdNds(tree, tids)
   spns <- getNdsSlt(tree, slt_nm="spn", ids=prids,
                     parallel=parallel, progress=progress)
   sum(spns)
@@ -292,7 +290,7 @@ calcPhyDv <- function(tree, tids,
 #' Isaac, N.J.B., Turvey, S.T., Collen, B., Waterman, C. and Baillie, J.E.M. (2007). 
 #'  Mammals on the EDGE: conservation priorities based on threat and phylogeny. PLoS ONE, 2, e296.
 #' @seealso
-#' \code{\link{calcPhyDv}}
+#' \code{\link{calcPhyDv}}, \code{\link{calcPrtFrPrp}},
 #' \url{https://github.com/DomBennett/treeman/wiki/calc-methods}
 #' @export
 #' @examples
@@ -303,10 +301,10 @@ calcFrPrp <- function(tree, tids, progress="none") {
   .calc <- function(i) {
     id <- tree@all[i]
     spn <- getNdSlt(tree, "spn", id)
-    kids <- getNdKids(tree, id)
-    if(length(kids) == 0) {
+    if(id %in% tree@tips) {
       spn_shres[i, id] <<- spn
     } else {
+      kids <- getNdKids(tree, id)
       spn_shre <- spn/length(kids)
       spn_shres[i, kids] <<- spn_shre
     }
@@ -318,6 +316,51 @@ calcFrPrp <- function(tree, tids, progress="none") {
   colSums(spn_shres[, tids])
 }
 
+#' @name calcPrtFrPrp
+#' @title Calculate evolutionary distinctness for part of tree
+#' @description Returns the evolutationary distinctness of ids using the fair proportion metric.
+#' @details Extension of \code{calcFrPrp()} but with ignore argument.
+#' Use \code{ignr} to ignore certain tips from calculation. For example, if any of tips
+#' are extinct you may wish to ignore these.
+#' @param tree \code{TreeMan} object
+#' @param tids tip IDs
+#' @param ignr tips to ignore in calculation
+#' @param progress name of the progress bar to use, see \code{\link{create_progress_bar}}
+#' @references
+#' Isaac, N.J.B., Turvey, S.T., Collen, B., Waterman, C. and Baillie, J.E.M. (2007). 
+#'  Mammals on the EDGE: conservation priorities based on threat and phylogeny. PLoS ONE, 2, e296.
+#' @seealso
+#' \code{\link{calcFrPrp}}
+#' \url{https://github.com/DomBennett/treeman/wiki/calc-methods}
+#' @export
+#' @examples
+#' library(treeman)
+#' tree <- randTree(10)
+#' calcPrtFrPrp(tree, c('t1','t3'), ignr='t2')
+calcPrtFrPrp <- function(tree, tids, ignr=NULL, progress="none") {
+  .calc <- function(i) {
+    id <- allnds[i]
+    spn <- getNdSlt(tree, "spn", id)
+    if(id %in% tips) {
+      spn_shres[i, id] <<- spn
+    } else {
+      kids <- getNdKids(tree, id)
+      kids <- kids[!kids %in% ignr]
+      if(length(kids) > 0) {
+        spn_shre <- spn/length(kids)
+        spn_shres[i, kids] <<- spn_shre
+      }
+    }
+  }
+  tips <- tree@tips[!tree@tips %in% ignr]
+  allnds <- tree@all[!tree@all %in% ignr]
+  spn_shres <- matrix(0, ncol=length(tips), nrow=length(allnds))
+  colnames(spn_shres) <- tips
+  plyr::m_ply(.data=data.frame(i=1:length(allnds)), .fun = .calc,
+              .progress=progress)
+  colSums(spn_shres[, tids])
+}
+
 #' @name calcDstMtrx
 #' @title Calculate the distance matrix
 #' @description Returns a distance matrix for specified ids of a tree.

R/get-spcl-methods.R

@@ -152,6 +152,49 @@ getOtgrp <- function(tree, ids) {
 
 # SPECIAL
 
+#' @name getUnqNds
+#' @title Get unique nodes represented by tips
+#' @description Return a list of IDs for any node that are represented by tip IDs given.
+#' @details Returns a vector.
+#' @param tree \code{TreeMan} object
+#' @param tids vector of tip IDs
+#' @seealso
+#' \code{\link{getCnnctdNds}}, \code{\link{calcFrPrp}},
+#' \code{\link{calcPhyDv}}
+#' @export
+#' @examples
+#' library(treeman)
+#' tree <- randTree(10)
+#' unqnds <- getUnqNds(tree, c('t1', 't2'))
+getUnqNds <- function(tree, tids) {
+  rmng <- tree@tips[!tree@tips %in% tids]
+  ignr <- c(unique(unlist(getNdsPrids(tree, tids))), tids)
+  rmng <- c(unique(unlist(getNdsPrids(tree, rmng))), rmng)
+  ignr[!ignr %in% rmng]
+}
+
+#' @name getCnnctdNds
+#' @title Get all nodes connected by given tips
+#' @description Return a vector of IDs of all nodes that are connected to tip IDs given.
+#' @details Returns a vector. This function is the basis for \code{calcPhyDv()}, it determines
+#' the unique set of nodes connected for a set of tips.
+#' @param tree \code{TreeMan} object
+#' @param tids vector of tip IDs
+#' @seealso
+#' \code{\link{getUnqNds}}, \code{\link{calcFrPrp}},
+#' \code{\link{calcPhyDv}}
+#' @export
+#' @examples
+#' library(treeman)
+#' tree <- randTree(10)
+#' cnntdnds <- getCnnctdNds(tree, c('t1', 't2'))
+getCnnctdNds <- function(tree, tids) {
+  prids <- c(unlist(getNdsPrids(tree, tids)),
+             tids)
+  counts <- table(prids)
+  names(counts)[counts < length(tids)]
+}
+
 #' @name getTxnyms
 #' @title Get IDs for nodes represented txnyms
 #' @description Return a list of IDs for any node that contains the given txnyms.

R/update-methods.R

@@ -62,6 +62,7 @@ updateSlts <- function(tree) {
 #' The matrix is generated with bigmemory's `as.big.matrix()`.
 #' @param tree \code{TreeMan} object
 #' @param shared T/F, should the bigmatrix be shared? See bigmemory documentation.
+#' @param ... \code{as.big.matrix()} additional arguments
 #' @seealso
 #' \code{\link{updateSlts}}, \code{\link{rmNdmtrx}},
 #' \link{https://cran.r-project.org/web/packages/bigmemory/index.html}

other/2_simulation.R

@@ -7,52 +7,51 @@ library(treeman)
 # PARAMETER
 # balanced tree to start
 tree_string <- "((A:1.0,B:1.0):1.0,(C:1.0,D:1.0):1.0);"
+tree_age <- 2
 tree <- readTree(text=tree_string)
 iterations <- 200
 burnin <- 10
 d <- 1
-b <- 2  # birth for burnin
+b <- 3  # birth for burnin
 b_true <- 1  # birth after burnin
-ext <- tree["tips"]
+extnt <- tree["tips"]
+extnct <- NULL
 
 # LOOP
 cat('Simulating ....\n')
 for(i in 1:iterations) {
-  if(length(ext) < 3) {
-    stop('Too few tips remaining!')
+  if(length(extnt) < 3) {
+    stop('Too few extant tips remaining!\nRun again or increase birth before burnin.')
   }
   if(i > burnin) {
     b <- b_true
   }
   cat('.... i=[', i, ']\n', sep='')
   # calculate fair proportion
-  fps <- calcFrPrp(tree, ext)
+  fps <- calcPrtFrPrp(tree, tids=extnt, ignr=extnct)
   # add/remove based on b and d
   to_add <- sample(c(TRUE, FALSE), size=1, prob=c(b,d))
   if(to_add) {
-    sid <- sample(ext, prob=1/fps, size=1)  # sister ID of the new tip
+    sid <- sample(extnt, prob=1/fps, size=1)  # sister ID of the new tip
     tid <- paste0('t', i)  # new tip ID
-    tree <- treeman::addTip(tree, tid=tid, sid=sid, start=0, end=0)
+    tree <- addTip(tree, tid=tid, sid=sid, strt_age=0, end_age=0,
+                   tree_age=tree_age)
+    extnt <- c(extnt, tid)
   } else {
-    tid <- sample(ext, prob=1/fps, size=1)
-    tree <- rmTip(tree, tid=tid)
+    tid <- sample(extnt, prob=1/fps, size=1)
+    extnct <- c(extnct, tid)
+    extnt <- extnt[extnt != tid]
   }
   # grow tree
-  ext <- tree['tips']
-  spns <- getNdsSlt(tree, slt_nm="spn", ids=ext)
-  tree <- setNdsSpn(tree, ids=ext, vals=spns+1)
+  spns <- getNdsSlt(tree, slt_nm="spn", ids=extnt)
+  tree <- setNdsSpn(tree, ids=extnt, vals=spns+1)
+  tree_age <- tree_age + 1
 }
 cat('Done.\n')
 
 # VIZ
-library(MoreTreeTools)
-tree_phylo <- as(tree, 'phylo')
-# plot simulated tree with edges coloured by proximate diversity
-tree_phylo$edge.label <- paste0 ('edge_', 1:nrow(tree_phylo$edge))
-# intervals are used to calculate the colour of the branch
-# diversity is the number of descedents of a branch with an interval
-ed <- calcEdgeDiversity(tree_phylo, n.intervals=10)
-ed$col <- (log(ed$count) - mean(log(ed$count))) /
-  sd(log(ed$count))
-p <- chromatophylo(tree_phylo, edge.cols=ed, legend.title='Diversity')
-print(p)
+# plot the fossil record tree, and the reconstructed phylogeny
+extnt_tree <- rmTips(tree, tids=getDcsd(tree))
+par(mfrow=c(1,2))
+plot(as(tree, 'phylo'), show.tip.label=FALSE, main='Fossil record')
+plot(as(extnt_tree, 'phylo'), main='Reconstructed')