Skip to content

Commit

Permalink
Simplified code, progress on ropensci/software-review#209
Browse files Browse the repository at this point in the history
  • Loading branch information
richelbilderbeek committed Apr 13, 2018
1 parent 87a7495 commit 29fe259
Show file tree
Hide file tree
Showing 2 changed files with 1 addition and 69 deletions.
2 changes: 0 additions & 2 deletions NAMESPACE
Original file line number Diff line number Diff line change
@@ -1,6 +1,5 @@
# Generated by roxygen2: do not edit by hand

export(beast2out.read.trees)
export(calc_act)
export(calc_act_r)
export(calc_ess)
Expand All @@ -20,7 +19,6 @@ export(parse_beast_output_files)
export(parse_beast_posterior)
export(parse_beast_state_operators)
export(parse_beast_trees)
export(parse_beast_trees_impl_2)
export(remove_burn_in)
export(remove_burn_ins)
importFrom(Rcpp,sourceCpp)
Expand Down
68 changes: 1 addition & 67 deletions R/parse_beast_trees.R
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,7 @@ parse_beast_trees <- function(
}

trees <- tryCatch({
parse_beast_trees_impl_2(filename)
ape::read.nexus(filename)
},
error = function(cond) {
stop("invalid file")
Expand All @@ -32,69 +32,3 @@ parse_beast_trees <- function(
class(trees) <- "multiPhylo"
trees
}

#' Parses a BEAST2 .trees output file
#' @param filename name of the BEAST2 .trees output file
#' @return the phylogenies in the posterior
#' @export
#' @examples
#' trees_filename <- get_tracerer_path("beast2_example_output.trees")
#' posterior <- parse_beast_trees(trees_filename)
#' testit::assert(is_trees_posterior(posterior))
#' @author Richel J.C. Bilderbeek
#' @noRd
parse_beast_trees_impl_2 <- function(
filename
) {
ape::read.nexus(filename)
}


#' Extract a list of phylogenies from a BEAST2 posterior file
#' @param file name of the BEAST2 posterior filename, usually ends with '.trees'
#' @param opt.rescale.edge.length Rescaling factor of the phylogenies,
#' a value of 1.0 denotes that the original edge lengths are maintained
#' @param opt.burnin how many phylogenies to discard, an value of will keep all trees
#' @return a list of phylogenies of type 'phylo'
#' @examples
#' trees_file <- "vignettes/example.trees"
#' testit::assert(file.exists(trees_file))
#' posterior <- beast2out.read.trees(trees_file)
#' testit::assert(length(posterior) == 10)
#' testit::assert(class(posterior[[1]]) == "phylo")
#' @export
#' @author Oliver Ratmann
#' @noRd
beast2out.read.trees <- function(
file,
opt.rescale.edge.length= 1.,
opt.burnin=0
)
{
tmp <- readLines(file, n=2e3, warn = FALSE)
tmp <- which( grepl('#NEXUS', tmp) )
if(length(tmp)>1)
{
cat(paste('\nFound #NEXUS headers, n=',length(tmp),'.\nDiscard all lines before last entry on line', tail(tmp,1)))
cmd <- paste('sed -i".bak" 1,',tail(tmp,1)-1,'d ', file, sep='')
system(cmd)
cmd <- paste('sed -i".bak2" 1s/\\;// ', file, sep='')
system(cmd)
cmd <- list.files(paste(rev(rev(strsplit(file, '/')[[1]])[-1]),collapse='/'), pattern='*bak*', full.names=TRUE)
cat(paste('\nrm files\n', paste(cmd, collapse='\n')))
file.remove(cmd)
}
mph <- ape::read.nexus(file)
# remove burn in
tmp <- regexpr('[0-9]+',names(mph))
if(any(tmp<0)) stop('unexpected nexus file without STATE iteration numbers')
mph.it <- as.numeric( regmatches( names(mph), tmp) )
mph <- lapply( which( mph.it>opt.burnin), function(j) mph[[j]] )
mph.it <- mph.it[ mph.it > opt.burnin ]
names(mph) <- paste('STATE_',mph.it,sep='')
# rescale edge lengths
if(opt.rescale.edge.length!=1.)
for(j in seq_along(mph))
mph[[j]]$edge.length <- mph[[j]]$edge.length * opt.rescale.edge.length
mph
}

0 comments on commit 29fe259

Please sign in to comment.