R/main.R

Defines functions sls_main

Documented in sls_main

#' @title Run the full sls routine
#' @description Simulate an sls process and infer parameters maximizing the
#' likelihood(s) for given function(s)
#' @inheritParams default_params_doc
#' @details mle inference
#' @export
sls_main <- function(
  sim_pars,
  cond = 1,
  l_2 = sim_get_standard_l_2(
    crown_age = 5,
    shift_time = 2
  ),
  seed,
  start_pars = c(0.2, 0.1, 0.2, 0.1),
  optim_ids = rep(TRUE, length(start_pars)),
  models = get_logliks(),
  project_folder = NULL,
  verbose = FALSE
) {
  # specific set up
  lambdas <- sim_pars[c(1, 3)]
  mus <- sim_pars[c(2, 4)]
  ks <- c(Inf, Inf)
  t_0s <- l_2$birth_time
  n_0 <- l_2$n_0[1]
  n_0s <- l_2$n_0

  # generic set up
  function_names <- get_function_names( # nolint internal function
    models = models
  )
  model_names <- get_model_names( # nolint internal function
    function_names = function_names,
    verbose = verbose
  )

  # simulate
  set.seed(seed)
  sim <- sls_sim(
    lambdas = lambdas,
    mus = mus,
    ks = ks,
    cond = cond,
    l_2 = l_2
  )
  brts <- sim$brts
  print_info(brts = brts, n_0 = n_0, cond = cond, verbose = verbose) # nolint internal function
  if (!is.list(brts)) {
    tips <- (n_0s[1] - 1) + length(brts)
  } else {
    tips <- rep(NA, length(brts))
    for (i in seq_along(brts)) {
      tips[i] <- n_0s[i] - 1 + length(brts[[i]])
    }
  }

  # maximum likelihood
  for (m in seq_along(models)) {
    if (verbose == FALSE) {
      if (rappdirs::app_dir()$os != "win") {
        sink(file.path(rappdirs::user_cache_dir(), "ddd"))
      } else {
        sink(rappdirs::user_cache_dir())
      }
    }
    mle <- sls_ml(
      loglik_function = get(function_names[m]),
      brts = brts,
      cond = cond,
      n_0 = n_0,
      start_pars = start_pars,
      optim_ids = optim_ids,
      true_pars = sim_pars,
      verbose = verbose
    )
    if (verbose == FALSE) {
      sink()
    }
  }

  # format results
  results <- cbind(
    matrix(
      sim_pars,
      nrow = length(models),
      ncol = length(start_pars),
      byrow = TRUE
    ),
    mle,
    seed,
    cond,
    n_0,
    matrix(
      t_0s,
      nrow = length(models),
      ncol = length(t_0s),
      byrow = TRUE
    ),
    matrix(
      tips,
      nrow = length(models),
      ncol = length(tips),
      byrow = TRUE
    ),
    matrix(
      optim_ids,
      nrow = length(models),
      ncol = length(optim_ids),
      byrow = TRUE
    ),
    model = model_names
  )
  if (length(t_0s) > 1) {
    t_0s_label <- paste0("t_0_", 1:length(t_0s))
  } else {
    t_0s_label <- "t_0"
  }
  if (length(tips) > 1) {
    tips_label <- paste0("tips_", 1:length(tips))
  } else {
    tips_label <- "tips"
  }
  colnames(results) <- c(
    paste0("sim_", colnames(mle[1:length(start_pars)])),
    colnames(mle),
    "seed",
    "cond",
    "n_0",
    t_0s_label,
    tips_label,
    paste0("optim_", colnames(mle[1:length(start_pars)])),
    "model"
  )
  rownames(results) <- NULL
  results <- data.frame(results)

  # save data
  main_save_files( # nolint internal function
    project_folder = project_folder,
    sim_pars = sim_pars,
    optim_ids = optim_ids,
    cond = cond,
    n_0 = n_0,
    t_0s = t_0s,
    seed = seed,
    sim = sim,
    results = results
  )
  print_info(brts = brts, n_0 = n_0, cond = cond, verbose = verbose) # nolint internal function
  results
}
richelbilderbeek/splendid documentation built on May 20, 2019, 9:42 a.m.