suppressPackageStartupMessages({
  library(tidyverse)
  library(RERconverge)
})

master_tree <- read.tree(file="<<<dep master_tree>>>")
gene_trees <- readTrees("<<<dep gene_tree_set>>>", max.read=<<<max_read>>>, masterTree=master_tree, minSpecs=<<<min_specs>>>)
is_continuous = <<<template_bool continuous>>>

relative_rates <- getAllResiduals(
  gene_trees,
  transform = "<<<template_transform transform>>>",
  weighted = <<<template_bool weighted>>>,
  scale = <<<template_bool scale>>>
)

load_species_file <- function(path) {
    read_table(path, col_names = "species") %>% pull(species)
}

load_traits_file <- function(path, is_continuous) {
  if (is_continuous)
    read_table(path, col_names = c('species', 'value')) %>% deframe()
  else
    load_species_file(path)
}

phenotypes <- load_traits_file("<<<dep phenotypes>>>", is_continuous)

build_traits_paths <-
  function(phenotypes,
           gene_trees,
           is_continuous,
           clade)
  {
    if (is_continuous) {
      char2Paths(phenotypes, gene_trees)
    } else {
      foreground2Tree(phenotypes,
                      gene_trees,
                      clade = clade,
                      plot=F # FIXME : plot=T raises cryptic error
                      ) %>%
        tree2Paths(gene_trees)
    }
  }

clades = c("ancestral", "terminal", "all")

results = lapply(clades, FUN=function(clade){

  traits_paths <- build_traits_paths(
    phenotypes,
    gene_trees,
    is_continuous,
    clade
  )

  res <- if (is_continuous) {
    correlateWithContinuousPhenotype(
      relative_rates,
      traits_paths,
    )
  } else {
    correlateWithBinaryPhenotype(
      relative_rates,
      traits_paths,
    )
  }

  table <- res %>% 
    mutate(gene_id=names(gene_trees$trees)) %>% 
    select(gene_id, everything())

  list(
    table=table,
    traits_paths=traits_paths
  )
})

names(results) = clades

results = results %>% 
  list_modify(
    relative_rates = relative_rates,
    master_tree = master_tree,
    gene_trees = gene_trees,
    phenotypes = phenotypes
  )

saveRDS(results, file="<<<dest>>>")