diff --git a/DESCRIPTION b/DESCRIPTION index e891573..2a7e18d 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -32,14 +32,15 @@ Imports: purrr, tibble, dplyr, - coda + coda, + crayon, + tidyr Suggests: devtools, ggplot2, ape, knitr, testthat, - crayon, BiocManager, ggtree (>= 3.16.0), evolnets (>= 0.0.0.9000), diff --git a/NAMESPACE b/NAMESPACE index 80ef8b5..aba6fc7 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -9,6 +9,7 @@ export(tp_installing_treeppl) export(tp_json_to_phylo) export(tp_list) export(tp_map_tree) +export(tp_mcmc_convergence) export(tp_model_library) export(tp_parse_host_rep) export(tp_parse_mcmc) @@ -23,6 +24,9 @@ export(tp_write_data) export(tp_write_model) importFrom(bnpsd,tree_reorder) importFrom(methods,is) -importFrom(phangorn,add_edge_length) -importFrom(phangorn,allCompat) +importFrom(methods,new) +importFrom(phangorn, + add_edge_length, + allCompat +) importFrom(rlang,.data) diff --git a/R/model.R b/R/model.R index c1927e4..7002aa0 100644 --- a/R/model.R +++ b/R/model.R @@ -258,8 +258,9 @@ tp_compile <- function(model, method = "mcmc", ...) { user_list <- append(list(method = method), tp_list(...)) tmp <- list_to_options(user_list) - sampler$compile_options <- tmp[["compile"]] + sampler$compile_options <- c(tmp[["compile"]], "model" = model) full_options = append(tmp[["compile"]], tmp[["runtime"]]) sampler$exe_path <- compilation(sampler$model_path, options_to_string(full_options)) + return(sampler) } diff --git a/R/post_treatment.R b/R/post_treatment.R index f824f5e..6f3fdf8 100644 --- a/R/post_treatment.R +++ b/R/post_treatment.R @@ -1,17 +1,20 @@ #' Parse TreePPL SMC output into a tidy data frame #' -#' Converts a list of parsed SMC sweeps (from \code{tp_run()}) into a +#' Converts a JSON file from an SMC analysis produced by TreePPL into a #' single tidy tibble of particles, their samples, and normalized weights. #' The function internally removes sweeps with an undefined normalizing constant. #' -#' @param treeppl_out A list of sweeps parsed from a SMC JSON output: i.e., -#' the output object of \code{tp_run()}. +#' @param json_path The full path to the (SMC) JSON file produced by TreePPL. +#' @param wide Logical. If \code{TRUE} (default), return the data frame in wide format, +#' with one column per parameter. If \code{FALSE}, return the data frame in long format, +#' with parameter names and values stored in \code{parameter} +#' and \code{sample} columns. #' #' @return A tibble with one row per particle, containing: #' \describe{ #' \item{sweep}{Sweep index.} #' \item{parameter}{Parameter name, if present in the input JSON.} -#' \item{samples}{Sampled value.} +#' \item{sample}{Sampled value.} #' \item{log_weight}{Log weight of the particle.} #' \item{norm_constant}{Log normalizing constant for the sweep.} #' \item{norm_weight}{Normalized weight, rescaled so the maximum @@ -20,21 +23,35 @@ #' #' @examples #' \dontrun{ -#' # Fit a quick CRBD model: -#' path_data <- tp_data(data_input = "crbd") -#' sampler_smc <- tp_compile( -#' model = "crbd", -#' method = "smc-apf", -#' sweeps = 2, -#' particles = 10 +#' # fit a CRBD model +#' run_smc <- tp_run( +#' data = tp_data(data_input = "crbd"), +#' sampler = tp_compile( +#' model = "crbd", +#' method = "smc-apf", +#' sweeps = 2, +#' particles = 10 +#' ) #' ) -#' mod_smc <- tp_run(sampler = sampler_smc, data = path_data) #' -#' tp_parse_smc(mod_smc) +#' # get the path to the output JSON file: +#' out_file <- list.files( +#' path = tp_tempdir(), +#' pattern = "out", +#' full.names = TRUE +#' ) +#' +#' # parse JSON to a tidy data frame +#' tp_parse_smc(json_path = out_file) #' } #' #' @export -tp_parse_smc <- function(treeppl_out) { +tp_parse_smc <- function(json_path, wide = TRUE) { + # read in the JSON file(s) + treeppl_out <- readr::read_lines(json_path) + treeppl_out <- lapply(treeppl_out, jsonlite::fromJSON, simplifyVector = FALSE) + + # parse sweeps parse_sweep <- function(sweep, sweep_id) { # remove sweeps with nan norm const if (identical(sweep$normConst, "nan")) { @@ -80,14 +97,16 @@ tp_parse_smc <- function(treeppl_out) { tibble::tibble( particle = particle_id, parameter = names(param_values), - samples = as.numeric(unlist(param_values)) + sample = as.numeric(unlist(param_values)) ) }) } else { samples_df <- purrr::imap_dfr(samples, function(s, particle_id) { + values <- as.numeric(unlist(s)) tibble::tibble( particle = particle_id, - samples = as.numeric(unlist(s)) + parameter = paste0("x", seq_along(values)), + sample = values ) }) } @@ -117,7 +136,7 @@ tp_parse_smc <- function(treeppl_out) { } # calculate the normalized weight and return - result_df |> + result_df <- result_df |> dplyr::filter(!is.infinite(.data$log_weight)) |> dplyr::mutate( total_lweight = .data$log_weight + .data$norm_constant, @@ -125,73 +144,129 @@ tp_parse_smc <- function(treeppl_out) { ) |> dplyr::select(-"total_lweight") |> dplyr::relocate("sweep") + + if (!wide) { + return(result_df) + } + result_df |> + dplyr::group_by(.data$sweep, .data$parameter) |> + dplyr::mutate(particle = dplyr::row_number()) |> + dplyr::ungroup() |> + tidyr::pivot_wider( + id_cols = c("sweep", "particle", "log_weight", "norm_constant", "norm_weight"), + names_from = "parameter", + values_from = "sample" + ) } #' Parse TreePPL MCMC output into a tidy data frame #' -#' Converts a list of parsed MCMC runs (from \code{tp_run()}) into a +#' Converts JSON file(s) produced by an MCMC analysis in TreePPL into a #' single tidy tibble of samples, one row per iteration. #' -#' @param treeppl_out A list of MCMC runs parsed from MCMC JSON output files: i.e., -#' the output object of \code{tp_run()}. +#' @param json_path The full path to the (MCMC) JSON file(s) produced by TreePPL. +#' @param wide Logical. If \code{TRUE} (default), return the data frame in wide format, +#' with one column per parameter. If \code{FALSE}, return the data frame in long format, +#' with parameter names and values stored in \code{parameter} +#' and \code{sample} columns. #' #' @return A tibble with one row per iteration, containing: #' \describe{ -#' \item{run}{Run index, corresponding to the position of the run in -#' `treeppl_out`.} +#' \item{run}{Run index.} #' \item{parameter}{Parameter name, if present in the input JSON.} -#' \item{samples}{Sampled value.} +#' \item{sample}{Sampled value.} #' } #' #' @examples #' \dontrun{ -#' # example using a CRBD model with two MCMC chains -#' path_data <- tp_data(data_input = "crbd") -#' sampler_mcmc <- tp_compile(model = "crbd", method = "mcmc", iterations = 10) -#' mod_mcmc <- tp_run( -#' sampler = sampler_mcmc, -#' data = path_data, -#' n_runs = 2 +#' +#' # Let's use a quick CRBD model as example +#' run_mcmc <- tp_run( +#' sampler = tp_compile(model = "crbd", method = "mcmc", iterations = 10), +#' data = tp_data(data_input = "crbd"), +#' n_runs = 2, # this will produce two JSON files as output +#' n_processes = 2 #' ) #' -#' tp_parse_mcmc(mod_mcmc) +#' # get the path to the output JSON file; note that the number of JSON +#' # files produced is equal to n_runs specified above +#' out_file <- list.files( +#' path = tp_tempdir(), +#' pattern = "out", +#' full.names = TRUE +#' ) +#' +#' # parse JSON to a tidy data frame +#' tp_parse_mcmc(json_path = out_file) #' } #' #' @export -tp_parse_mcmc <- function(treeppl_out) { - parse_run <- function(run, run_id) { - samples <- run$samples +tp_parse_mcmc <- function(json_path, wide = TRUE) { + parse_file <- function(file, run_id) { + lines <- readr::read_lines(file) + parsed <- lapply(lines, jsonlite::fromJSON, simplifyVector = FALSE) + + # detect if JSON was created using --incremental-printing + is_incr_print <- identical(parsed[[1]][["__constructor__"]], "ReturnType") + if (is_incr_print) { + is_sample <- purrr::map_lgl( + parsed, ~ identical(.x[["__constructor__"]], "ReturnType") + ) + n_dropped <- sum(!is_sample) + samples <- parsed[is_sample] + } else { + # regular JSON format, i.e., no --incremental-printing + samples <- parsed[[1]]$samples + } + + if (length(samples) == 0) { + stop("No samples found in file: ", file) + } + + # parameter names has_parameter_names <- is.list(samples[[1]]) && !is.null(samples[[1]][["__data__"]]) if (has_parameter_names) { + # use parameter names, if they are present in the JSON purrr::imap_dfr(samples, function(s, iteration_id) { param_values <- s[["__data__"]] tibble::tibble( run = run_id, iteration = iteration_id, parameter = names(param_values), - samples = as.numeric(unlist(param_values)) + sample = as.numeric(unlist(param_values)) ) }) } else { + # build data frame with arbitrary parameter names (e.g., x1, x2,...xn) purrr::imap_dfr(samples, function(s, iteration_id) { + values <- as.numeric(unlist(s)) tibble::tibble( run = run_id, iteration = iteration_id, - samples = as.numeric(unlist(s)) + parameter = paste0("x", seq_along(values)), + sample = values ) }) } } - result_df <- purrr::imap_dfr(treeppl_out, parse_run) - + result_df <- purrr::imap_dfr(json_path, parse_file) if (nrow(result_df) == 0) { stop("All runs failed") } - return(result_df) + if (!wide) { + return(result_df) + } + + result_df |> + tidyr::pivot_wider( + id_cols = c("run", "iteration"), + names_from = "parameter", + values_from = "sample" + ) } @@ -407,11 +482,13 @@ tp_smc_convergence <- function(treeppl_out) { #' runs are available, the upper limit of the Gelman-Rubin potential scale #' reduction factor (R-hat), using the \pkg{coda} package. #' -#' @param treeppl_out A tibble produced by \code{tp_parse_mcmc()}, with -#' columns `run`, `iteration`, `parameter`, and `samples`. The input must -#' have a `parameter` column; output produced from an unnamed TreePPL -#' return type is not supported and will raise an error, since there -#' is no reliable way to distinguish multiple parameters. +#' @param treeppl_out A tibble produced by \code{tp_parse_mcmc()}, in either long +#' or wide format. +#' +#' @details +#' Output produced from an unnamed return type in TreePPL (`.tppl` file) is not +#' supported and will raise an error, since there is no reliable way to +#' distinguish multiple parameters. #' #' @return A tibble with one row per parameter, containing: #' \describe{ @@ -424,15 +501,34 @@ tp_smc_convergence <- function(treeppl_out) { #' #' @examples #' \dontrun{ -#' d <- tp_parse_mcmc(mod_mcmc) -#' tp_mcmc_convergence(d) +#' +#' # CRBD model using MCMC +#' run_mcmc <- tp_run( +#' sampler = tp_compile(model = "crbd", method = "mcmc", iterations = 10), +#' data = tp_data(data_input = "crbd"), +#' n_runs = 2, +#' n_processes = 2 +#' ) +#' +#' # tp_run() already returns the output of tp_parse_mcmc(), so we can call +#' # tp_mcmc_convergence() directly: +#' tp_mcmc_convergence(run_mcmc) #' } #' #' @export tp_mcmc_convergence <- function(treeppl_out) { - has_parameter <- "parameter" %in% names(treeppl_out) + # if the input is in wide format + if (!"sample" %in% colnames(treeppl_out)) { + treeppl_out <- treeppl_out |> + tidyr::pivot_longer( + cols = -c("run", "iteration"), + names_to = "parameter", + values_to = "sample" + ) + } - # stop if the JSON has no parameter names + # sanity check + has_parameter <- "parameter" %in% names(treeppl_out) if (!has_parameter) { stop( "Output JSON has no parameter names.\n", @@ -441,6 +537,7 @@ tp_mcmc_convergence <- function(treeppl_out) { ) } + # coda::mcmc objects runs <- sort(unique(treeppl_out$run)) parameters <- unique(treeppl_out$parameter) @@ -449,17 +546,17 @@ tp_mcmc_convergence <- function(treeppl_out) { vals <- treeppl_out |> dplyr::filter(.data$run == r, .data$parameter == p) |> dplyr::arrange(.data$iteration) |> - dplyr::pull(.data$samples) + dplyr::pull(.data$sample) coda::mcmc(vals) }) coda::mcmc.list(chains) }) names(chains_by_parameter) <- parameters - + # ESS ess <- purrr::map_dbl(chains_by_parameter, coda::effectiveSize) - result <- tibble::tibble(parameter = names(ess), ess = ess) + # Gelman and Rubin's R-hat if (length(runs) > 1) { rhat_upper <- purrr::map_dbl(chains_by_parameter, function(chain_list) { coda::gelman.diag(chain_list)$psrf[, "Upper C.I."] @@ -475,7 +572,6 @@ tp_mcmc_convergence <- function(treeppl_out) { if (!has_parameter) { result <- dplyr::select(result, -"parameter") } - return(result) } diff --git a/R/run.R b/R/run.R index c953a00..66756f8 100644 --- a/R/run.R +++ b/R/run.R @@ -1,22 +1,36 @@ -#' Run a TreePPL program +#' Run a TreePPL sampler #' #' @description -#' Run TreePPL and return output. +#' Executes a compiled TreePPL sampler on the given data, saves the raw JSON +#' output to disk, prints a run summary, and, when a parser is available for +#' the model/method combination, parses the output into tidy data frames. #' -#' @param sampler a [treepplr::sampler_T] outputted by [treepplr::tp_compile]. -#' @param data a [base::character] with the full path to the data file in TreePPL -#' JSON format (as outputted by [treepplr::tp_data]). -#' @param dir a [base::character] with the full path to the directory where you -#' want to save the output. Default is [base::tempdir()]. -#' @param out_file_name a [base::character] with the name of the output file in -#' JSON format. Default is "out". -#' @param n_runs a [base::numeric] giving the numbers of sweeps(SMC)/runs(MCMC). -#' @param n_processes a [base::numeric], number of parallel processes to use. -#' Can't be superior to n_runs. -#' @param ... See [treepplr::tp_runtime_options] for all supported arguments. +#' @param sampler a sampler produced by \code{tp_compile()}. +#' @param data input data, produced by \code{tp_data()}. +#' @param dir the full path to the directory where +#' you want to save the output. Defaults to \code{tp_tempdir()}. +#' @param out_file_name the name of the output file in JSON format. Defaults to `"out"`. +#' @param n_runs (\code{integer}) the number of sweeps (SMC) or runs +#' (MCMC). +#' @param n_processes (\code{integer}) the number of parallel processes to use. +#' Cannot be greater than `n_runs`. +#' @param verbose (\code{logical}) Whether to print a run summary on completion. +#' Default is `TRUE`. +#' @param ... See [treepplr::tp_runtime_options()] for all supported arguments. #' +#' @details +#' If the model belongs to a category without an available parser (e.g. +#' `"host-repertoire-evolution"`, `"tree-inference"`), or if the inference +#' method is neither SMC nor MCMC, no parsing is attempted: a message is +#' printed pointing to the output directory, and the raw output file path(s) +#' are returned instead. +#' +#' @return +#' If a parser is available for the model/method combination, a parsed +#' tidy data frame of TreePPL output via `tp_parse_smc()` or `tp_parse_mcmc()`. +#' Otherwise, the full path(s) to the raw JSON output file(s), +#' along with a console message explaining that no parser is available. #' -#' @return A list of TreePPL output in parsed JSON format. #' @export #' #' @examples @@ -42,14 +56,18 @@ #' # run TreePPL #' result <- tp_run(exe_path, data_path) #' } - -tp_run <- function(sampler, - data, - dir = NULL, - out_file_name = "out", - n_runs = 1, - n_processes = 3, - ...) { +tp_run <- function( + sampler, + data, + dir = NULL, + out_file_name = "out", + n_runs = 1, + n_processes = 3, + verbose = TRUE, + ... +) { + # start time + tt <- Sys.time() if (is.null(dir)) { dir_path <- tp_tempdir() @@ -57,16 +75,18 @@ tp_run <- function(sampler, dir_path <- dir } - listFiles <- list.files(path = dir_path, - pattern = out_file_name, - full.names = TRUE) - if(length(listFiles) != 0) { + listFiles <- list.files( + path = dir_path, + pattern = out_file_name, + full.names = TRUE + ) + if (length(listFiles) != 0) { file.remove(listFiles) } - output_path <- paste0(dir_path, out_file_name, ".json") + output_path <- file.path(dir_path, paste0(out_file_name, ".json")) - #If a list have multiple time the same key + # If a list have multiple time the same key # list[[key]] will return the first key # Exemple #> lis <- list(method = "mcmc", method = "smc") @@ -96,19 +116,86 @@ tp_run <- function(sampler, future.apply::future_sapply( 1:n_runs, FUN = function(i) { - system(paste0(command, i)) + output_path_i <- file.path(dir_path, paste0(out_file_name, i, ".json")) + command_i <- paste( + "LD_LIBRARY_PATH= ", + sampler$exe_path, + data, + options_to_string(tpplc_options[["runtime"]]), + paste(">", output_path_i) + ) + system(command_i) } ) } else { system(command) } - listFiles <- list.files(path = dir_path, - pattern = out_file_name, - full.names = TRUE) + # the output (JSON) files + listFiles <- list.files( + path = dir_path, + pattern = out_file_name, + full.names = TRUE + ) + + # Run Info # + # elapsed time + et <- round(Sys.time() - tt, digits = 2) + # take method & model from sampler + mtd <- sampler$compile_options$method + mod <- sampler$compile_options$model + # output files + of <- list.files( + path = dir_path, + pattern = out_file_name, + full.names = FALSE + ) + of <- paste(of, collapse = ", ") + + # run info summary + run_info <- paste0( + crayon::bold("Analysis Summary\n"), + "-----------------------------\n", + crayon::green("Status: "), "Completed\n", + crayon::cyan("Elapsed time: "), et, " ", units(et), "\n", + crayon::cyan("Model: "), mod, "\n", + crayon::cyan("Method: "), mtd, "\n", + crayon::cyan("Output directory: "), dir_path, "\n", + crayon::cyan("Output file(s): "), of, "\n" + ) - json_out <- readr::read_lines(listFiles) |> - lapply(jsonlite::fromJSON, simplifyVector = FALSE) + # print run info summary + if (verbose) cat(run_info) - return(json_out) + # parse JSON to tidy data frames & return # + # get model category: this is needed because at the moment, we do not have parsers + # for models that return trees. NB: This is a temporary solution while we come up + # with new parsers. + mc <- tp_model_library() + mod_cat <- mc$category[match(mod, mc$model_name)] + + # model categories with unavailable parsers + no_parsers <- c( + "host-repertoire-evolution", + "tree-inference" + ) + + if (mod_cat %in% no_parsers) { + message( + "Sorry, we don't have a parser for this model and/or inference method yet.\n", + paste0("The output file(s) can be found in: ", dir_path) + ) + res <- listFiles + } else if (grepl("smc", mtd, ignore.case = TRUE)) { + res <- tp_parse_smc(listFiles) + } else if (grepl("mcmc", mtd, ignore.case = TRUE)) { + res <- tp_parse_mcmc(listFiles) + } else { + message( + "Sorry, we don't have a parser for method '", mtd, "' yet.\n", + paste0("The output file(s) can be found in: ", dir_path) + ) + res <- listFiles + } + return(res) } diff --git a/R/utils.R b/R/utils.R index dd4bdd8..fff52ad 100644 --- a/R/utils.R +++ b/R/utils.R @@ -153,23 +153,9 @@ sep <- function() { #' @return A list of model names. #' @export tp_model_library <- function() { - tp_installing_treeppl() - # make sure you get the appropriate version if you have more than one treeppl folder in the tmp - fd <- list.files("/tmp", - pattern = paste0("treeppl-", TPPLC_VERSION), - full.names = TRUE) - # go to the right treeppl folder, whatever it is called - fd <- list.files(fd, pattern = "treeppl", full.names = TRUE) - # add the rest of the path - fd <- paste0(fd, "/lib/mcore/treeppl/models") - # model names - mn <- list.files(fd, - full.names = TRUE, - recursive = TRUE, - pattern = "\\.tppl$") + mn <- tp_find("","tppl") subcategory <- grepl(".*models/([^/]+)/([^/]+)/([^/]+)\\.tppl$", mn) no_sub <- mn[!subcategory] - # results in a data frame rs <- data.frame( "category" = sub(".*models/([^/]+)/.*", "\\1", no_sub), diff --git a/man/tp_mcmc_convergence.Rd b/man/tp_mcmc_convergence.Rd index 1ec9210..21e6105 100644 --- a/man/tp_mcmc_convergence.Rd +++ b/man/tp_mcmc_convergence.Rd @@ -7,11 +7,8 @@ tp_mcmc_convergence(treeppl_out) } \arguments{ -\item{treeppl_out}{A tibble produced by \code{tp_parse_mcmc()}, with -columns \code{run}, \code{iteration}, \code{parameter}, and \code{samples}. The input must -have a \code{parameter} column; output produced from an unnamed TreePPL -return type is not supported and will raise an error, since there -is no reliable way to distinguish multiple parameters.} +\item{treeppl_out}{A tibble produced by \code{tp_parse_mcmc()}, in either long +or wide format.} } \value{ A tibble with one row per parameter, containing: @@ -28,10 +25,25 @@ Computes per-parameter effective sample size (ESS) and, when multiple runs are available, the upper limit of the Gelman-Rubin potential scale reduction factor (R-hat), using the \pkg{coda} package. } +\details{ +Output produced from an unnamed return type in TreePPL (\code{.tppl} file) is not +supported and will raise an error, since there is no reliable way to +distinguish multiple parameters. +} \examples{ \dontrun{ -d <- tp_parse_mcmc(mod_mcmc) -tp_mcmc_convergence(d) + +# CRBD model using MCMC +run_mcmc <- tp_run( +sampler = tp_compile(model = "crbd", method = "mcmc", iterations = 10), +data = tp_data(data_input = "crbd"), +n_runs = 2, +n_processes = 2 +) + +# tp_run() already returns the output of tp_parse_mcmc(), so we can call +# tp_mcmc_convergence() directly: +tp_mcmc_convergence(run_mcmc) } } diff --git a/man/tp_parse_mcmc.Rd b/man/tp_parse_mcmc.Rd index 7dba5e2..e40912e 100644 --- a/man/tp_parse_mcmc.Rd +++ b/man/tp_parse_mcmc.Rd @@ -4,37 +4,49 @@ \alias{tp_parse_mcmc} \title{Parse TreePPL MCMC output into a tidy data frame} \usage{ -tp_parse_mcmc(treeppl_out) +tp_parse_mcmc(json_path, wide = TRUE) } \arguments{ -\item{treeppl_out}{A list of MCMC runs parsed from MCMC JSON output files: i.e., -the output object of \code{tp_run()}.} +\item{json_path}{The full path to the (MCMC) JSON file(s) produced by TreePPL.} + +\item{wide}{Logical. If \code{TRUE} (default), return the data frame in wide format, +with one column per parameter. If \code{FALSE}, return the data frame in long format, +with parameter names and values stored in \code{parameter} +and \code{sample} columns.} } \value{ A tibble with one row per iteration, containing: \describe{ -\item{run}{Run index, corresponding to the position of the run in -\code{treeppl_out}.} +\item{run}{Run index.} \item{parameter}{Parameter name, if present in the input JSON.} -\item{samples}{Sampled value.} +\item{sample}{Sampled value.} } } \description{ -Converts a list of parsed MCMC runs (from \code{tp_run()}) into a +Converts JSON file(s) produced by an MCMC analysis in TreePPL into a single tidy tibble of samples, one row per iteration. } \examples{ \dontrun{ -# example using a CRBD model with two MCMC chains -path_data <- tp_data(data_input = "crbd") -sampler_mcmc <- tp_compile(model = "crbd", method = "mcmc", iterations = 10) -mod_mcmc <- tp_run( - sampler = sampler_mcmc, - data = path_data, - n_runs = 2 + +# Let's use a quick CRBD model as example +run_mcmc <- tp_run( +sampler = tp_compile(model = "crbd", method = "mcmc", iterations = 10), +data = tp_data(data_input = "crbd"), +n_runs = 2, # this will produce two JSON files as output +n_processes = 2 +) + +# get the path to the output JSON file; note that the number of JSON +# files produced is equal to n_runs specified above +out_file <- list.files( + path = tp_tempdir(), + pattern = "out", + full.names = TRUE ) -tp_parse_mcmc(mod_mcmc) +# parse JSON to a tidy data frame +tp_parse_mcmc(json_path = out_file) } } diff --git a/man/tp_parse_smc.Rd b/man/tp_parse_smc.Rd index a4003ef..eac4dec 100644 --- a/man/tp_parse_smc.Rd +++ b/man/tp_parse_smc.Rd @@ -4,18 +4,22 @@ \alias{tp_parse_smc} \title{Parse TreePPL SMC output into a tidy data frame} \usage{ -tp_parse_smc(treeppl_out) +tp_parse_smc(json_path, wide = TRUE) } \arguments{ -\item{treeppl_out}{A list of sweeps parsed from a SMC JSON output: i.e., -the output object of \code{tp_run()}.} +\item{json_path}{The full path to the (SMC) JSON file produced by TreePPL.} + +\item{wide}{Logical. If \code{TRUE} (default), return the data frame in wide format, +with one column per parameter. If \code{FALSE}, return the data frame in long format, +with parameter names and values stored in \code{parameter} +and \code{sample} columns.} } \value{ A tibble with one row per particle, containing: \describe{ \item{sweep}{Sweep index.} \item{parameter}{Parameter name, if present in the input JSON.} -\item{samples}{Sampled value.} +\item{sample}{Sampled value.} \item{log_weight}{Log weight of the particle.} \item{norm_constant}{Log normalizing constant for the sweep.} \item{norm_weight}{Normalized weight, rescaled so the maximum @@ -23,23 +27,32 @@ total log weight across all particles is 1.} } } \description{ -Converts a list of parsed SMC sweeps (from \code{tp_run()}) into a +Converts a JSON file from an SMC analysis produced by TreePPL into a single tidy tibble of particles, their samples, and normalized weights. The function internally removes sweeps with an undefined normalizing constant. } \examples{ \dontrun{ -# Fit a quick CRBD model: -path_data <- tp_data(data_input = "crbd") -sampler_smc <- tp_compile( - model = "crbd", - method = "smc-apf", - sweeps = 2, - particles = 10 +# fit a CRBD model +run_smc <- tp_run( + data = tp_data(data_input = "crbd"), + sampler = tp_compile( + model = "crbd", + method = "smc-apf", + sweeps = 2, + particles = 10 + ) +) + +# get the path to the output JSON file: +out_file <- list.files( + path = tp_tempdir(), + pattern = "out", + full.names = TRUE ) -mod_smc <- tp_run(sampler = sampler_smc, data = path_data) -tp_parse_smc(mod_smc) +# parse JSON to a tidy data frame +tp_parse_smc(json_path = out_file) } } diff --git a/man/tp_run.Rd b/man/tp_run.Rd index b417891..84436c4 100644 --- a/man/tp_run.Rd +++ b/man/tp_run.Rd @@ -2,7 +2,7 @@ % Please edit documentation in R/run.R \name{tp_run} \alias{tp_run} -\title{Run a TreePPL program} +\title{Run a TreePPL sampler} \usage{ tp_run( sampler, @@ -11,33 +11,48 @@ tp_run( out_file_name = "out", n_runs = 1, n_processes = 3, + verbose = TRUE, ... ) } \arguments{ -\item{sampler}{a \link{sampler_T} outputted by \link{tp_compile}.} +\item{sampler}{a sampler produced by \code{tp_compile()}.} -\item{data}{a \link[base:character]{base::character} with the full path to the data file in TreePPL -JSON format (as outputted by \link{tp_data}).} +\item{data}{input data, produced by \code{tp_data()}.} -\item{dir}{a \link[base:character]{base::character} with the full path to the directory where you -want to save the output. Default is \code{\link[base:tempdir]{base::tempdir()}}.} +\item{dir}{the full path to the directory where +you want to save the output. Defaults to \code{tp_tempdir()}.} -\item{out_file_name}{a \link[base:character]{base::character} with the name of the output file in -JSON format. Default is "out".} +\item{out_file_name}{the name of the output file in JSON format. Defaults to \code{"out"}.} -\item{n_runs}{a \link[base:numeric]{base::numeric} giving the numbers of sweeps(SMC)/runs(MCMC).} +\item{n_runs}{(\code{integer}) the number of sweeps (SMC) or runs +(MCMC).} -\item{n_processes}{a \link[base:numeric]{base::numeric}, number of parallel processes to use. -Can't be superior to n_runs.} +\item{n_processes}{(\code{integer}) the number of parallel processes to use. +Cannot be greater than \code{n_runs}.} -\item{...}{See \link{tp_runtime_options} for all supported arguments.} +\item{verbose}{(\code{logical}) Whether to print a run summary on completion. +Default is \code{TRUE}.} + +\item{...}{See \code{\link[=tp_runtime_options]{tp_runtime_options()}} for all supported arguments.} } \value{ -A list of TreePPL output in parsed JSON format. +If a parser is available for the model/method combination, a parsed +tidy data frame of TreePPL output via \code{tp_parse_smc()} or \code{tp_parse_mcmc()}. +Otherwise, the full path(s) to the raw JSON output file(s), +along with a console message explaining that no parser is available. } \description{ -Run TreePPL and return output. +Executes a compiled TreePPL sampler on the given data, saves the raw JSON +output to disk, prints a run summary, and, when a parser is available for +the model/method combination, parses the output into tidy data frames. +} +\details{ +If the model belongs to a category without an available parser (e.g. +\code{"host-repertoire-evolution"}, \code{"tree-inference"}), or if the inference +method is neither SMC nor MCMC, no parsing is attempted: a message is +printed pointing to the output directory, and the raw output file path(s) +are returned instead. } \examples{ \dontrun{ diff --git a/man/treepplr-package.Rd b/man/treepplr-package.Rd index 4306e7c..aa86ce6 100644 --- a/man/treepplr-package.Rd +++ b/man/treepplr-package.Rd @@ -24,7 +24,7 @@ Authors: \itemize{ \item Mariana P Braga \email{mpiresbr@gmail.com} (\href{https://orcid.org/0000-0002-1253-2536}{ORCID}) \item Tim Virgoulay - \item Stenio I A Foerster + \item Stenio Foerster } } diff --git a/tests/testthat/test-run.R b/tests/testthat/test-run.R index aea9cbc..ccde83c 100644 --- a/tests/testthat/test-run.R +++ b/tests/testthat/test-run.R @@ -5,36 +5,57 @@ require(crayon) cat(crayon::yellow("\nTest-run : Running TreePPL.\n")) -test_that("Test-run_1a : tp_run", { - cat("\tTest-run_1a : tp_run \n") - - sampler <- treepplr::tp_compile("crbd", method = "smc-apf", particles = 2) - data <- treepplr::tp_data("crbd") - - result <- treepplr::tp_run(sampler, data, sweeps = 1) - - expect_equal(2, length(result[[1]]$samples)) - +test_that("Test-run_1a : tp_run SMC", { + cat("\tTest-run_1a : tp_run SMC \n") + run_smc <- tp_run( + sampler = tp_compile(model = "crbd", method = "smc-apf", sweeps = 2, particles = 5), + data = tp_data(data_input = "crbd"), + verbose = FALSE + ) + expect_equal(2, length(unique(run_smc$sweep))) }) -test_that("Test-run_1a : tp_run custom name", { - cat("\tTest-run_1a : tp_run \n") - - sampler <- treepplr::tp_compile("crbd", method = "smc-apf", particles = 2) - data <- treepplr::tp_data("crbd") - - result <-treepplr::tp_run(sampler, data, sweeps = 1, out_file_name = "test_out", particles = 5) - - expect_equal(5, length(result[[1]]$samples)) +test_that("Test-run_1b : tp_run MCMC", { + cat("\tTest-run_1b : tp_run MCMC \n") + run_mcmc <- tp_run( + sampler = tp_compile(model = "crbd", method = "mcmc", iterations = 10), + data = tp_data(data_input = "crbd"), + n_runs = 2, + n_processes = 2, + verbose = FALSE + ) + expect_equal(2, length(unique(run_mcmc$run))) }) -test_that("Test-run_1c : tp_run threading", { - cat("\tTest-run_1c : tp_run \n") - - sampler <- treepplr::tp_compile("crbd", method = "smc-apf", particles = 2) - data <- treepplr::tp_data("crbd") +test_that("Test-run_1c : tp_run custom_name", { + cat("\tTest-run_1c : tp_run custom_name \n") + run_smc <- tp_run( + sampler = tp_compile(model = "crbd", method = "smc-apf", sweeps = 2, particles = 5), + data = tp_data(data_input = "crbd"), + out_file_name = "test_out", + verbose = FALSE + ) + expect_equal(2, length(unique(run_smc$sweep))) +}) - result <-treepplr::tp_run(sampler, data, sweeps = 1, out_file_name = "test_out", particles = 5, , n_runs = 10) +test_that("Test-run_1d : tp_run threading", { + cat("\tTest-run_1d : tp_run threading \n") + run_smc <- tp_run( + sampler = tp_compile(model = "crbd", method = "smc-apf", sweeps = 2, particles = 5), + data = tp_data(data_input = "crbd"), + n_processes = 2, + verbose = FALSE + ) + expect_equal(2, length(unique(run_smc$sweep))) +}) - expect_equal(10, length(result)) +test_that("Test-run_1e : tp_run no_parser", { + cat("\tTest-run_1e : tp_run no_parser \n") + run_smc <- tp_run( + sampler = tp_compile(model = "tree_inference", method = "smc-apf", sweeps = 2, particles = 5), + data = tp_data(data_input = "tree_inference"), + n_processes = 2, + verbose = FALSE + ) + expect_true(is.character(run_smc)) }) diff --git a/vignettes/coin-example.Rmd b/vignettes/coin-example.Rmd index 41faa82..8989ca5 100644 --- a/vignettes/coin-example.Rmd +++ b/vignettes/coin-example.Rmd @@ -14,17 +14,17 @@ knitr::opts_chunk$set( collapse = TRUE, comment = "#>", warning = FALSE, - message = FALSE + message = FALSE, eval = FALSE ) options(rmarkdown.html_vignette.check_title = FALSE) ``` + This tutorial describes how to analyze a simple coin-flipping model using *treepplr*. We assume that the probability of obtaining heads with our coin is `p`, and that we have a set of flips informing us about the value of `p`. -In a Bayesian analysis of this problem, we need to specify a prior probability distribution for the value of `p`. -Here, we will assume that `p` is drawn from a `Beta(2,2)` distribution. +In a Bayesian analysis of this problem, we need to specify a prior probability distribution for the value of `p`. Here, we will assume that `p` is drawn from a `Beta(2,2)` distribution. ## Load the required R packages @@ -39,8 +39,7 @@ library(readr) ## Understanding the coin model -The coin model is one of the most basic models in the TreePPL model library. You -can find all available models and some information about them [here](https://treeppl.org/docs/model-library). +The coin model is one of the most basic models in the TreePPL model library. You can find all available models and some information about them [here](https://treeppl.org/docs/model-library). If you want to look at the TreePPL code here in R, you can use the following functions: @@ -54,7 +53,7 @@ readr::read_file(sampler$model_path) The main part of the model is defined in a function called `coinModel`: -``` +``` model function coinModel(coinflips: Bool[]) => Real { // Uncomment if you want to test the input //printLn("Input:"); @@ -69,66 +68,43 @@ model function coinModel(coinflips: Bool[]) => Real { } ``` -The definition of this function is preceded by the keywords `model function`. -All model scripts must have exactly one model function. +The definition of this function is preceded by the keywords `model function`. All model scripts must have exactly one model function. -The model function takes as input argument(s) the observed data that we wish to condition on. -In our case, the data are represented by a sequence (vector or array) of Boolean (`TRUE`/`FALSE`) values. +The model function takes as input argument(s) the observed data that we wish to condition on. In our case, the data are represented by a sequence (vector or array) of Boolean (`TRUE`/`FALSE`) values. -Note that TreePPL uses type annotation of input variables in the form of `: `. -The square brackets `[]` are used to denote a sequence type, that is, `coinflips: Bool[]` tells us that -the model function takes a single argument by the name of `coinflips`, and the data type is a sequence of Booleans. +Note that TreePPL uses type annotation of input variables in the form of `: `. The square brackets `[]` are used to denote a sequence type, that is, `coinflips: Bool[]` tells us that the model function takes a single argument by the name of `coinflips`, and the data type is a sequence of Booleans. -The return type of a function is specified using the format `=> `. -Our model function returns the value of `p`. -In general, the model function should return the model parameters for which we are interested in inferring the posterior distribution. +The return type of a function is specified using the format `=> `. Our model function returns the value of `p`. In general, the model function should return the model parameters for which we are interested in inferring the posterior distribution. -The model function starts with a few statements that are commented out by having `//` put in front of them. -The code in these lines can be used to print the value of the `coinflips` argument. +The model function starts with a few statements that are commented out by having `//` put in front of them. The code in these lines can be used to print the value of the `coinflips` argument. -The statement `assume p ~ Beta(2.0,2.0)` specifies the prior probability distribution for the parameter `p` in our model. -TreePPL provides a number of built-in probability distributions; they all have names starting with a capital letter, like the `Beta` distribution used here. +The statement `assume p ~ Beta(2.0,2.0)` specifies the prior probability distribution for the parameter `p` in our model. TreePPL provides a number of built-in probability distributions; they all have names starting with a capital letter, like the `Beta` distribution used here. -The `assume` statement is followed by a loop over the input data sequence, conditioning the simulation on each observed value (heads/`TRUE` or tails/`FALSE`). -Specifically, this is achieved by calling the help function `flip`. +The `assume` statement is followed by a loop over the input data sequence, conditioning the simulation on each observed value (heads/`TRUE` or tails/`FALSE`). Specifically, this is achieved by calling the help function `flip`. The `flip` function is defined as follows: -``` +``` function flip(datapoint: Bool, probability: Real) { observe datapoint ~ Bernoulli(probability); } ``` -`Bernoulli` is a probability distribution on a binary outcome space, -represented as a Boolean variable taking the values `TRUE` or `FALSE`. -This is the appropriate probability distribution for coin flipping. +`Bernoulli` is a probability distribution on a binary outcome space, represented as a Boolean variable taking the values `TRUE` or `FALSE`. This is the appropriate probability distribution for coin flipping. -The single parameter of the `Bernoulli` distribution, called `probability` in the `flip` function, -is assumed by convention to be the probability of obtaining the outcome `TRUE`. -In our case, this would be the probability of obtaining heads. +The single parameter of the `Bernoulli` distribution, called `probability` in the `flip` function, is assumed by convention to be the probability of obtaining the outcome `TRUE`. In our case, this would be the probability of obtaining heads. -The `observe` statement weights the simulation with the likelihood of observing -a particlar data point from the Bernoulli distribution. - -This completes the TreePPL description of the coin flipping model. -All TreePPL model descriptions essentially follow the same format. -For a more complete coverage of the TreePPL language, see the language overview in the [TreePPL online documentation](https://treeppl.org/docs/). +The `observe` statement weights the simulation with the likelihood of observing a particlar data point from the Bernoulli distribution. +This completes the TreePPL description of the coin flipping model. All TreePPL model descriptions essentially follow the same format. For a more complete coverage of the TreePPL language, see the language overview in the [TreePPL online documentation](https://treeppl.org/docs/). ## Model compilation and inference strategy -TreePPL offers a variety of inference methods. Different methods work best for -different models. Here, we will use sequential Monte Carlo (SMC), specifically -the bootstrap particle filter version (the `method=smc-bpf` option). +TreePPL offers a variety of inference methods. Different methods work best for different models. Here, we will use sequential Monte Carlo (SMC), specifically the bootstrap particle filter version (the `method=smc-bpf` option). -The function `tp_compile()` has many optional arguments that allow you to select among the inference -methods supported by TreePPL, and setting relevant options for each one of them. -For an up-to-date description of available inference strategies supported, -see `tp_compile_options()`. +The function `tp_compile()` has many optional arguments that allow you to select among the inference methods supported by TreePPL, and setting relevant options for each one of them. For an up-to-date description of available inference strategies supported, see `tp_compile_options()`. -Now let's compile the model to en executable that also contains the necessary machinery -to run the chosen inference method. +Now let's compile the model to en executable that also contains the necessary machinery to run the chosen inference method. ```{r, eval = FALSE} exe_path <- tp_compile(model = "coin", method = "smc-bpf", particles = 5000) @@ -136,9 +112,7 @@ exe_path <- tp_compile(model = "coin", method = "smc-bpf", particles = 5000) ## Data -Now we can start analyzing the model by inferring the value of `p` given some observed sequence of coin flips. -To do this, we need to provide the observations in a suitable format. -To load the example data provided in the *treepplr* package, use: +Now we can start analyzing the model by inferring the value of `p` given some observed sequence of coin flips. To do this, we need to provide the observations in a suitable format. To load the example data provided in the *treepplr* package, use: ```{r, eval=FALSE} data <- tp_data(data_input = "coin") @@ -150,31 +124,18 @@ We can look at the structure of the input data using: jsonlite::fromJSON(data) ``` -```` +``` $coinflips [1] TRUE TRUE TRUE FALSE TRUE FALSE FALSE TRUE TRUE FALSE FALSE FALSE TRUE FALSE TRUE FALSE FALSE TRUE FALSE FALSE -```` - - -The TreePPL compiler requires data in JSON format. -The conversion from R variables to appropriate JSON code understood by TreePPL -is done automatically by *treepplr* for supported data types. -For instance, logical vectors in R are converted to sequences of Booleans (`Bool[]`) -in TreePPL. - -In R, the data are collected in a list. -Each element in the list is named using the corresponding argument name expected by the TreePPL model function. -The arguments can be given in any order in the list. -In our case, the model function only takes one argument called `coinflips`, -so the R list only contains one element named `coinflips`. +``` +The TreePPL compiler requires data in JSON format. The conversion from R variables to appropriate JSON code understood by TreePPL is done automatically by *treepplr* for supported data types. For instance, logical vectors in R are converted to sequences of Booleans (`Bool[]`) in TreePPL. +In R, the data are collected in a list. Each element in the list is named using the corresponding argument name expected by the TreePPL model function. The arguments can be given in any order in the list. In our case, the model function only takes one argument called `coinflips`, so the R list only contains one element named `coinflips`. ## Run TreePPL -Now we can run the TreePPL program, inferring the posterior distribution -of `p` conditioned on the input data. -This is done using the `tp_run` function. +Now we can run the TreePPL program, inferring the posterior distribution of `p` conditioned on the input data. This is done using the `tp_run` function. Let's run 10 sweeps (10 SMC runs, if you wish). @@ -188,40 +149,26 @@ output_list <- readRDS("rdata/coin/output_coin.rds") The run should take a few seconds to complete depending on your machine. -In general, the TreePPL inference strategies can be described as nested approaches where -the outermost shell defines the character of the inference output. If the outer shell is -SMC, as is the case here, the returned object will be a nested R list. -Each entry in the outermost list layer will correspond to a sweep. -Each sweep will contain three values: the normalizing constant estimated in that sweep, -each of the returned values from the model function (one returned value for each SMC particle), -and the likelihood weight of each returned value (one weight for each particle). +In general, the TreePPL inference strategies can be described as nested approaches where the outermost shell defines the character of the inference output. If the outer shell is SMC, as is the case here, the returned object will be a nested R list. Each entry in the outermost list layer will correspond to a sweep. Each sweep will contain three values: the normalizing constant estimated in that sweep, each of the returned values from the model function (one returned value for each SMC particle), and the likelihood weight of each returned value (one weight for each particle). ## Plot the posterior distribution -In SMC we can assess the quality of the inference by running several sweeps -and comparing their normalizing constants to test if they generated similar -estimates of the posterior distribution. A popular argument from the SMC -literature suggests that the SMC estimate is accurate if the variance of the -estimates of the normalizing constant across sweeps is lower than 1. +In SMC we can assess the quality of the inference by running several sweeps and comparing their normalizing constants to test if they generated similar estimates of the posterior distribution. A popular argument from the SMC literature suggests that the SMC estimate is accurate if the variance of the estimates of the normalizing constant across sweeps is lower than 1. -To check the variance of the normalizing constant, use the `tp_smc_convergence()` function. -But first, we'll use the `tp_parse_smc()` function to convert the results returned by TreePPL into an R data frame -of appropriately weighted values, taking both the particle weights and normalizing constants -(the sweep weights, if you wish) into account. Note that both the particle weights and -normalizing constants are given in log units. +To check the variance of the normalizing constant, use the `tp_smc_convergence()` function. But first, we'll use the `tp_parse_smc()` function to convert the results returned by TreePPL into an R data frame of appropriately weighted values, taking both the particle weights and normalizing constants (the sweep weights, if you wish) into account. Note that both the particle weights and normalizing constants are given in log units. -```{r} -output <- tp_parse_smc(output_list) +```{r, eval=FALSE} +output <- tp_parse_smc(output_list) # update this according to new parsers tp_smc_convergence(output) ``` -It seems that our run provides a quite accurate estimate of the posterior distribution, -as the variance is much smaller than 1.0. +It seems that our run provides a quite accurate estimate of the posterior distribution, as the variance is much smaller than 1.0. It is also easy to plot the sampled values. -```{r, fig.height=5, fig.width=5} +```{r, fig.height=5, fig.width=5, eval=FALSE} +# update this according to new parsers ggplot2::ggplot(output, ggplot2::aes(samples, weight = norm_weight)) + ggplot2::geom_histogram(ggplot2::aes(y = ggplot2::after_stat(density)), col = "white", fill = "lightblue", binwidth=0.01) + diff --git a/vignettes/crbd-example.Rmd b/vignettes/crbd-example.Rmd index e93827d..d8b63aa 100644 --- a/vignettes/crbd-example.Rmd +++ b/vignettes/crbd-example.Rmd @@ -14,7 +14,7 @@ knitr::opts_chunk$set( collapse = TRUE, comment = "#>", warning = FALSE, - message = FALSE + message = FALSE, eval = FALSE # update this according new parsers ) options(rmarkdown.html_vignette.check_title = FALSE) @@ -35,13 +35,8 @@ output_list <- tp_run(exe_path, data_path, sweeps = 4) output_list <- readRDS("rdata/crbd/output_crbd.rds") ``` - ```{r} output <- tp_parse_smc(output_list) tp_smc_convergence(output) ``` - - - -