From d98be09279f836aa3f4b05b864d5a3088d8195c7 Mon Sep 17 00:00:00 2001 From: foersterst Date: Tue, 1 Sep 2026 16:04:41 +0200 Subject: [PATCH 01/10] nicer output for tp_run(); includes updates to tests and documentation of all functions involved. --- DESCRIPTION | 4 +- NAMESPACE | 7 +- R/model.R | 3 +- R/post_treatment.R | 79 +++++++++++++-------- R/run.R | 142 +++++++++++++++++++++++++++++--------- man/tp_parse_mcmc.Rd | 35 ++++++---- man/tp_parse_smc.Rd | 34 +++++---- man/tp_run.Rd | 39 +++++++---- man/treepplr-package.Rd | 2 +- tests/testthat/test-run.R | 68 +++++++++++------- 10 files changed, 280 insertions(+), 133 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index e891573..1543074 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -32,14 +32,14 @@ Imports: purrr, tibble, dplyr, - coda + coda, + crayon 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..3d2ac0f 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,8 @@ export(tp_write_data) export(tp_write_model) importFrom(bnpsd,tree_reorder) importFrom(methods,is) -importFrom(phangorn,add_edge_length) -importFrom(phangorn,allCompat) +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..4fb7659 100644 --- a/R/post_treatment.R +++ b/R/post_treatment.R @@ -1,11 +1,10 @@ #' 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. #' #' @return A tibble with one row per particle, containing: #' \describe{ @@ -20,21 +19,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) { + # 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")) { @@ -130,36 +143,48 @@ tp_parse_smc <- function(treeppl_out) { #' 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. #' #' @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.} #' } #' #' @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) { +tp_parse_mcmc <- function(json_path) { + # read in the JSON file(s) + treeppl_out <- readr::read_lines(json_path) + treeppl_out <- lapply(treeppl_out, jsonlite::fromJSON, simplifyVector = FALSE) + + # parse mcmc runs parse_run <- function(run, run_id) { samples <- run$samples has_parameter_names <- is.list(samples[[1]]) && !is.null(samples[[1]][["__data__"]]) diff --git a/R/run.R b/R/run.R index c953a00..46b4b0d 100644 --- a/R/run.R +++ b/R/run.R @@ -1,22 +1,34 @@ -#' 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 ... 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 +54,17 @@ #' # 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, + ... +) { + # start time + tt <- Sys.time() if (is.null(dir)) { dir_path <- tp_tempdir() @@ -57,16 +72,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") - #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") @@ -103,12 +120,71 @@ tp_run <- function(sampler, 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("Time elapsed: "), 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 + 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[mc$model_name == mod, ]$category + + # 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/man/tp_parse_mcmc.Rd b/man/tp_parse_mcmc.Rd index 7dba5e2..53f934b 100644 --- a/man/tp_parse_mcmc.Rd +++ b/man/tp_parse_mcmc.Rd @@ -4,37 +4,44 @@ \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) } \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.} } \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.} } } \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..f930bf1 100644 --- a/man/tp_parse_smc.Rd +++ b/man/tp_parse_smc.Rd @@ -4,11 +4,10 @@ \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) } \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.} } \value{ A tibble with one row per particle, containing: @@ -23,23 +22,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 + ) ) -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) } } diff --git a/man/tp_run.Rd b/man/tp_run.Rd index b417891..c046d6f 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, @@ -15,29 +15,40 @@ tp_run( ) } \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{...}{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..2f38b42 100644 --- a/tests/testthat/test-run.R +++ b/tests/testthat/test-run.R @@ -5,36 +5,52 @@ 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") + ) + 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 + ) + 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" + ) + 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 + ) + 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 + ) + expect_true(is.character(run_smc)) }) From e37cd454c6819f716aa32bc26a8d78e1ab0e83bd Mon Sep 17 00:00:00 2001 From: foersterst Date: Tue, 1 Sep 2026 18:28:52 +0200 Subject: [PATCH 02/10] set eval = FALSE in opts_chunk$set() for vignettes, just to pass the check. --- vignettes/coin-example.Rmd | 121 +++++++++++-------------------------- vignettes/crbd-example.Rmd | 7 +-- 2 files changed, 35 insertions(+), 93 deletions(-) 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) ``` - - - - From cb3569ea422989419298145c752503de4daf1f00 Mon Sep 17 00:00:00 2001 From: foersterst Date: Wed, 2 Sep 2026 15:53:54 +0200 Subject: [PATCH 03/10] wide df by default; inlcudes updates to tp_mcmc_convergence & updated docs --- DESCRIPTION | 3 +- NAMESPACE | 1 + R/post_treatment.R | 108 +++++++++++++++++++++++++++++-------- man/tp_mcmc_convergence.Rd | 26 ++++++--- man/tp_parse_mcmc.Rd | 9 +++- man/tp_parse_smc.Rd | 9 +++- 6 files changed, 121 insertions(+), 35 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 1543074..2a7e18d 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -33,7 +33,8 @@ Imports: tibble, dplyr, coda, - crayon + crayon, + tidyr Suggests: devtools, ggplot2, diff --git a/NAMESPACE b/NAMESPACE index 3d2ac0f..088e370 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -22,6 +22,7 @@ export(tp_tempdir) export(tp_treeppl_json) export(tp_write_data) export(tp_write_model) +importFrom("methods", "new") importFrom(bnpsd,tree_reorder) importFrom(methods,is) importFrom(phangorn, diff --git a/R/post_treatment.R b/R/post_treatment.R index 4fb7659..62943d8 100644 --- a/R/post_treatment.R +++ b/R/post_treatment.R @@ -5,12 +5,16 @@ #' The function internally removes sweeps with an undefined normalizing constant. #' #' @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 @@ -42,7 +46,7 @@ #' } #' #' @export -tp_parse_smc <- function(json_path) { +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) @@ -93,14 +97,14 @@ tp_parse_smc <- function(json_path) { 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) { tibble::tibble( particle = particle_id, - samples = as.numeric(unlist(s)) + sample = as.numeric(unlist(s)) ) }) } @@ -130,7 +134,7 @@ tp_parse_smc <- function(json_path) { } # 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, @@ -138,6 +142,25 @@ tp_parse_smc <- function(json_path) { ) |> dplyr::select(-"total_lweight") |> dplyr::relocate("sweep") + + + if (!wide) { + return(result_df) + } else if (wide && !"parameter" %in% colnames(result_df)) { + message("Parameter names could not be found; returning data frame in long format.") + return(result_df) + } else if ("parameter" %in% colnames(result_df) && wide) { + result_df <- result_df |> + dplyr::group_by(sweep, 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 + ) + return(result_df) + } } @@ -147,12 +170,16 @@ tp_parse_smc <- function(json_path) { #' single tidy tibble of samples, one row per iteration. #' #' @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.} #' \item{parameter}{Parameter name, if present in the input JSON.} -#' \item{samples}{Sampled value.} +#' \item{sample}{Sampled value.} #' } #' #' @examples @@ -179,7 +206,7 @@ tp_parse_smc <- function(json_path) { #' } #' #' @export -tp_parse_mcmc <- function(json_path) { +tp_parse_mcmc <- 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) @@ -196,7 +223,7 @@ tp_parse_mcmc <- function(json_path) { run = run_id, iteration = iteration_id, parameter = names(param_values), - samples = as.numeric(unlist(param_values)) + sample = as.numeric(unlist(param_values)) ) }) } else { @@ -204,7 +231,7 @@ tp_parse_mcmc <- function(json_path) { tibble::tibble( run = run_id, iteration = iteration_id, - samples = as.numeric(unlist(s)) + sample = as.numeric(unlist(s)) ) }) } @@ -216,7 +243,21 @@ tp_parse_mcmc <- function(json_path) { stop("All runs failed") } - return(result_df) + # long or wide + if (!wide) { + return(result_df) + } else if (wide && !"parameter" %in% colnames(result_df)) { + message("Parameter names could not be found; returning data frame in long format.") + return(result_df) + } else if ("parameter" %in% colnames(result_df) && wide) { + result_df <- result_df |> + tidyr::pivot_wider( + id_cols = c(run, iteration), + names_from = parameter, + values_from = sample + ) + return(result_df) + } } @@ -432,11 +473,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{ @@ -449,15 +492,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", @@ -466,6 +528,7 @@ tp_mcmc_convergence <- function(treeppl_out) { ) } + # coda::mcmc objects runs <- sort(unique(treeppl_out$run)) parameters <- unique(treeppl_out$parameter) @@ -474,17 +537,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."] @@ -500,7 +563,6 @@ tp_mcmc_convergence <- function(treeppl_out) { if (!has_parameter) { result <- dplyr::select(result, -"parameter") } - return(result) } 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 53f934b..e40912e 100644 --- a/man/tp_parse_mcmc.Rd +++ b/man/tp_parse_mcmc.Rd @@ -4,17 +4,22 @@ \alias{tp_parse_mcmc} \title{Parse TreePPL MCMC output into a tidy data frame} \usage{ -tp_parse_mcmc(json_path) +tp_parse_mcmc(json_path, wide = TRUE) } \arguments{ \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.} \item{parameter}{Parameter name, if present in the input JSON.} -\item{samples}{Sampled value.} +\item{sample}{Sampled value.} } } \description{ diff --git a/man/tp_parse_smc.Rd b/man/tp_parse_smc.Rd index f930bf1..eac4dec 100644 --- a/man/tp_parse_smc.Rd +++ b/man/tp_parse_smc.Rd @@ -4,17 +4,22 @@ \alias{tp_parse_smc} \title{Parse TreePPL SMC output into a tidy data frame} \usage{ -tp_parse_smc(json_path) +tp_parse_smc(json_path, wide = TRUE) } \arguments{ \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 From 2103cd8cbb1d4184ceb9bf16da4e8c82c7a13f83 Mon Sep 17 00:00:00 2001 From: foersterst Date: Tue, 8 Sep 2026 15:37:50 +0200 Subject: [PATCH 04/10] Always wide: arbitrary parameter names when the json file has no parameter names. Applies to tp_parse_smc() and tp_parse_mcmc(). --- R/post_treatment.R | 51 +++++++++++++++++++++------------------------- 1 file changed, 23 insertions(+), 28 deletions(-) diff --git a/R/post_treatment.R b/R/post_treatment.R index 62943d8..e230733 100644 --- a/R/post_treatment.R +++ b/R/post_treatment.R @@ -102,9 +102,11 @@ tp_parse_smc <- function(json_path, wide = TRUE) { }) } else { samples_df <- purrr::imap_dfr(samples, function(s, particle_id) { + values <- as.numeric(unlist(s)) tibble::tibble( particle = particle_id, - sample = as.numeric(unlist(s)) + parameter = paste0("x", seq_along(values)), + sample = values ) }) } @@ -143,24 +145,18 @@ tp_parse_smc <- function(json_path, wide = TRUE) { dplyr::select(-"total_lweight") |> dplyr::relocate("sweep") - if (!wide) { - return(result_df) - } else if (wide && !"parameter" %in% colnames(result_df)) { - message("Parameter names could not be found; returning data frame in long format.") - return(result_df) - } else if ("parameter" %in% colnames(result_df) && wide) { - result_df <- result_df |> - dplyr::group_by(sweep, 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 - ) return(result_df) } + result_df |> + dplyr::group_by(sweep, 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 + ) } @@ -228,10 +224,12 @@ tp_parse_mcmc <- function(json_path, wide = TRUE) { }) } else { purrr::imap_dfr(samples, function(s, iteration_id) { + values <- as.numeric(unlist(s)) tibble::tibble( run = run_id, iteration = iteration_id, - sample = as.numeric(unlist(s)) + parameter = paste0("x", seq_along(values)), + sample = values ) }) } @@ -246,21 +244,18 @@ tp_parse_mcmc <- function(json_path, wide = TRUE) { # long or wide if (!wide) { return(result_df) - } else if (wide && !"parameter" %in% colnames(result_df)) { - message("Parameter names could not be found; returning data frame in long format.") - return(result_df) - } else if ("parameter" %in% colnames(result_df) && wide) { - result_df <- result_df |> - tidyr::pivot_wider( - id_cols = c(run, iteration), - names_from = parameter, - values_from = sample - ) - return(result_df) } + + result_df |> + tidyr::pivot_wider( + id_cols = c(run, iteration), + names_from = parameter, + values_from = sample + ) } + #' Parse TreePPL json output for host repertoire model #' #' @description From 13d64f9153ed2a14e1a39ef31dcdfc187423783e Mon Sep 17 00:00:00 2001 From: foersterst Date: Tue, 8 Sep 2026 16:51:54 +0200 Subject: [PATCH 05/10] fix out names for mcmc with n_runs > 1; added time units to printed summary at the end of tp_run(). --- R/run.R | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/R/run.R b/R/run.R index 46b4b0d..a9a46f1 100644 --- a/R/run.R +++ b/R/run.R @@ -81,7 +81,7 @@ tp_run <- function( 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 # list[[key]] will return the first key @@ -113,7 +113,15 @@ tp_run <- function( 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 { @@ -146,7 +154,7 @@ tp_run <- function( crayon::bold("Analysis Summary\n"), "-----------------------------\n", crayon::green("Status: "), "Completed\n", - crayon::cyan("Time elapsed: "), et, "\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", From e18730c2bed4450bbf533b42d33bf6b330d717f0 Mon Sep 17 00:00:00 2001 From: foersterst Date: Tue, 8 Sep 2026 17:32:04 +0200 Subject: [PATCH 06/10] tp_parse_mcmc() parses JSON files created using --incremental-printing in TreePPL. --- R/post_treatment.R | 34 ++++++++++++++++++++++++---------- 1 file changed, 24 insertions(+), 10 deletions(-) diff --git a/R/post_treatment.R b/R/post_treatment.R index e230733..b5c1121 100644 --- a/R/post_treatment.R +++ b/R/post_treatment.R @@ -203,16 +203,32 @@ tp_parse_smc <- function(json_path, wide = TRUE) { #' #' @export tp_parse_mcmc <- 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_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) + } - # parse mcmc runs - parse_run <- function(run, run_id) { - samples <- run$samples + # 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( @@ -223,6 +239,7 @@ tp_parse_mcmc <- function(json_path, wide = TRUE) { ) }) } 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( @@ -235,13 +252,11 @@ tp_parse_mcmc <- function(json_path, wide = TRUE) { } } - 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") } - # long or wide if (!wide) { return(result_df) } @@ -255,7 +270,6 @@ tp_parse_mcmc <- function(json_path, wide = TRUE) { } - #' Parse TreePPL json output for host repertoire model #' #' @description From fb40dfefdef4b4557357c439c26a2be102b3605c Mon Sep 17 00:00:00 2001 From: foersterst Date: Wed, 23 Sep 2026 17:58:48 +0200 Subject: [PATCH 07/10] added argument verbose to tp_run(); may set it to FALSE in the test suite. --- R/run.R | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/R/run.R b/R/run.R index a9a46f1..c3f78b2 100644 --- a/R/run.R +++ b/R/run.R @@ -61,6 +61,7 @@ tp_run <- function( out_file_name = "out", n_runs = 1, n_processes = 3, + verbose = TRUE, ... ) { # start time @@ -162,7 +163,7 @@ tp_run <- function( ) # print run info summary - cat(run_info) + if (verbose) cat(run_info) # parse JSON to tidy data frames & return # # get model category: this is needed because at the moment, we do not have parsers From bb61fa1f71993d9e462bfd0c5760955bb9ddac92 Mon Sep 17 00:00:00 2001 From: foersterst Date: Wed, 30 Sep 2026 15:21:26 +0200 Subject: [PATCH 08/10] set tp_run(..., verbose = FALSE) in the tests. --- tests/testthat/test-run.R | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/tests/testthat/test-run.R b/tests/testthat/test-run.R index 2f38b42..ccde83c 100644 --- a/tests/testthat/test-run.R +++ b/tests/testthat/test-run.R @@ -9,7 +9,8 @@ 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") + data = tp_data(data_input = "crbd"), + verbose = FALSE ) expect_equal(2, length(unique(run_smc$sweep))) }) @@ -20,7 +21,8 @@ test_that("Test-run_1b : tp_run MCMC", { sampler = tp_compile(model = "crbd", method = "mcmc", iterations = 10), data = tp_data(data_input = "crbd"), n_runs = 2, - n_processes = 2 + n_processes = 2, + verbose = FALSE ) expect_equal(2, length(unique(run_mcmc$run))) }) @@ -30,7 +32,8 @@ test_that("Test-run_1c : tp_run custom_name", { 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" + out_file_name = "test_out", + verbose = FALSE ) expect_equal(2, length(unique(run_smc$sweep))) }) @@ -40,7 +43,8 @@ test_that("Test-run_1d : tp_run threading", { 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 + n_processes = 2, + verbose = FALSE ) expect_equal(2, length(unique(run_smc$sweep))) }) @@ -50,7 +54,8 @@ test_that("Test-run_1e : tp_run no_parser", { 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 + n_processes = 2, + verbose = FALSE ) expect_true(is.character(run_smc)) }) From 42653a217c4fc641c1c1874f3f2e4e0ddcf89857 Mon Sep 17 00:00:00 2001 From: foersterst Date: Fri, 2 Oct 2026 17:39:50 +0200 Subject: [PATCH 09/10] small fixes to resolve warnings and notes during CMD Check. --- NAMESPACE | 2 +- R/post_treatment.R | 16 ++++++++-------- R/run.R | 4 +++- man/tp_run.Rd | 4 ++++ 4 files changed, 16 insertions(+), 10 deletions(-) diff --git a/NAMESPACE b/NAMESPACE index 088e370..aba6fc7 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -22,9 +22,9 @@ export(tp_tempdir) export(tp_treeppl_json) export(tp_write_data) export(tp_write_model) -importFrom("methods", "new") importFrom(bnpsd,tree_reorder) importFrom(methods,is) +importFrom(methods,new) importFrom(phangorn, add_edge_length, allCompat diff --git a/R/post_treatment.R b/R/post_treatment.R index b5c1121..6f3fdf8 100644 --- a/R/post_treatment.R +++ b/R/post_treatment.R @@ -149,13 +149,13 @@ tp_parse_smc <- function(json_path, wide = TRUE) { return(result_df) } result_df |> - dplyr::group_by(sweep, parameter) |> + 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 + id_cols = c("sweep", "particle", "log_weight", "norm_constant", "norm_weight"), + names_from = "parameter", + values_from = "sample" ) } @@ -263,9 +263,9 @@ tp_parse_mcmc <- function(json_path, wide = TRUE) { result_df |> tidyr::pivot_wider( - id_cols = c(run, iteration), - names_from = parameter, - values_from = sample + id_cols = c("run", "iteration"), + names_from = "parameter", + values_from = "sample" ) } @@ -521,7 +521,7 @@ tp_mcmc_convergence <- function(treeppl_out) { if (!"sample" %in% colnames(treeppl_out)) { treeppl_out <- treeppl_out |> tidyr::pivot_longer( - cols = -c(run, iteration), + cols = -c("run", "iteration"), names_to = "parameter", values_to = "sample" ) diff --git a/R/run.R b/R/run.R index c3f78b2..66756f8 100644 --- a/R/run.R +++ b/R/run.R @@ -14,6 +14,8 @@ #' (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 @@ -170,7 +172,7 @@ tp_run <- function( # 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[mc$model_name == mod, ]$category + mod_cat <- mc$category[match(mod, mc$model_name)] # model categories with unavailable parsers no_parsers <- c( diff --git a/man/tp_run.Rd b/man/tp_run.Rd index c046d6f..84436c4 100644 --- a/man/tp_run.Rd +++ b/man/tp_run.Rd @@ -11,6 +11,7 @@ tp_run( out_file_name = "out", n_runs = 1, n_processes = 3, + verbose = TRUE, ... ) } @@ -30,6 +31,9 @@ you want to save the output. Defaults to \code{tp_tempdir()}.} \item{n_processes}{(\code{integer}) the number of parallel processes to use. Cannot be greater than \code{n_runs}.} +\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{ From d1cf89ba1700ef8a2908ee7d37bfea08f4bdab68 Mon Sep 17 00:00:00 2001 From: foersterst Date: Wed, 7 Oct 2026 11:27:23 +0200 Subject: [PATCH 10/10] fix to tp_model_library(): it now works with the TPPLC installation. --- R/utils.R | 16 +--------------- 1 file changed, 1 insertion(+), 15 deletions(-) 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),