R/conversions.R

Defines functions fossils.to.pyrate fossils.to.paleotree.record paleotree.record.to.fossils beast.fbd.format prune.fossils reconstructed.tree.fossils.objects sampled.tree.from.combined

Documented in beast.fbd.format fossils.to.paleotree.record fossils.to.pyrate paleotree.record.to.fossils prune.fossils reconstructed.tree.fossils.objects sampled.tree.from.combined

#' Removes all unsampled lineages from a combined tree.
#' Extinct tips are only sampled if they are fossils. With default settings all extant tips are sampled.
#'
#' @param tree Combined tree with fossils.
#' @param rho Sampling probability of extant tips. Default 1, will be disregarded if sampled_tips is not null.
#' @param sampled_tips List of tip labels corresponding to sampled extant tips.
#' @return Sampled tree with fossils.
#' @examples
#' # simulate tree
#' t = ape::rtree(6)
#'
#' # simulate fossils
#' f = sim.fossils.poisson(rate = 2, tree = t)
#'
#' # transform format
#' t2 = SAtree.from.fossils(t,f)$tree
#'
#' # transform to sampled tree
#' t3 = sampled.tree.from.combined(t2)
#' plot(t3)
#' @export
sampled.tree.from.combined = function(tree, rho = 1, sampled_tips = NULL) {
  if(!("SAtree" %in% class(tree)) ){
    if("phylo" %in% class(tree)) tree = SAtree(tree)
    else stop(paste('object "',class(tree),'" is not of class "SAtree"',sep=""))
  }
  if(!tree$complete && rho == 1 && is.null(sampled_tips)) stop("Tree is already sampled")

  remove_tips = c()

  depths = ape::node.depth.edgelength(tree)
  times = max(depths) - depths

  for(i in 1:length(tree$tip.label)) {
    if(times[i] < 1e-5) { #extant tip
      if((!is.null(sampled_tips) && !tree$tip.label[i] %in% sampled_tips) || #tip not sampled from sampled_tips
         (is.null(sampled_tips) && runif(1) > rho)) { #tip not sampled from rho
        remove_tips = c(remove_tips, i)
      }
    }
    else if(tree$complete) { #extinct tip
      edge = which(tree$edge[,2]==i)
      if(tree$edge.length[edge] > 1e-5) { #not on zero-length edge = not a fossil
        remove_tips = c(remove_tips, i)
      }
    }
  }

  tree = ape::drop.tip(tree, remove_tips)
  tree$complete = FALSE
  tree
}

#' Returns tree and fossil objects that you can use to plot the reconstructed tree.
#'
#' Note that for datasets containing extinct only samples (& rho = 0) the ages output are scaled so that the youngest sample = 0.
#'
#' @param tree Tree object.
#' @param fossils Fossils object.
#' @param rho Extant species sampling probability. Default = 1, will be disregarded if fossils object already contains extant samples.
#'
#' @return A list containing the tree and fossil objects.
#'
#' @examples
#' # simulate tree
#' birth = 0.1
#' death = 0.05
#' tips = 10
#' t = TreeSim::sim.bd.taxa(tips, 1, birth, death)[[1]]
#'
#' # simulate fossils
#' f = sim.fossils.poisson(rate = 0.3, tree = t)
#'
#' # simulate extant samples
#' f = sim.extant.samples(f, tree = t, rho = 0.5)
#'
#' # plot the complete tree
#' plot(f,t)
#'
#' # generate tree & fossil objects corresponding to the reconstructed tree
#' out = reconstructed.tree.fossils.objects(f, t)
#' f.reconst = out$fossils
#' t.reconst = out$tree
#'
#' # plot the reconstructed tree
#' plot(f.reconst, t.reconst)
#'
#' @export
# function to generate tree and corresponding fossil object for the reconstructed tree
reconstructed.tree.fossils.objects = function(fossils, tree, rho = 1){

  if(!is.null(tree) && !"phylo" %in% class(tree))
    stop("tree must be an object of class \"phylo\"")

  if(!is.null(fossils) && !"fossils" %in% class(fossils))
    stop("fossils must be an object of class \"fossils\"")

  if(!(rho >= 0 && rho <= 1))
    stop("rho must be a probability between 0 and 1")

  tol = min((min(tree$edge.length)/100),1e-8)
  samp_tips = NULL

  if(any( abs(fossils$hmax) < tol )){
    # identify extant sampled tips
    samp_tips = tree$tip.label[fossils$edge[which(abs(fossils$hmax) < tol)]]
    # eliminate them from the dataframe (required for SAtree.from.fossils)
    fossils = fossils[-c(which(abs(fossils$hmax) < tol)),]
  }

  # create SAtree object
  sa.tree = SAtree.from.fossils(tree, fossils)$tree

  # match samp_tips and sa.tree tip labels
  if(!is.null(samp_tips)){
    for(i in 1:length(samp_tips)){
      samp_tips[i] = sa.tree$tip.label [ grep(paste0(samp_tips[i],"_"), sa.tree$tip.label)[1] ]
    }
  }

  # removes all unsampled lineages
  samp.tree = sampled.tree.from.combined(sa.tree, rho = rho, sampled_tips = samp_tips)

  # identify sampled ancestors and drop from the sampled tree
  sa = c()
  for(i in 1:length(samp.tree$tip.label)){
    blength = samp.tree$edge.length[which(samp.tree$edge[,2]==i)]
    if(abs(blength) < tol) sa = c(sa, samp.tree$tip.label[i])
  }
  no.sa.tree = ape::drop.tip(samp.tree, sa)

  # create new fossils object based on the reconstructed tree
  # & deal with sampled ancestors
  if(length(fossils$sp) > 0){
    f.new = data.frame()
    nages = n.ages(samp.tree)
    for(i in sa){
      anc = ancestor(which(samp.tree$tip.label==i), samp.tree)

      s1 = ape::extract.clade(samp.tree,anc)$tip.label
      s2 = no.sa.tree$tip.label[which(no.sa.tree$tip.label %in% s1)]

      # assign fossils to nearest descendent node in the tree
      if(length(s2) > 1){
        j = ape::getMRCA(no.sa.tree,s2)
      } else { # i is SA on terminal branch
        j = which(no.sa.tree$tip.label==s2)
      }
      h = nages[[which(samp.tree$tip.label==i)]]
      f.new = rbind(f.new, data.frame(sp = j,
                                      edge = j,
                                      hmin = h,
                                      hmax = h))
    }
    f.new = fossils(f.new)
  } else {
    f.new = fossils()
  }

  # deal with non SA samples
  f.new = FossilSim::sim.tip.samples(f.new, no.sa.tree)

  # add a root edge for any stem fossils
  if( any(f.new$edge == (length(no.sa.tree$tip.label)+1) )){
    h = f.new$hmin[which(f.new$edge == (length(no.sa.tree$tip.label)+1))]
    no.sa.tree$root.edge = max(h) - max(n.ages(no.sa.tree))
  } else { no.sa.tree$root.edge = NULL }

  return(list(tree = no.sa.tree, fossils = f.new))

}

#' Removes all intermediate fossils from a combined tree and labels the first and last fossils of each lineage.
#' Can be used with sampled or complete trees. If only one fossil is present for a particular species it is labelled as first.
#'
#' @param tree Combined tree with fossils.
#' @return Tree with pruned fossils.
#' @examples
#' # simulate tree
#' t = ape::rtree(6)
#'
#' # simulate fossils
#' f = sim.fossils.poisson(rate = 2, tree = t)
#'
#' # transform format
#' t2 = SAtree.from.fossils(t,f)$tree
#'
#' # prune fossils
#' t4 = prune.fossils(t2)
#'
#' # or transform to sampled tree first
#' t3 = sampled.tree.from.combined(t2)
#' t4 = prune.fossils(t3)
#' plot(t4)
#' @export
prune.fossils = function(tree) {
  if(!("SAtree" %in% class(tree)) ){
    if("phylo" %in% class(tree)) tree = SAtree(tree)
    else stop(paste('object "',class(tree),'" is not of class "SAtree"',sep=""))
  }

  remove_tips = c()

  split_names = cbind(sub("_[^_]*$","",tree$tip.label),sub("^.+_","",tree$tip.label))
  for(name in unique(split_names[,1])) {
    idx = which(split_names[,1] == name)
    mx = max(split_names[idx,2])
    for(id in idx) {
      if(split_names[id,2] == 1) tree$tip.label[id] = paste0(name,"_first") # 1 corresponds to oldest sample in that lineage
      else if(mx >1 && split_names[id,2] == mx) tree$tip.label[id] = paste0(name,"_last") # earliest sample
      else remove_tips = c(remove_tips, id) # intermediate sample, to remove
    }
  }

  tree = ape::drop.tip(tree, remove_tips)
  tree
}

#' Transforms a tree and fossils into a sampled tree in beast-usable format and writes it in Newick format.
#' Designed to work with FBD.
#'
#' @param tree Complete tree.
#' @param fossils fossils dataframe.
#' @param rho Sampling probability of extant tips. Default 1, will be disregarded if sampled_tips is not null.
#' @param sampled_tips List of tip labels corresponding to sampled extant tips.
#' @param ... Additional parameters will be passed to ape::write.tree
#' @return Output of write.tree.
#' @examples
#' # simulate tree
#' t = ape::rtree(6)
#'
#' # simulate fossils
#' f = sim.fossils.poisson(rate = 2, tree = t)
#'
#' # output for BEAST
#' beast.fbd.format(t, f) # output on the console
#' \dontrun{
#' beast.fbd.format(t, f, file="example.tre") # output in file
#' }
#' @export
beast.fbd.format = function(tree, fossils, rho = 1, sampled_tips = NULL, ...) {
  proc_tree = prune.fossils(sampled.tree.from.combined(SAtree.from.fossils(tree,fossils)$tree, rho = rho, sampled_tips = sampled_tips))
  ape::write.tree(proc_tree, ...)
}

#' Transforms a fossilRecordSimulation object from package paleotree to a tree and taxonomy and fossils objects.
#'
#' The returned tree is in paleotree format, with zero-length edges leading to tips at bifurcation and anagenetic events.
#' Fossils and taxonomy are only specified on non-zero-length edges.
#' The label assigned to the parent of the origin or root will be zero.
#'
#' @param record fossilRecordSimulation object.
#' @param alphanumeric If TRUE function will return alphanumeric species labels (i.e. species labels contain the "t" prefix) (default). If FALSE function will return numeric only species labels.
#' @return A list containing the converted tree, taxonomy and fossils
#' @examples
#' if (requireNamespace("paleotree", quietly = TRUE)) {
#' # simulate record
#' record = paleotree::simFossilRecord(p=0.1, q=0.1,r=0.1, nruns=1, nTotalTaxa=c(30,40),
#'     nExtant=0, nSamp = c(5,25))
#'
#' # transform format
#' l_tf = paleotree.record.to.fossils(record)
#' l_tf$tree
#' l_tf$taxonomy
#' l_tf$fossils
#' }
#'
#' @export
#' @seealso \code{\link{taxonomy}}, \code{\link{fossils}}, \code{\link{fossils.to.paleotree.record}}
# NB: modes not branch specific
# NB: cryptic speciation: parent id is true parent, ie can be = cryptic id
paleotree.record.to.fossils = function(record, alphanumeric = TRUE) {
  # check that paleotree is installed - should be but you never know
  if (!requireNamespace("paleotree", quietly = TRUE)) {
    stop("Paleotree needed for this function to work. Please install it.", call. = FALSE)
  }

  tree = paleotree::taxa2phylo(paleotree::fossilRecord2fossilTaxa(record))
  # recording node labels to keep track after changing the phylogeny
  tree$node.label = (length(tree$tip.label)+1):(length(tree$tip.label)+tree$Nnode)
  fossildf = fossils()
  taxonomy = data.frame(edge = numeric(), sp = numeric(), start = numeric(), end = numeric(), mode = character(),
                        cryptic = numeric(), cryptic.id = numeric(), parent = numeric(), stringsAsFactors = F)
  ages = n.ages(tree)

  root_time = 0

  # this is a hack so the mode checking will work
  # basically paleotree outputs bifurcation nodes as 2 zero-length edges with the tip randomly assigned
  # here I'm reordering so the 2 edges always have the same ancestor
  for(e in 1:length(tree$edge.length)) {
    if(tree$edge.length[e] == 0 && tree$edge[e,2] <= length(record)) {
      parent = which(tree$edge[,2] == tree$edge[e,1])
      if(length(parent) > 0 && tree$edge.length[parent] == 0) {
        other_edge = which(tree$edge[,1] == tree$edge[parent,1])
        other_edge = other_edge[which(other_edge != parent)]
        tip = tree$edge[e,2]
        tree$edge[e,2] = tree$edge[other_edge,2]
        tree$edge.length[e] = tree$edge.length[other_edge]
        tree$edge[other_edge,2] = tip
        tree$edge.length[other_edge] = 0
      }
    }
  }

  for(i in 1:length(record)) {
    if(record[[i]]$taxa.data['orig.time'] > root_time) root_time = record[[i]]$taxa.data[['orig.time']]

    # recording positions of species on phylogeny
    node_idx = which(tree$tip.label == names(record)[i])
    sampled_nodes = c()
    if(tree$edge.length[which(tree$edge[,2] == node_idx)] == 0) {
      node_idx = tree$edge[which(tree$edge[,2] == node_idx),1]
    }
    all_nodes = node_idx
    age = ages[node_idx] + tree$edge.length[which(tree$edge[,2] == node_idx)] #age of the parent
    ends = record[[i]]$taxa.data[['ext.time']]
    starts = min(age, record[[i]]$taxa.data[["orig.time"]])

    for(t in sort(record[[i]]$sampling.times)) {
      while(node_idx != length(record) + 1 && t > age) {
        ends = c(ends, age)
        node_idx = tree$edge[which(tree$edge[,2] == node_idx),1]
        age = age + tree$edge.length[which(tree$edge[,2] == node_idx)]
        all_nodes = c(all_nodes, node_idx)
        starts = c(starts, min(age, record[[i]]$taxa.data[["orig.time"]]))
      }
      sampled_nodes = c(sampled_nodes, node_idx)
    }

    while(node_idx != length(record) + 1 && record[[i]]$taxa.data[["orig.time"]] > age) {
      prev = tree$edge[which(tree$edge[,2] == node_idx),1]
      if(prev %in% taxonomy$edge) break()
      lgth = tree$edge.length[which(tree$edge[,2] == prev)]
      if(length(lgth) > 0 && lgth == 0) break()

      ends = c(ends, age)
      node_idx = prev
      age = age + lgth
      all_nodes = c(all_nodes, node_idx)
      starts = c(starts, min(age, record[[i]]$taxa.data[["orig.time"]]))
    }

    if(is.na(record[[i]]$taxa.data[["ancestor.id"]])) { #root species
      parent = "t0"
      mode = 'r'
    } else {
      parent = names(record)[record[[i]]$taxa.data[['ancestor.id']]]
      origin = tree$edge[which(tree$edge[,2] == node_idx),1]
      orig_edge = which(tree$edge[,2] == origin)
      if(length(orig_edge) == 0 || tree$edge.length[orig_edge] > 0) {
        desc_anc = which(tree$edge[,1] == origin)
        if(any(tree$edge.length[desc_anc] == 0)) mode = 'a' # anagenetic event
        else mode = 'b' # budding
      } else {
        mode = 's' # bifurcation event
      }
    }

    taxonomy = rbind(taxonomy, data.frame(sp = names(record)[i], edge = all_nodes, parent = parent,
                                          start = starts, end = ends,
                                          mode = mode, cryptic = !(record[[i]]$taxa.data[['taxon.id']] == record[[i]]$taxa.data[['looks.like']]),
                                          cryptic.id = names(record)[record[[i]]$taxa.data[['looks.like']]],
                                          stringsAsFactors = F))
    if(length(record[[i]]$sampling.times)>0)
      fossildf = rbind(fossildf, data.frame(sp = names(record)[i], edge = sampled_nodes,
                                            hmin = sort(record[[i]]$sampling.times), hmax = sort(record[[i]]$sampling.times),
                                            stringsAsFactors = F))
  }

  row.names(taxonomy) = NULL
  row.names(fossildf) = NULL
  fossildf = as.fossils(fossildf, TRUE)
  taxonomy = as.taxonomy(taxonomy)

  if(!alphanumeric){
    fossildf$sp = gsub("t", "", fossildf$sp)
    taxonomy$sp = gsub("t", "", taxonomy$sp)
    taxonomy$cryptic.id = gsub("t", "", taxonomy$cryptic.id)
    taxonomy$parent = gsub("t", "", taxonomy$parent)
  }

  tree$root.edge = root_time - tree$root.time
  tree$origin.time = root_time

  #removing extant samples (if present) from fossils
  ext = which(fossildf$hmax < 1e-8)
  if(length(ext > 0)) fossildf = fossildf[-ext,]

  return(list(tree = tree, fossils = fossildf, taxonomy = taxonomy))
}

#' Transforms a fossils dataframe and either taxonomy or tree into a fossilRecordSimulation object from package paleotree.
#'
#' @param fossils fossils object
#' @param tree phylo object containing the tree. If provided and taxonomy = NULL, all speciation is assumed symmetric
#' @param taxonomy taxonomy object. If both tree and taxonomy are provided, only taxonomy will be used.
#' @return The converted paleotree record
#' @examples
#' # simulate tree
#' t = ape::rtree(6)
#' # simulate fossils using taxonomy
#' s = sim.taxonomy(t, 0.5, 1, 0.5)
#' f = sim.fossils.poisson(2, taxonomy = s)
#' # transform format
#' record = fossils.to.paleotree.record(f, taxonomy = s)
#' @seealso \code{\link{taxonomy}}, \code{\link{fossils}}, \code{\link{paleotree.record.to.fossils}}
#' @export
fossils.to.paleotree.record = function(fossils, tree = NULL, taxonomy = NULL) {
  if(is.null(taxonomy) && is.null(tree)) stop("Either tree or taxonomy needs to be provided")
  fossils = sim.extant.samples(fossils, tree = tree, taxonomy = taxonomy)

  rec_names = c("taxon.id","ancestor.id","orig.time","ext.time", "still.alive","looks.like")

  if(length(fossils$sp) > 0 & !any(grepl("t", fossils$sp))){
    fossils$sp = paste0("t", fossils$sp)
  }

  if(!is.null(taxonomy)) {
    # then record based purely on taxonomy

    # add "t" prefix if missing
    if(length(taxonomy$sp) > 0 & !any(grepl("t", taxonomy$sp))){
      taxonomy$sp = paste0("t", taxonomy$sp)
      taxonomy$parent = paste0("t", taxonomy$parent)
      taxonomy$cryptic.id = paste0("t", taxonomy$cryptic.id)
    }

    # order by species id
    taxonomy = taxonomy[order(as.numeric(gsub("t", "", taxonomy$sp))),]

    species = species.record.from.taxonomy(taxonomy)
    record = vector("list", length = length(species$sp))
    names(record) = species$sp
    for(i in 1:length(species$sp)) {
      if(is.na(species$parent[i]) | species$parent[i] == 0 | species$parent[i] == "t0") anc = NA else anc = which(names(record) == species$parent[i])
      record[[i]] = list(taxa.data = c(i, anc, species$species.start[i], species$species.end[i], (species$species.end[i] < 1e-5), which(names(record) == species$cryptic.id[i])),
                         sampling.times = numeric())
      names(record[[i]]$taxa.data) = rec_names

      f = which(fossils$sp == species$sp[i])
      for(fid in f) {
        record[[i]]$sampling.times = c(record[[i]]$sampling.times, (fossils$hmin[fid]+fossils$hmax[fid])/2)
      }
      record[[i]]$sampling.times = sort(record[[i]]$sampling.times, decreasing = T)
    }
  } else {
    #auxiliary function for handling one edge
    .convert_one_species = function(current_node, ancestor, record) {
      if(current_node <= length(tree$tip.label)) current_species = as.character(tree$tip.label[current_node])
      else if(!is.null(tree$node.label)) current_species = as.character(tree$node.label[current_node - length(tree$tip.label)])
      else current_species = paste0("t",current_node)

      above_node = tree$edge[which(tree$edge[,2] == current_node),1]
      record[[current_species]] = list(taxa.data = c(length(record) +1, ancestor, node.ages[above_node], node.ages[current_node], (node.ages[current_node] < 1e-3), length(record) +1),
                                       sampling.times = numeric())
      names(record[[current_species]]$taxa.data) = rec_names

      f = which(fossils$edge == current_node) #no taxonomy, all fossils on this branch belong to one species
      for(fid in f) {
        record[[current_species]]$sampling.times = c(record[[current_species]]$sampling.times, (fossils$hmin[fid]+fossils$hmax[fid])/2)
      }

      desc = tree$edge[which(tree$edge[,1] == current_node),2]
      for(d in desc) {
        record = .convert_one_species(d, which(names(record) == current_species), record)
      }
      record
    }

    # no taxonomy => parse tree assuming symmetric speciation everywhere
    node.ages = ape::node.depth.edgelength(tree)
    node.ages = max(node.ages) - node.ages

    current_node = length(tree$tip.label) + 1
    if(is.null(tree$root.edge)) {
      # no root edge => assuming that the tree started with 2 new species created at the root time
      if(!is.null(tree$origin.time)) node.ages = node.ages + (tree$origin.time - max(node.ages))
      record = list()
      desc = tree$edge[which(tree$edge[,1] == current_node),2]
      for(d in desc) {
        record = .convert_one_species(d, NA, record)
      }
    } else {
      # root edge present => assuming that the tree started with one species which existed only on the root edge
      newn = length(tree$tip.label) + tree$Nnode + 1
      tree$edge = rbind(c(newn, current_node), tree$edge)
      tree$edge.length = c(tree$root.edge, tree$edge.length)
      node.ages = c(node.ages, max(node.ages) + tree$root.edge)
      if(!is.null(tree$origin.time)) node.ages = node.ages + (tree$origin.time - max(node.ages))
      record = .convert_one_species(current_node, NA, list())
    }
  }

  class(record) = "fossilRecordSimulation"
  record
}

#' Generate output in the format used by the program PyRate
#'
#' @param fossils Fossils object.
#' @param python If TRUE the function outputs the data in the python format used by PyRate (default).
#' If FALSE the function outputs a tab-delimited table used by tools associated with PyRate.
#' @param traits Vector of trait values equal to the number of unique species in the fossils dataframe.
#' The order should correspond to the order in which they appear in \code{unique(fossils$sp)}.
#' @param cutoff Exclude occurrences with age uncertainty greater than this value i.e. \code{hmax - hmin > cutoff}.
#' @param random If TRUE use a random number from within the interval U(hmin, hmax) for specimen ages,
#' otherwise use the midpoint of this interval (default). Applicable only when \code{python = TRUE} and for specimens with \code{hmin != hmax}.
#' @param min Value used to represent the minimum possible interval age of extinct specimens with \code{hmin = 0}. By default \code{min = NULL} and the function will use the sampling times in the fossils dataframe.
#' @param exclude.extant.singletons If TRUE exclude species that have extant samples only (default = TRUE).
#' @param file Output file name.
#' @param use.sp.names If TRUE use the value in fossils$sp as the complete taxon name, otherwise the function adds the prefix "taxa" (default = FALSE).
#'
#' @examples
#'
#' set.seed(123)
#'
#' # simulate tree
#' t = ape::rtree(6)
#'
#' # assign a max age based on tree height
#' max.age = tree.max(t)
#'
#' # define a set of non-uniform length intervals
#' times = c(0, sort(runif(3, min = 0, max = max.age)), max.age)
#' rates = c(1,2,3,4)
#'
#' # simulate fossils reflect age uncertainty
#' f = sim.fossils.intervals(tree = t, interval.ages = times, rates = rates,
#'     use.exact.times = FALSE)
#'
#' # simulate extant samples
#' rho = 1
#' f = sim.extant.samples(f, t, rho = 1)
#'
#' plot(f, t)
#'
#' # generate input files for pyrate
#' fossils.to.pyrate(f)
#' fossils.to.pyrate(f, python = FALSE)
#'
#' # add trait values
#' traits = runif(length(unique(f$sp)))
#' fossils.to.pyrate(f, traits = traits)
#'
#' @export
fossils.to.pyrate = function(fossils, python = TRUE, traits = NULL, cutoff = NULL, random = FALSE, min = NULL,
                             exclude.extant.singletons = TRUE, file = "", use.sp.names = FALSE){

  if(length(fossils$sp) < 1) stop("Number of specimens must be > 0")

  if(!is.null(traits) && length(unique(fossils$sp)) != length(traits))
    stop("Number of trait values must equal the number of species")

  if(!is.null(traits)) traits = data.frame(sp = unique(fossils$sp), trait = round(traits, digits = 6))

  if(!is.null(cutoff) && !(cutoff > 0)) stop("Cutoff must be > 0")

  # tolerance for extant specimens age bounds = 0
  tol = 1e-8

  if(!is.null(min)){
    fossils$hmin[which(fossils$hmin != fossils$hmax && fossils$hmin < tol)] = min
    if(any(fossils$hmin > fossils$hmax)) stop("min value has generated intervals with hmax < hmin")
  }

  # deal with specimen age options
  if(!is.null(cutoff) && cutoff > 0){
    cat(paste0("Excluding ", length(which(fossils$hmax - fossils$hmin > cutoff)), " occurrences based on cutoff...\n"))
    fossils = fossils[which(fossils$hmax - fossils$hmin < cutoff),]
    if(length(fossils$sp) < 1) stop("Number of specimens after cutoff must be > 0")
    else if(exclude.extant.singletons && length(fossils[-which(fossils$hmin == fossils$hmax & fossils$hmin < tol),]$sp) < 1) stop("Number of specimens after cutoff must be > 0")
  }

  if(python){
    # use random ages
    if(random && any(fossils$hmax != fossils$hmin)){
      fossils$hmin = round(unlist(lapply(1:dim(fossils)[1], function(x) { runif(1, fossils$hmin[x], fossils$hmax[x]) })), digits = 6)
      fossils$hmax = fossils$hmin
    } else if (any(fossils$hmax != fossils$hmin)) { # use median ages
      fossils$hmin = round(unlist(lapply(1:dim(fossils)[1], function(x) { mean( c(fossils$hmin[x], fossils$hmax[x])) })), digits = 6)
      fossils$hmax = fossils$hmin
    } else {
      fossils$hmin = round(fossils$hmin, digits = 6)
      fossils$hmax = round(fossils$hmax, digits = 6)
    }
  }

  total = length(unique(fossils$sp))

  if(python){
    cat("#!/usr/bin/env python", "from numpy import * ", "",  sep = "\n", file = file, append = FALSE)

    cat("data_1=[", sep = "", file = file, append = TRUE)

    data_1 = c()
    trait1 = c()
    taxa_names = c()
    for(i in 1:total){
      sp = unique(fossils$sp)[i]
      occs = fossils[which(fossils$sp == sp),]
      times = c()
      values = c()
      for(j in 1:length(occs$sp)){
        times = c(times, occs$hmin[j])
      }
      if(!is.null(traits)) values = c(values, traits$trait[which(traits$sp == sp)])

      if(exclude.extant.singletons && length(times) == 1 && times[1] == 0) next

      if(length(times) > 0){
        data_1 = c(data_1, paste0("array([", paste(times, collapse = ","), "])"))
        if(use.sp.names)
          taxa_names = c(taxa_names, paste0("'", sp, "'"))
        else
          taxa_names = c(taxa_names, paste0("'taxa", sp, "'"))
        if(!is.null(traits)) trait1 = c(trait1, paste0("array([", paste(values, collapse = ","), "])"))
      }
    }

    # print fossil ages
    cat(paste(data_1, collapse = ",\n"), "]\n", sep = "\n", file = file, append = TRUE)

    cat("d=[data_1]", paste0("names=['", file ,"']"), "def get_data(i): return d[i]", "def get_out_name(i): return  names[i]",
        sep = "\n", file = file, append = TRUE)

    # print species names
    cat("taxa_names=[", paste(taxa_names, collapse = ","), "]\n", sep = "", file = file, append = TRUE)

    cat("def get_taxa_names(): return taxa_names", sep = "\n", file = file, append = TRUE)

    # print trait values
    if(!is.null(traits)){
      cat("\ntrait1=[", paste(trait1, collapse = ",\n"), "\n]\n\ntraits=[trait1]\ndef get_continuous(i): return traits[i]\n", sep = "", file = file, append = TRUE)
    }

  } else {

    if(is.null(traits))
      cat("Species\tStatus\tmin_age\tmax_age\n", sep = "", file = file, append = FALSE)
    else
      cat("Species\tStatus\tmin_age\tmax_age\ttraits\n", file = file, append = FALSE)

    for(i in 1:total){
      status = NULL
      sp = unique(fossils$sp)[i]
      occs = fossils[which(fossils$sp == sp),]
      # identify extant lineages
      if(any( (occs$hmin - occs$hmax < tol) & occs$hmin < tol))
        status = "extant"
      else status = "extinct"
      if(!is.null(traits)) tr = traits$trait[which(traits$sp == sp)]

      if(exclude.extant.singletons && length(occs$sp) == 1 && occs$hmax == 0) next

      for(j in 1:length(occs$sp)){
        # skip extant samples
        if(occs$hmin[j] == occs$hmax[j] && occs$hmin[j] < tol) next
        if(is.null(traits))
          cat(paste0("taxa",i), "\t", status, "\t", occs$hmin[j], "\t", occs$hmax[j], "\n", file = file, append = TRUE)
        else
          cat(paste0("taxa",i), "\t", status, "\t", occs$hmin[j], "\t", occs$hmax[j], "\t", tr, "\n", file = file, append = TRUE)
      }
    }
  }
}
rachelwarnock/fossilsim documentation built on Jan. 31, 2024, 1:04 a.m.