Get lambda and mu estimates of twin tree

[richelbilderbeek]

Needed for #117.

[richelbilderbeek]

#' Simulate a Birth-Death (BD) twin tree from the true phylogeny
# ...
sim_bd_twin_tree <- function(
  true_phylogeny,
  method = "random_tree",
  n_replicates = 1e4,
  os = rappdirs::app_dir()$os
) {
  beastier::check_os(os)
  methods <- c("random_tree", "max_clade_cred", "max_likelihood")
  if (!method %in% methods) {
    stop(
      "'method' not in the supported methods. \n",
      "Supported methods: ", methods, ". \n",
      "Actual value: ", method
    )
  }
  phylogeny <- true_phylogeny

  age  <- beautier::get_crown_age(phylogeny)
  phylo_brts <- sort(
    pirouette::convert_tree2brts(phylogeny),
    decreasing = TRUE
  )
  n_tips <- ape::Ntip(phylogeny)
  soc <- 1 + n_tips - length(phylo_brts)
  testit::assert(soc == 1 | soc == 2)
  difference <- (log(n_tips) - log(soc)) / age
  mu <- 0.1
  lambda <- mu + difference

  if (os != "win") {
    sink(tempfile())
    bd_pars <- DDD::bd_ML(
      brts = sort(phylo_brts, decreasing = TRUE),
      cond = 1, # conditioning on stem or crown age # nolint
      initparsopt = c(lambda, mu),
      idparsopt = 1:2,
      missnumspec = 0,
      tdmodel = 0,
      btorph = 1,
      soc = soc
    )
    sink()
  } else {
    x <- utils::capture.output({
        bd_pars <- DDD::bd_ML(
          brts = sort(phylo_brts, decreasing = TRUE),
          cond = 1,
          initparsopt = c(lambda, mu),
          idparsopt = 1:2,
          missnumspec = 0,
          tdmodel = 0,
          btorph = 1,
          soc = soc
        )
    }
    )
    rm(x)
  }

  lambda_bd <- as.numeric(unname(bd_pars[1]))
  mu_bd <- as.numeric(unname(bd_pars[2]))
  testit::assert(beautier::is_one_double(lambda_bd))
  testit::assert(beautier::is_one_double(mu_bd))

Then added this:

  # Write to temp file
  readr::write_lines(
    paste("lambda_bd: ", lambda_bd, "mu_bd: ", mu_bd),
    "~/estimates.txt"
  )

Then business as usual:

  # generate bd branching times from the inferred parameters
  bd_brts0 <- pirouette::create_twin_branching_times(
    phylogeny = phylogeny,
    lambda = lambda_bd,
    mu = mu_bd,
    n_replicates = n_replicates,
    method = method
  )

  pirouette::combine_brts_and_topology(
    brts = bd_brts0,
    tree = phylogeny
  )
}

[richelbilderbeek]

richel@N141CU:~$ cat estimates.txt 
lambda_bd:  0.199967443181565 mu_bd:  0.000157467752061131

[richelbilderbeek]

Done!