From 68b48b7708243b359f8d8f70e478fa2837a47b5c Mon Sep 17 00:00:00 2001 From: foersterst Date: Tue, 25 Aug 2026 16:19:43 +0200 Subject: [PATCH 1/7] fix to tp_parse_smc() to ensure that weights are lined up with samples. Update to the function's documentation. --- DESCRIPTION | 10 +- R/post_treatment.R | 241 ++++++++++++++++++++++++++++---------------- man/tp_parse_smc.Rd | 38 +++++-- 3 files changed, 191 insertions(+), 98 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index fb62dfe..9b4d7ad 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -28,7 +28,10 @@ Imports: curl, bnpsd, rlang, - phangorn + phangorn, + purrr, + tibble, + dplyr Suggests: devtools, ggplot2, @@ -42,12 +45,11 @@ Suggests: treeio (>= 1.32.0), rmarkdown, readr, - pandoc, - dplyr + pandoc VignetteBuilder: knitr Config/Needs/website: rmarkdown Remotes: github::maribraga/evolnets, bioc::ggtree, bioc::treeio -Config/roxygen2/version: 8.0.0 +Config/roxygen2/version: 8.1.0 diff --git a/R/post_treatment.R b/R/post_treatment.R index 03631bc..e18172b 100644 --- a/R/post_treatment.R +++ b/R/post_treatment.R @@ -1,64 +1,130 @@ -#' Parse simple TreePPL json SMC output +#' Parse TreePPL SMC output into a tidy data frame #' -#' @description -#' `tp_parse_smc` takes TreePPL json SMC output and returns a data.frame +#' Converts a list of parsed SMC sweeps (from \code{tp_run()}) 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 character vector giving the TreePPL json output -#' produced by [tp_run] using an SMC method. +#' @param treeppl_out A list of sweeps parsed from a SMC JSON output: i.e., +#' the output object of \code{tp_run()}. #' -#' @details -#' Particles with -Inf weight are removed. +#' @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{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 +#' total log weight across all particles is 1.} +#' } #' +#' @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 +#' ) +#' mod_smc <- tp_run(sampler = sampler_smc, data = path_data) +#' +#' tp_parse_smc(mod_smc) +#' } #' -#' @return A data frame with the output from inference in TreePPL. #' @export tp_parse_smc <- function(treeppl_out) { + parse_sweep <- function(sweep, sweep_id) { + # remove sweeps with nan norm const + if (identical(sweep$normConst, "nan")) { + message("Removing sweep without normalizing constant (sweep ", sweep_id, ")") + return(NULL) + } - result_df <- list() + # sanity check to detect mismatches in the number of samples & weights + if (length(sweep$samples) != length(sweep$weights)) { + stop( + "Sweep ", sweep_id, ": samples (n=", length(sweep$samples), + ") and weights (n=", length(sweep$weights), ") have different lengths." + ) + } - for (i in seq_along(treeppl_out)) { + norm_const <- sweep$normConst + + # convert weight lists to numeric vectors, considering that some weights in the list + # may not be numeric, e.g., {"__float__": "-inf"} + log_weights <- purrr::map_dbl(sweep$weights, function(w) { + if (is.list(w) && !is.null(w[["__float__"]])) { + # use R convention for "nan" and "inf" (i.e., NaN, Inf, -Inf) + val <- tolower(as.character(w[["__float__"]])) + dplyr::case_when( + val == "-inf" ~ -Inf, + val == "inf" ~ Inf, + val == "nan" ~ NaN, + TRUE ~ as.numeric(val) + ) + } else { + as.numeric(w) + } + }) - # remove sweeps with nan norm const - if (treeppl_out[[i]]$normConst == "nan"){ - print("Removing sweep without normalizing constant") - } else { + samples <- sweep$samples - samples_c <- unlist(treeppl_out[[i]]$samples) - log_weight_c <- unlist(treeppl_out[[i]]$weights) + has_parameter_names <- is.list(samples[[1]]) && !is.null(samples[[1]][["__data__"]]) - if(is.null(names(samples_c))){ - result_df <- rbind(result_df, - data.frame(sweep = i, - samples = samples_c, - log_weight = log_weight_c, - norm_constant = treeppl_out[[i]]$normConst) + # if the json has parameter names + if (has_parameter_names) { + samples_df <- purrr::imap_dfr(samples, function(s, particle_id) { + param_values <- s[["__data__"]] + tibble::tibble( + particle = particle_id, + parameter = names(param_values), + samples = as.numeric(unlist(param_values)) ) - } else { - result_df <- rbind(result_df, - data.frame(sweep = i, - parameter = names(samples_c), - samples = samples_c, - log_weight = log_weight_c, - norm_constant = treeppl_out[[i]]$normConst) + }) + } else { + samples_df <- purrr::imap_dfr(samples, function(s, particle_id) { + tibble::tibble( + particle = particle_id, + samples = as.numeric(unlist(s)) ) - } + }) } + + # here we create indexes for the particles so that we can use them to + # match weights with samples using left_join() + weights_df <- tibble::tibble( + particle = seq_along(log_weights), + log_weight = log_weights + ) + + # left_join() ensures the correct alignment of samples with their weights + samples_df |> + dplyr::left_join(weights_df, by = "particle") |> + dplyr::mutate( + sweep = sweep_id, + norm_constant = norm_const + ) |> + dplyr::select(-"particle") } - # check if all sweeps were removed - if (length(result_df) == 0) { + # parse sweeps one at time and rbind the results + result_df <- purrr::imap_dfr(treeppl_out, parse_sweep) + + if (nrow(result_df) == 0) { stop("All sweeps failed") - } else{ - - result_df <- result_df |> - # remove particles with -Inf weight - dplyr::mutate(log_weight = as.numeric(.data$log_weight)) |> - dplyr::filter(!is.infinite(.data$log_weight)) |> - dplyr::mutate(total_lweight = .data$log_weight + .data$norm_constant) |> - dplyr::mutate(norm_weight = exp(.data$total_lweight - max(.data$total_lweight))) |> - dplyr::select(-"total_lweight") } - return(result_df) + + # calculate the normalized weight and return + result_df |> + dplyr::filter(!is.infinite(.data$log_weight)) |> + dplyr::mutate( + total_lweight = .data$log_weight + .data$norm_constant, + norm_weight = exp(.data$total_lweight - max(.data$total_lweight)) + ) |> + dplyr::select(-"total_lweight") |> + dplyr::relocate("sweep") } #' Parse simple TreePPL json MCMC output @@ -72,23 +138,27 @@ tp_parse_smc <- function(treeppl_out) { #' @return A data frame with the output from inference in TreePPL. #' @export tp_parse_mcmc <- function(treeppl_out) { - result_df <- list() for (i in seq_along(treeppl_out)) { - samples_c <- unlist(treeppl_out[[i]]$samples) - if(is.null(names(samples_c))){ - result_df <- rbind(result_df, - data.frame(run = i, - samples = samples_c) + if (is.null(names(samples_c))) { + result_df <- rbind( + result_df, + data.frame( + run = i, + samples = samples_c + ) ) } else { - result_df <- rbind(result_df, - data.frame(run = i, - parameter = names(samples_c), - samples = samples_c) + result_df <- rbind( + result_df, + data.frame( + run = i, + parameter = names(samples_c), + samples = samples_c + ) ) } } @@ -111,7 +181,6 @@ tp_parse_mcmc <- function(treeppl_out) { #' @export tp_parse_host_rep <- function(treeppl_out) { - result_list <- list() for (index in seq_along(treeppl_out)) { @@ -145,7 +214,7 @@ tp_parse_host_rep <- function(treeppl_out) { "child2_index" ) - #for (i in seq_along(output_trppl[1][[1]])) { + # for (i in seq_along(output_trppl[1][[1]])) { for (i in seq_along(output_trppl$samples)) { res <- data.frame(matrix(ncol = nbr_col, nrow = 0)) colnames(res) <- c( @@ -194,18 +263,18 @@ tp_parse_host_rep <- function(treeppl_out) { return(result_list) } -#Recursive function to go deep in the tree +# Recursive function to go deep in the tree peel_tree <- function(subtree, - index, - pindex, - lweight, - lnorm_const, - mu, - beta, - lambda, - prev_age, - start_state, - result) { + index, + pindex, + lweight, + lnorm_const, + mu, + beta, + lambda, + prev_age, + start_state, + result) { base <- c( iteration = as.numeric(index - 1), log_weight = as.numeric(lweight), @@ -227,28 +296,29 @@ peel_tree <- function(subtree, if (!is.null(subtree$left)) { base[["child1_index"]] <- as.numeric(subtree$left$`__data__`$label - 1) - base[["child2_index"]] <- + base[["child2_index"]] <- as.numeric(subtree$right$`__data__`$label - 1) } - base[["end_state"]] <- base[["start_state"]] + base[["end_state"]] <- base[["start_state"]] chang_nbr <- length(subtree$history) if (chang_nbr != 0) { df <- data.frame(matrix(ncol = 2, nrow = chang_nbr)) for (i in 1:chang_nbr) { - #"end_state" + # "end_state" df[i, 1] <- as.numeric(paste(subtree$history[[i]]$`__data__`$repertoire, - collapse = "")) - #"transition_time" + collapse = "" + )) + # "transition_time" df[i, 2] <- as.numeric(subtree$history[[i]]$`__data__`$age) } df <- df[order(-df$X2), ] for (j in 1:chang_nbr) { - base[["start_state"]] <- base[["end_state"]] - base[["end_state"]] <- df[j, 1] - base[["transition_time"]] <- df[j, 2] + base[["start_state"]] <- base[["end_state"]] + base[["end_state"]] <- df[j, 1] + base[["transition_time"]] <- df[j, 2] result[nrow(result) + 1, ] <- base } } else { @@ -295,7 +365,6 @@ peel_tree <- function(subtree, #' @export #' tp_smc_convergence <- function(treeppl_out) { - zs <- treeppl_out |> dplyr::slice_head(n = 1, by = .data$sweep) |> dplyr::pull(.data$norm_constant) @@ -311,11 +380,9 @@ tp_smc_convergence <- function(treeppl_out) { #' @returns Gelman and Rubin's convergence diagnostic #' tp_mcmc_convergence <- function(treeppl_out) { - # create coda::mcmc objects for each run # list(mcmc objects) - #coda::gelman.diag - + # coda::gelman.diag } @@ -328,7 +395,6 @@ tp_mcmc_convergence <- function(treeppl_out) { #' @export #' tp_map_tree <- function(trees_out) { - trees <- trees_out$trees weights <- trees_out$weights @@ -342,9 +408,8 @@ tp_map_tree <- function(trees_out) { # Identify unique topologies trees_ready <- lapply(trees, function(tree) { - # normalize edge lengths for the tip reordering - tree$edge.length <- tree$edge.length/max(tree$edge.length) + tree$edge.length <- tree$edge.length / max(tree$edge.length) # Ladderize to fix edge indices tree_lad <- ladderize_tree(tree) # Order tip labels as similarly as possible @@ -352,14 +417,13 @@ tp_map_tree <- function(trees_out) { # Remove edge lengths to only focus on topology tree_ord$edge.length <- NULL return(tree_ord) - }) # This compresses the list into unique tree topologies unique_topologies <- ape::unique.multiPhylo(trees_ready, use.edge.length = FALSE) # Map every original tree to a unique topology index - #match_indices <- match(trees_ready, unique_topologies) + # match_indices <- match(trees_ready, unique_topologies) match_indices <- attr(unique_topologies, "old.index") # Sum weights for each unique topology @@ -375,20 +439,23 @@ tp_map_tree <- function(trees_out) { # Compute Mean Branch Lengths for the MAP Topology # We take all samples that matched the MAP topology... matching_indices <- which(match_indices == best_index) - matching_trees <- trees[matching_indices] + matching_trees <- trees[matching_indices] matching_weights <- weights[matching_indices] # ...and compute a consensus to average their branch lengths. # Ideally, we should do a weighted average of the lengths, # but ape::consensus uses simple mean. For most purposes, this is sufficient. - final_map <- map <- phangorn::allCompat(matching_trees, rooted=TRUE) |> + final_map <- map <- phangorn::allCompat(matching_trees, rooted = TRUE) |> phangorn::add_edge_length(matching_trees, - fun = function(x) stats::weighted.mean(x, matching_weights)) + fun = function(x) stats::weighted.mean(x, matching_weights) + ) print(paste("MAP Topology found")) print(paste("Posterior Probability:", round(map_prob, 4))) - print(paste("Based on the topology of", length(matching_indices), - "samples out of", length(trees))) + print(paste( + "Based on the topology of", length(matching_indices), + "samples out of", length(trees) + )) return(final_map) } diff --git a/man/tp_parse_smc.Rd b/man/tp_parse_smc.Rd index 6db9ed1..a4003ef 100644 --- a/man/tp_parse_smc.Rd +++ b/man/tp_parse_smc.Rd @@ -2,20 +2,44 @@ % Please edit documentation in R/post_treatment.R \name{tp_parse_smc} \alias{tp_parse_smc} -\title{Parse simple TreePPL json SMC output} +\title{Parse TreePPL SMC output into a tidy data frame} \usage{ tp_parse_smc(treeppl_out) } \arguments{ -\item{treeppl_out}{a character vector giving the TreePPL json output -produced by \link{tp_run} using an SMC method.} +\item{treeppl_out}{A list of sweeps parsed from a SMC JSON output: i.e., +the output object of \code{tp_run()}.} } \value{ -A data frame with the output from inference in TreePPL. +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{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 +total log weight across all particles is 1.} +} } \description{ -\code{tp_parse_smc} takes TreePPL json SMC output and returns a data.frame +Converts a list of parsed SMC sweeps (from \code{tp_run()}) 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 +) +mod_smc <- tp_run(sampler = sampler_smc, data = path_data) + +tp_parse_smc(mod_smc) } -\details{ -Particles with -Inf weight are removed. + } From 40661cd38c22fc4f4db4ed11e4dbb6d0d0ea4b82 Mon Sep 17 00:00:00 2001 From: foersterst Date: Wed, 26 Aug 2026 15:10:00 +0200 Subject: [PATCH 2/7] fix to tp_parse_mcmc() to ensure that values in samples match the utput JSON file(s). Includes updated documentation. --- R/post_treatment.R | 74 ++++++++++++++++++++++++++++++-------------- man/tp_parse_mcmc.Rd | 32 ++++++++++++++++--- 2 files changed, 77 insertions(+), 29 deletions(-) diff --git a/R/post_treatment.R b/R/post_treatment.R index e18172b..fa6d7a9 100644 --- a/R/post_treatment.R +++ b/R/post_treatment.R @@ -127,42 +127,68 @@ tp_parse_smc <- function(treeppl_out) { dplyr::relocate("sweep") } -#' Parse simple TreePPL json MCMC output + +#' Parse TreePPL MCMC output into a tidy data frame #' -#' @description -#' `tp_parse_mcmc` takes TreePPL json MCMC output and returns a data.frame +#' Converts a list of parsed MCMC runs (from \code{tp_run()}) into a +#' single tidy tibble of samples, one row per iteration. #' -#' @param treeppl_out a character vector giving the TreePPL json output -#' produced by [tp_run] using an MCMC method. +#' @param treeppl_out A list of MCMC runs parsed from MCMC JSON output files: i.e., +#' the output object of \code{tp_run()}. +#' +#' @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{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 +#' ) +#' +#' tp_parse_mcmc(mod_mcmc) +#' } #' -#' @return A data frame with the output from inference in TreePPL. #' @export tp_parse_mcmc <- function(treeppl_out) { - result_df <- list() - - for (i in seq_along(treeppl_out)) { - samples_c <- unlist(treeppl_out[[i]]$samples) + parse_run <- function(run, run_id) { + samples <- run$samples + has_parameter_names <- is.list(samples[[1]]) && !is.null(samples[[1]][["__data__"]]) - if (is.null(names(samples_c))) { - result_df <- rbind( - result_df, - data.frame( - run = i, - samples = samples_c + if (has_parameter_names) { + purrr::imap_dfr(samples, function(s, iteration_id) { + param_values <- s[["__data__"]] + tibble::tibble( + run = run_id, + parameter = names(param_values), + samples = as.numeric(unlist(param_values)) ) - ) + }) } else { - result_df <- rbind( - result_df, - data.frame( - run = i, - parameter = names(samples_c), - samples = samples_c + purrr::imap_dfr(samples, function(s, iteration_id) { + tibble::tibble( + run = run_id, + samples = as.numeric(unlist(s)) ) - ) + }) } } + result_df <- purrr::imap_dfr(treeppl_out, parse_run) + + if (nrow(result_df) == 0) { + stop("All runs failed") + } + return(result_df) } diff --git a/man/tp_parse_mcmc.Rd b/man/tp_parse_mcmc.Rd index 1922c9f..7dba5e2 100644 --- a/man/tp_parse_mcmc.Rd +++ b/man/tp_parse_mcmc.Rd @@ -2,17 +2,39 @@ % Please edit documentation in R/post_treatment.R \name{tp_parse_mcmc} \alias{tp_parse_mcmc} -\title{Parse simple TreePPL json MCMC output} +\title{Parse TreePPL MCMC output into a tidy data frame} \usage{ tp_parse_mcmc(treeppl_out) } \arguments{ -\item{treeppl_out}{a character vector giving the TreePPL json output -produced by \link{tp_run} using an MCMC method.} +\item{treeppl_out}{A list of MCMC runs parsed from MCMC JSON output files: i.e., +the output object of \code{tp_run()}.} } \value{ -A data frame with the output from inference in TreePPL. +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{parameter}{Parameter name, if present in the input JSON.} +\item{samples}{Sampled value.} +} } \description{ -\code{tp_parse_mcmc} takes TreePPL json MCMC output and returns a data.frame +Converts a list of parsed MCMC runs (from \code{tp_run()}) 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 +) + +tp_parse_mcmc(mod_mcmc) +} + } From 1fbdb6d713bc2941e891b19be7daf02f754f5761 Mon Sep 17 00:00:00 2001 From: foersterst Date: Thu, 27 Aug 2026 14:00:55 +0200 Subject: [PATCH 3/7] small fix to add column "iteration" to the output of tp_parse_mcmc(); we might need this column for post-processing. --- R/post_treatment.R | 2 ++ 1 file changed, 2 insertions(+) diff --git a/R/post_treatment.R b/R/post_treatment.R index fa6d7a9..24abec8 100644 --- a/R/post_treatment.R +++ b/R/post_treatment.R @@ -169,6 +169,7 @@ tp_parse_mcmc <- function(treeppl_out) { param_values <- s[["__data__"]] tibble::tibble( run = run_id, + iteration = iteration_id, parameter = names(param_values), samples = as.numeric(unlist(param_values)) ) @@ -177,6 +178,7 @@ tp_parse_mcmc <- function(treeppl_out) { purrr::imap_dfr(samples, function(s, iteration_id) { tibble::tibble( run = run_id, + iteration = iteration_id, samples = as.numeric(unlist(s)) ) }) From 72885e7a05800241aa51bdd34b217f2b12a4d66c Mon Sep 17 00:00:00 2001 From: foersterst Date: Mon, 31 Aug 2026 13:28:25 +0200 Subject: [PATCH 4/7] added tp_mcmc_convergence() --- DESCRIPTION | 3 +- R/post_treatment.R | 78 +++++++++++++++++++++++++++++++++++--- man/tp_mcmc_convergence.Rd | 28 ++++++++++++-- man/treepplr-package.Rd | 2 + 4 files changed, 100 insertions(+), 11 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 9b4d7ad..e891573 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -31,7 +31,8 @@ Imports: phangorn, purrr, tibble, - dplyr + dplyr, + coda Suggests: devtools, ggplot2, diff --git a/R/post_treatment.R b/R/post_treatment.R index 24abec8..f824f5e 100644 --- a/R/post_treatment.R +++ b/R/post_treatment.R @@ -401,16 +401,82 @@ tp_smc_convergence <- function(treeppl_out) { } -#' Check for convergence across multiple MCMC runs. +#' Assess MCMC convergence for a TreePPL analysis #' -#' @param treeppl_out a data frame outputted by [tp_parse_mcmc()]. +#' 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. #' -#' @returns Gelman and Rubin's convergence diagnostic +#' @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. #' +#' @return A tibble with one row per parameter, containing: +#' \describe{ +#' \item{parameter}{Parameter name.} +#' \item{ess}{Effective sample size, pooled across all runs.} +#' \item{rhat_upper}{Upper limit of the Gelman-Rubin R-hat statistic. +#' Only computed when `treeppl_out` contains more than one run; +#' otherwise `NA`, with a message explaining why.} +#' } +#' +#' @examples +#' \dontrun{ +#' d <- tp_parse_mcmc(mod_mcmc) +#' tp_mcmc_convergence(d) +#' } +#' +#' @export tp_mcmc_convergence <- function(treeppl_out) { - # create coda::mcmc objects for each run - # list(mcmc objects) - # coda::gelman.diag + has_parameter <- "parameter" %in% names(treeppl_out) + + # stop if the JSON has no parameter names + if (!has_parameter) { + stop( + "Output JSON has no parameter names.\n", + "Re-run the TreePPL model with a named return type (.tppl file) ", + "so that parameters can be identified." + ) + } + + runs <- sort(unique(treeppl_out$run)) + parameters <- unique(treeppl_out$parameter) + + chains_by_parameter <- purrr::map(parameters, function(p) { + chains <- purrr::map(runs, function(r) { + vals <- treeppl_out |> + dplyr::filter(.data$run == r, .data$parameter == p) |> + dplyr::arrange(.data$iteration) |> + dplyr::pull(.data$samples) + coda::mcmc(vals) + }) + coda::mcmc.list(chains) + }) + names(chains_by_parameter) <- parameters + + ess <- purrr::map_dbl(chains_by_parameter, coda::effectiveSize) + + result <- tibble::tibble(parameter = names(ess), ess = ess) + + if (length(runs) > 1) { + rhat_upper <- purrr::map_dbl(chains_by_parameter, function(chain_list) { + coda::gelman.diag(chain_list)$psrf[, "Upper C.I."] + }) + result$rhat_upper <- rhat_upper[result$parameter] + } else { + result$rhat_upper <- NA + message( + "Only one run detected; Gelman-Rubin R-hat requires >= 2 runs and was not computed." + ) + } + + 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 9ac120a..1ec9210 100644 --- a/man/tp_mcmc_convergence.Rd +++ b/man/tp_mcmc_convergence.Rd @@ -2,16 +2,36 @@ % Please edit documentation in R/post_treatment.R \name{tp_mcmc_convergence} \alias{tp_mcmc_convergence} -\title{Check for convergence across multiple MCMC runs.} +\title{Assess MCMC convergence for a TreePPL analysis} \usage{ tp_mcmc_convergence(treeppl_out) } \arguments{ -\item{treeppl_out}{a data frame outputted by \code{\link[=tp_parse_mcmc]{tp_parse_mcmc()}}.} +\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.} } \value{ -Gelman and Rubin's convergence diagnostic +A tibble with one row per parameter, containing: +\describe{ +\item{parameter}{Parameter name.} +\item{ess}{Effective sample size, pooled across all runs.} +\item{rhat_upper}{Upper limit of the Gelman-Rubin R-hat statistic. +Only computed when \code{treeppl_out} contains more than one run; +otherwise \code{NA}, with a message explaining why.} +} } \description{ -Check for convergence across multiple MCMC runs. +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. +} +\examples{ +\dontrun{ +d <- tp_parse_mcmc(mod_mcmc) +tp_mcmc_convergence(d) +} + } diff --git a/man/treepplr-package.Rd b/man/treepplr-package.Rd index 108d026..4306e7c 100644 --- a/man/treepplr-package.Rd +++ b/man/treepplr-package.Rd @@ -23,6 +23,8 @@ Useful links: 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 } } From 91a220bdab1870537d1abdae4986486d285537cc Mon Sep 17 00:00:00 2001 From: foersterst Date: Tue, 1 Sep 2026 16:04:41 +0200 Subject: [PATCH 5/7] 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 b9beb43aed7e2cd2b464bd7a38e1652ecdbd6f0a Mon Sep 17 00:00:00 2001 From: foersterst Date: Tue, 1 Sep 2026 18:28:52 +0200 Subject: [PATCH 6/7] 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 f402742135d22a5d096499afba858cffdfd995b5 Mon Sep 17 00:00:00 2001 From: foersterst Date: Wed, 2 Sep 2026 15:53:54 +0200 Subject: [PATCH 7/7] 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