From 02e2cd5ebcaa8e236b88780db26ad240af13d6ce Mon Sep 17 00:00:00 2001 From: njtierney Date: Tue, 11 Mar 2025 16:26:00 +1100 Subject: [PATCH 01/55] moving tests and utils over from add-snaper-hmc branch Merge branch 'add-snaper-hmc' into adaptive-hmc-v2-i765 # Conflicts: # R/inference_class.R # tests/testthat/test_posteriors_geweke.R --- R/utils.R | 82 +++++++++++++++ tests/testthat/helpers.R | 2 + tests/testthat/test-adaptive-hmc.R | 61 +++++++++++ tests/testthat/test_inference.R | 3 + tests/testthat/test_posteriors_binomial.R | 32 ++++++ tests/testthat/test_posteriors_chi_squared.R | 22 ++++ tests/testthat/test_posteriors_geweke.R | 100 ++++++++++++------- 7 files changed, 266 insertions(+), 36 deletions(-) create mode 100644 tests/testthat/test-adaptive-hmc.R diff --git a/R/utils.R b/R/utils.R index 679a479c..4ff9601a 100644 --- a/R/utils.R +++ b/R/utils.R @@ -1079,3 +1079,85 @@ n_warmup <- function(x) { x_info <- attr(x, "model_info") x_info$warmup } + +as_tensor_spec <- function(tensor) { + tf$TensorSpec(shape = tensor$shape, dtype = tensor$dtype) +} + +maybe_make_tensor_shape <- function(x) { + if (tf$is_tensor(x)) { + as_tensor_spec(x) + } else { + x + } +} + +# get the final model parameter state from a chain as returned in the all_states +# object from tfp$mcmc$sample_chain +get_last_state <- function(all_states) { + n_iter <- dim(all_states)[1] + tf$gather(all_states, n_iter - 1L, 0L) +} + +# find out if MCMC steps had non-finite acceptance probabilities +bad_steps <- function(kernel_results) { + log_accept_ratios <- recursive_get_log_accept_ratio(kernel_results) + !is.finite(log_accept_ratios) +} + + +# recursively extract the log accaptance ratio from the MCMC kernel +recursive_get_log_accept_ratio <- function(kernel_results) { + nm <- names(kernel_results) + if ("log_accept_ratio" %in% nm) { + log_accept_ratios <- kernel_results$log_accept_ratio + } else if ("inner_results" %in% nm) { + log_accept_ratios <- recursive_get_log_accept_ratio( + kernel_results$inner_results + ) + } else { + stop("non-standard kernel structure") + } + as.array(log_accept_ratios) +} + + +# given a warmed up sampler object, return a compiled TF function that generates +# a new burst of samples from samples from it +make_sampler_function <- function(warm_sampler) { + # make the uncompiled function (with curried arguments) + sample_raw <- function(current_state, n_samples) { + results <- tfp$mcmc$sample_chain( + # how many iterations + num_results = n_samples, + # where to start from + current_state = current_state, + # kernel + kernel = warm_sampler$kernel, + # tuned sampler settings + previous_kernel_results = warm_sampler$kernel_results, + # what to trace (nothing) + trace_fn = function(current_state, kernel_results) { + # could compute badness here to save memory? + # is.finite(kernel_results$inner_results$inner_results$inner_results$log_accept_ratio) + kernel_results + } + ) + # return the parameter states and the kernel results + list( + all_states = results$all_states, + kernel_results = results$trace + ) + } + + # compile it into a concrete function and return + sample <- tf_function( + sample_raw, + list( + as_tensorspec(warm_sampler$current_state), + tf$TensorSpec(shape = c(), dtype = tf$int32) + ) + ) + + sample +} diff --git a/tests/testthat/helpers.R b/tests/testthat/helpers.R index fe589749..5273e3d8 100644 --- a/tests/testthat/helpers.R +++ b/tests/testthat/helpers.R @@ -714,6 +714,7 @@ p_theta_greta <- function( data, p_theta, p_x_bar_theta, + # TODO note that we might want to change this to adaptive_hmc() sampler = hmc(), warmup = 1000 ) { @@ -934,6 +935,7 @@ get_distribution_name <- function(x) { check_samples <- function( x, iid_function, + # TODO note that we might want to change this to adaptive_hmc sampler = hmc(), n_effective = 3000, title = NULL, diff --git a/tests/testthat/test-adaptive-hmc.R b/tests/testthat/test-adaptive-hmc.R new file mode 100644 index 00000000..45874fa6 --- /dev/null +++ b/tests/testthat/test-adaptive-hmc.R @@ -0,0 +1,61 @@ +set.seed(2025 - 02 - 13) + +test_that("bad mcmc proposals are rejected", { + skip_if_not(check_tf_version()) + + # set up for numerical rejection of initial location + x <- rnorm(10000, 1e60, 1) + z <- normal(-1e60, 1e-60) + distribution(x) <- normal(z, 1e-60) + m <- model(z, precision = "single") + + # # catch badness in the progress bar + out <- get_output( + mcmc(m, n_samples = 10, warmup = 0, pb_update = 10) + ) + expect_match(out, "100% bad") + + expect_snapshot( + error = TRUE, + draws <- mcmc( + m, + chains = 1, + n_samples = 2, + warmup = 0, + verbose = FALSE, + sampler = adaptive_hmc(), + initial_values = initials(z = 1e120) + ) + ) + + # really bad proposals + x <- rnorm(100000, 1e120, 1) + z <- normal(-1e120, 1e-120) + distribution(x) <- normal(z, 1e-120) + m <- model(z, precision = "single") + expect_snapshot( + error = TRUE, + mcmc(m, + chains = 1, + n_samples = 1, + warmup = 0, + sampler = adaptive_hmc(), + verbose = FALSE) + ) + + # proposals that are fine, but rejected anyway + z <- normal(0, 1) + m <- model(z, precision = "single") + expect_ok(mcmc( + m, + adaptive_hmc( + epsilon = 100, + # Lmin = 1, + # Lmax = 1 + ), + chains = 1, + n_samples = 5, + warmup = 0, + verbose = FALSE + )) +}) diff --git a/tests/testthat/test_inference.R b/tests/testthat/test_inference.R index 4f883a04..5776572e 100644 --- a/tests/testthat/test_inference.R +++ b/tests/testthat/test_inference.R @@ -422,6 +422,9 @@ test_that("samplers print informatively", { expect_snapshot( hmc(Lmin = 1) ) + expect_snapshot( + adaptive_hmc(Lmin = 1) + ) # # check print sees changed parameters # out <- capture_output(hmc(Lmin = 1), TRUE) diff --git a/tests/testthat/test_posteriors_binomial.R b/tests/testthat/test_posteriors_binomial.R index 92d4a087..5772e7d7 100644 --- a/tests/testthat/test_posteriors_binomial.R +++ b/tests/testthat/test_posteriors_binomial.R @@ -29,3 +29,35 @@ test_that("posterior is correct (binomial)", { suppressWarnings(test <- ks.test(samples, comparison)) expect_gte(test$p.value, 0.01) }) + +test_that("posterior is correct (binomial) with adaptive_hmc", { + skip_if_not(check_tf_version()) + skip_on_os("windows") + # analytic solution to the posterior of the paramter of a binomial + # distribution, with uniform prior + n <- 100 + pos <- rbinom(1, n, runif(1)) + theta <- uniform(0, 1) + distribution(pos) <- binomial(n, theta) + m <- model(theta) + + draws <- get_enough_draws(m, adaptive_hmc(), 2000, verbose = FALSE) + + samples <- as.matrix(draws) + + # analytic solution to posterior is beta(1 + pos, 1 + N - pos) + shape1 <- 1 + pos + shape2 <- 1 + n - pos + + # qq plot against true quantiles + quants <- (1:99) / 100 + q_target <- qbeta(quants, shape1, shape2) + q_est <- quantile(samples, quants) + plot(q_target ~ q_est, main = "binomial posterior") + abline(0, 1) + + n_draws <- round(coda::effectiveSize(draws)) + comparison <- rbeta(n_draws, shape1, shape2) + suppressWarnings(test <- ks.test(samples, comparison)) + expect_gte(test$p.value, 0.01) +}) diff --git a/tests/testthat/test_posteriors_chi_squared.R b/tests/testthat/test_posteriors_chi_squared.R index 1332d347..63d9e842 100644 --- a/tests/testthat/test_posteriors_chi_squared.R +++ b/tests/testthat/test_posteriors_chi_squared.R @@ -19,3 +19,25 @@ test_that("samplers are unbiased for chi-squared", { expect_gte(stat$p.value, 0.01) }) + +test_that("samplers are unbiased for chi-squared", { + skip_if_not(check_tf_version()) + skip_on_os("windows") + df <- 5 + x <- chi_squared(df) + iid <- function(n) rchisq(n, df) + + chi_squared_checked <- check_samples( + x = x, + iid_function = iid, + sampler = adaptive_hmc() + ) + + # do the plotting + qqplot_checked_samples(chi_squared_checked) + + # do a formal hypothesis test + stat <- ks_test_mcmc_vs_iid(chi_squared_checked) + + expect_gte(stat$p.value, 0.01) +}) diff --git a/tests/testthat/test_posteriors_geweke.R b/tests/testthat/test_posteriors_geweke.R index 40a6cac5..2f2148ea 100644 --- a/tests/testthat/test_posteriors_geweke.R +++ b/tests/testthat/test_posteriors_geweke.R @@ -1,38 +1,35 @@ Sys.setenv("RELEASE_CANDIDATE" = "false") -test_that("samplers pass geweke tests", { +# run geweke tests on this model: +# theta ~ normal(mu1, sd1) +# x[i] ~ normal(theta, sd2) +# for i in N + +n <- 10 +mu1 <- rnorm(1, 0, 3) +sd1 <- rlnorm(1) +sd2 <- rlnorm(1) + +# prior (n draws) +p_theta <- function(n) { + rnorm(n, mu1, sd1) +} + +# likelihood +p_x_bar_theta <- function(theta) { + rnorm(n, theta, sd2) +} + +# define the greta model (single precision for slice sampler) +x <- as_data(rep(0, n)) +greta_theta <- normal(mu1, sd1) +distribution(x) <- normal(greta_theta, sd2) +model <- model(greta_theta, precision = "single") + +test_that("hmc sampler passes geweke tests", { skip_if_not(check_tf_version()) - skip_if_not_release() - # nolint start - # run geweke tests on this model: - # theta ~ normal(mu1, sd1) - # x[i] ~ normal(theta, sd2) - # for i in N - # nolint end - - n <- 10 - mu1 <- rnorm(1, 0, 3) - sd1 <- rlnorm(1) - sd2 <- rlnorm(1) - - # prior (n draws) - p_theta <- function(n) { - rnorm(n, mu1, sd1) - } - - # likelihood - p_x_bar_theta <- function(theta) { - rnorm(n, theta, sd2) - } - - # define the greta model (single precision for slice sampler) - x <- as_data(rep(0, n)) - greta_theta <- normal(mu1, sd1) - distribution(x) <- normal(greta_theta, sd2) - model <- model(greta_theta, precision = "single") - # run tests on all available samplers geweke_hmc <- check_geweke( sampler = hmc(), @@ -48,8 +45,13 @@ test_that("samplers pass geweke tests", { geweke_stat_hmc <- geweke_ks(geweke_hmc) testthat::expect_gte(geweke_stat_hmc$p.value, 0.005) +}) - geweke_hmc_rwmh <- check_geweke( +test_that("rwmh sampler passes geweke tests", { + skip_if_not(check_tf_version()) + skip_if_not_release() + + geweke_rwmh <- check_geweke( sampler = rwmh(), model = model, data = x, @@ -59,13 +61,19 @@ test_that("samplers pass geweke tests", { thin = 5 ) - geweke_qq(geweke_hmc_rwmh, title = "RWMH Geweke test") + geweke_qq(geweke_rwmh, title = "RWMH Geweke test") - geweke_stat_rwmh <- geweke_ks(geweke_hmc_rwmh) + geweke_stat_rwmh <- geweke_ks(geweke_rwmh) testthat::expect_gte(geweke_stat_rwmh$p.value, 0.005) +}) - geweke_hmc_slice <- check_geweke( +test_that("slice sampler passes geweke tests", { + skip_if_not(check_tf_version()) + + skip_if_not_release() + + geweke_slice <- check_geweke( sampler = slice(), model = model, data = x, @@ -73,7 +81,27 @@ test_that("samplers pass geweke tests", { p_x_bar_theta = p_x_bar_theta ) - geweke_qq(geweke_hmc_slice, title = "slice sampler Geweke test") + geweke_qq(geweke_slice, title = "slice sampler Geweke test") + + geweke_stat_slice <- geweke_ks(geweke_slice) + + testthat::expect_gte(geweke_stat_slice$p.value, 0.005) +}) + +test_that("adaptive hmc sampler passes geweke tests", { + skip_if_not(check_tf_version()) + + skip_if_not_release() + + geweke_adaptive_hmc <- check_geweke( + sampler = adaptive_hmc(), + model = model, + data = x, + p_theta = p_theta, + p_x_bar_theta = p_x_bar_theta + ) + + geweke_qq(geweke_adaptive_hmc, title = "adaptive hmc sampler Geweke test") - testthat::expect_gte(geweke_hmc_slice$p.value, 0.005) + testthat::expect_gte(geweke_adaptive_hmc$p.value, 0.005) }) From 40c5418e1b55bb2ae9c86299c3ac07e505c9a396 Mon Sep 17 00:00:00 2001 From: njtierney Date: Tue, 11 Mar 2025 16:34:41 +1100 Subject: [PATCH 02/55] add adaptive hmc sampler --- R/samplers.R | 66 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 66 insertions(+) diff --git a/R/samplers.R b/R/samplers.R index b1591e82..f5771953 100644 --- a/R/samplers.R +++ b/R/samplers.R @@ -309,3 +309,69 @@ slice_sampler <- R6Class( } ) ) + +adaptive_hmc_sampler <- R6Class( + "adaptive_hmc_sampler", + inherit = sampler, + public = list( + parameters = list( + # Lmin = 10, + # Lmax = 20, + max_leapfrog_steps = 1000, + # TODO clean up these parameter usage else where + # epsilon = 0.005, + # diag_sd = 1, + # TODO some kind of validity check of method? Currently this can only be + # "SNAPER". + method = "SNAPER" + ), + accept_target = 0.651, + + define_tf_kernel = function(sampler_param_vec) { + dag <- self$model$dag + tfe <- dag$tf_environment + + free_state_size <- length(sampler_param_vec) - 2 + + adaptive_hmc_max_leapfrog_steps <- tf$cast( + x = sampler_param_vec[0], + dtype = tf$int32 + ) + # TODO pipe that in properly + n_warmup <- sampler_param_vec[1] + # adaptive_hmc_epsilon <- sampler_param_vec[1] + # adaptive_hmc_diag_sd <- sampler_param_vec[2:(1+free_state_size)] + + kernel_base <- tfp$experimental$mcmc$SNAPERHamiltonianMonteCarlo( + target_log_prob_fn = dag$tf_log_prob_function_adjusted, + step_size = 1, + num_adaptation_steps = as.integer(self$warmup), + max_leapfrog_steps = adaptive_hmc_max_leapfrog_steps + ) + + sampler_kernel <- tfp$mcmc$DualAveragingStepSizeAdaptation( + inner_kernel = kernel_base, + num_adaptation_steps = as.integer(self$warmup) + ) + + return( + sampler_kernel + ) + }, + sampler_parameter_values = function() { + # random number of integration steps + max_leapfrog_steps <- self$parameters$max_leapfrog_steps + epsilon <- self$parameters$epsilon + diag_sd <- matrix(self$parameters$diag_sd) + method <- self$parameters$method + + # return named list for replacing tensors + list( + adaptive_hmc_max_leapfrog_steps = max_leapfrog_steps, + # adaptive_hmc_epsilon = epsilon, + # adaptive_hmc_diag_sd = diag_sd, + method = method + ) + } + ) +) From 2016dc4745dd2e1251136f7e5ea9344cfe05891b Mon Sep 17 00:00:00 2001 From: njtierney Date: Tue, 11 Mar 2025 16:57:51 +1100 Subject: [PATCH 03/55] brief scratchings of implementing adaptive HMC --- R/samplers.R | 196 ++++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 195 insertions(+), 1 deletion(-) diff --git a/R/samplers.R b/R/samplers.R index f5771953..b8d2acb6 100644 --- a/R/samplers.R +++ b/R/samplers.R @@ -346,6 +346,7 @@ adaptive_hmc_sampler <- R6Class( target_log_prob_fn = dag$tf_log_prob_function_adjusted, step_size = 1, num_adaptation_steps = as.integer(self$warmup), + # TODO potentially remove this line max_leapfrog_steps = adaptive_hmc_max_leapfrog_steps ) @@ -358,6 +359,45 @@ adaptive_hmc_sampler <- R6Class( sampler_kernel ) }, + + # given MCMC kernel `kernel` and initial model parameter state `init`, adapt + # the kernel tuning parameters whilst simultaneously burning-in the model + # parameter state. Return both finalised kernel tuning parameters and the + # burned-in model parameter state + warm_up_sampler = function(kernel, init) { + + # get the predetermined adaptation period of the kernel + n_adapt <- kernel$num_adaptation_steps + + # make the uncompiled function (with curried arguments) + warmup_raw <- function() { + tfp$mcmc$sample_chain( + num_results = n_adapt, + current_state = init, + kernel = kernel, + return_final_kernel_results = TRUE, + trace_fn = function(current_state, kernel_results) { + kernel_results$step #kernel_results + } + ) + } + + # compile it into a concrete function + warmup <- tf_function(warmup_raw) + + # execute it + result <- warmup() + + # return the last (burned-in) state of the model parameters and the final + # (tuned) kernel parameters + list( + kernel = kernel, + kernel_results = result$final_kernel_results, + current_state = get_last_state(result$all_states) + ) + + } + sampler_parameter_values = function() { # random number of integration steps max_leapfrog_steps <- self$parameters$max_leapfrog_steps @@ -372,6 +412,160 @@ adaptive_hmc_sampler <- R6Class( # adaptive_hmc_diag_sd = diag_sd, method = method ) + }, + + # given a warmed up sampler object, return a compiled TF function + # that generates a new burst of samples from samples from it + make_sampler_function = function(warm_sampler) { + + # make the uncompiled function (with curried arguments) + sample_raw <- function(current_state, n_samples) { + results <- tfp$mcmc$sample_chain( + # how many iterations + num_results = n_samples, + # where to start from + current_state = current_state, + # kernel + kernel = warm_sampler$kernel, + # tuned sampler settings + previous_kernel_results = warm_sampler$kernel_results, + # what to trace (nothing) + trace_fn = function(current_state, kernel_results) { + # could compute badness here to save memory? + # is.finite(kernel_results$inner_results$inner_results$inner_results$log_accept_ratio) + kernel_results + } + ) + # return the parameter states and the kernel results + list( + all_states = results$all_states, + kernel_results = results$trace + ) } - ) + + # compile it into a concrete function and return + sample <- tf_function(sample_raw, + list( + as_tensorspec(warm_sampler$current_state), + tf$TensorSpec(shape = c(), + dtype = tf$int32) + )) + + sample + + }, + + run_warmup = function( + n_samples, + pb_update, + ideal_burst_size, + verbose + ) { + perform_warmup <- self$warmup > 0 + if (perform_warmup) { + # adapt and warm up + # self$kernel? + # self$init? + self$warm_up_sampler(kernel, init) + } + + }, + + run_sampling = function( + n_samples, + pb_update, + ideal_burst_size, + trace_batch_size, + thin, + verbose + ) { + perform_sampling <- n_samples > 0 + if (perform_sampling) { + # on exiting during the main sampling period (even if killed by the + # user) trace the free state values + + on.exit(self$trace_values(trace_batch_size), add = TRUE) + + # main sampling + if (verbose) { + pb_sampling <- create_progress_bar( + phase = "sampling", + iter = c(self$warmup, n_samples), + pb_update = pb_update, + width = self$pb_width + ) + iterate_progress_bar( + pb = pb_sampling, + it = 0, + rejects = 0, + chains = self$n_chains, + file = self$pb_file + ) + } else { + pb_sampling <- NULL + } + + ### Adaptive start + print("Sampling parameters") + for (burst in seq_len(n_bursts)) { + burst_result <- sample( + current_state = current_state, + n_samples = burst_size + ) + + # trace the MCMC results from this burst + burst_idx <- (burst - 1) * burst_size + seq_len(burst_size) + trace[burst_idx, , ] <- as.array(burst_result$all_states) + + # overwrite the current state + current_state <- get_last_state(burst_result$all_states) + + # accumulate and report on the badness + new_badness <- sum(bad_steps(burst_result$kernel_results)) + n_bad <- n_bad + new_badness + n_evaluations <- burst * burst_size * n_chains + perc_badness <- round(100 * n_bad / n_evaluations) + + # report on progress + print(sprintf("burst %i of %i (%i%s bad)", + burst, + n_bursts, + perc_badness, + "%")) + + } + ### Adaptive end + + # split up warmup iterations into bursts of sampling + burst_lengths <- self$burst_lengths(n_samples, ideal_burst_size) + completed_iterations <- cumsum(burst_lengths) + + for (burst in seq_along(burst_lengths)) { + # so these bursts are R objects being passed through to python + # and how often to return them + # TF1/2 check todo + # replace with define_tf_draws + self$run_burst(n_samples = burst_lengths[burst], thin = thin) + # trace is it receiving the python + self$trace() + + if (verbose) { + # update the progress bar/percentage log + iterate_progress_bar( + pb = pb_sampling, + it = completed_iterations[burst], + rejects = self$numerical_rejections, + chains = self$n_chains, + file = self$pb_file + ) + + self$write_percentage_log( + total = n_samples, + completed = completed_iterations[burst], + stage = "sampling" + ) + } + } + } # end sampling + }, ) From 435fc07edfde74f04b32d0f7655a2062555421fc Mon Sep 17 00:00:00 2001 From: njtierney Date: Tue, 11 Mar 2025 17:00:46 +1100 Subject: [PATCH 04/55] some more scratchings --- R/samplers.R | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/R/samplers.R b/R/samplers.R index b8d2acb6..629808fd 100644 --- a/R/samplers.R +++ b/R/samplers.R @@ -466,9 +466,10 @@ adaptive_hmc_sampler <- R6Class( # adapt and warm up # self$kernel? # self$init? - self$warm_up_sampler(kernel, init) + result <- self$warm_up_sampler(kernel, init) } + result }, run_sampling = function( From 89efdb4e5fdc7bc3f6504ecb0f42616788594a5f Mon Sep 17 00:00:00 2001 From: njtierney Date: Wed, 12 Mar 2025 14:06:04 +1100 Subject: [PATCH 05/55] style test-adaptive-hmc --- tests/testthat/test-adaptive-hmc.R | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/tests/testthat/test-adaptive-hmc.R b/tests/testthat/test-adaptive-hmc.R index 45874fa6..cac35242 100644 --- a/tests/testthat/test-adaptive-hmc.R +++ b/tests/testthat/test-adaptive-hmc.R @@ -35,12 +35,14 @@ test_that("bad mcmc proposals are rejected", { m <- model(z, precision = "single") expect_snapshot( error = TRUE, - mcmc(m, - chains = 1, - n_samples = 1, - warmup = 0, - sampler = adaptive_hmc(), - verbose = FALSE) + mcmc( + m, + chains = 1, + n_samples = 1, + warmup = 0, + sampler = adaptive_hmc(), + verbose = FALSE + ) ) # proposals that are fine, but rejected anyway From f7746751379254910af4daf8f84a70e870e92710 Mon Sep 17 00:00:00 2001 From: njtierney Date: Wed, 12 Mar 2025 14:06:21 +1100 Subject: [PATCH 06/55] warmup now works for adaptive HMC!! --- R/samplers.R | 349 +++++++++++++++++++++++++++++---------------------- 1 file changed, 197 insertions(+), 152 deletions(-) diff --git a/R/samplers.R b/R/samplers.R index 629808fd..99a70974 100644 --- a/R/samplers.R +++ b/R/samplers.R @@ -94,6 +94,45 @@ slice <- function(max_doublings = 5) { obj } +#' @rdname samplers +#' @export +#' +#' @param epsilon leapfrog stepsize hyperparameter (positive, will be tuned) +#' @param diag_sd estimate of the posterior marginal standard deviations +#' (positive, will be tuned). +#' @param max_leapfrog_steps numeric. Default 1000. Maximum number of leapfrog +#' steps used. The algorithm will determine the optimal number less than this. +#' @param method character length one. Currently can only be "SNAPER" but in +#' the future this may expand to other adaptive samplers. +#' @details For `adaptive_hmc()`. The Lmin and Lmax parameters are learnt and so +#' not provided in this. The number of chains cannot be less than 2, due to +#' how adaptive HMC works. `diag_sd` is used to rescale the parameter space to +#' make it more uniform, and make sampling more efficient. +adaptive_hmc <- function( + max_leapfrog_steps = 1000, + epsilon = 0.1, + diag_sd = 1, + method = "SNAPER" +) { + method <- rlang::arg_match( + arg = method, + values = "SNAPER" + ) + + # nolint end + obj <- list( + parameters = list( + max_leapfrog_steps = max_leapfrog_steps, + epsilon = epsilon, + diag_sd = diag_sd + ), + class = adaptive_hmc_sampler + ) + class(obj) <- c("adaptive_hmc_sampler", "sampler") + obj +} + + #' @noRd #' @export print.sampler <- function(x, ...) { @@ -120,7 +159,6 @@ print.sampler <- function(x, ...) { cat(msg) } - hmc_sampler <- R6Class( "hmc_sampler", inherit = sampler, @@ -331,10 +369,11 @@ adaptive_hmc_sampler <- R6Class( dag <- self$model$dag tfe <- dag$tf_environment - free_state_size <- length(sampler_param_vec) - 2 + # free_state_size <- length(sampler_param_vec) - 2 + free_state_size <- length(sampler_param_vec) adaptive_hmc_max_leapfrog_steps <- tf$cast( - x = sampler_param_vec[0], + x = sampler_param_vec[1], dtype = tf$int32 ) # TODO pipe that in properly @@ -360,21 +399,19 @@ adaptive_hmc_sampler <- R6Class( ) }, - # given MCMC kernel `kernel` and initial model parameter state `init`, adapt - # the kernel tuning parameters whilst simultaneously burning-in the model - # parameter state. Return both finalised kernel tuning parameters and the - # burned-in model parameter state - warm_up_sampler = function(kernel, init) { - + # given MCMC kernel `sampler_kernel` and initial model parameter state + # `free_state`, adapt the kernel tuning parameters whilst simultaneously + # burning-in the model parameter state. Return both finalised kernel + # tuning parameters and the burned-in model parameter state + warm_up_sampler = function(sampler_kernel, n_adapt, free_state) { # get the predetermined adaptation period of the kernel - n_adapt <- kernel$num_adaptation_steps # make the uncompiled function (with curried arguments) warmup_raw <- function() { tfp$mcmc$sample_chain( num_results = n_adapt, - current_state = init, - kernel = kernel, + current_state = free_state, + kernel = sampler_kernel, return_final_kernel_results = TRUE, trace_fn = function(current_state, kernel_results) { kernel_results$step #kernel_results @@ -383,7 +420,7 @@ adaptive_hmc_sampler <- R6Class( } # compile it into a concrete function - warmup <- tf_function(warmup_raw) + warmup <- tensorflow::tf_function(warmup_raw) # execute it result <- warmup() @@ -391,12 +428,11 @@ adaptive_hmc_sampler <- R6Class( # return the last (burned-in) state of the model parameters and the final # (tuned) kernel parameters list( - kernel = kernel, + kernel = sampler_kernel, kernel_results = result$final_kernel_results, current_state = get_last_state(result$all_states) ) - - } + }, sampler_parameter_values = function() { # random number of integration steps @@ -414,159 +450,168 @@ adaptive_hmc_sampler <- R6Class( ) }, - # given a warmed up sampler object, return a compiled TF function - # that generates a new burst of samples from samples from it - make_sampler_function = function(warm_sampler) { - - # make the uncompiled function (with curried arguments) - sample_raw <- function(current_state, n_samples) { - results <- tfp$mcmc$sample_chain( - # how many iterations - num_results = n_samples, - # where to start from - current_state = current_state, - # kernel - kernel = warm_sampler$kernel, - # tuned sampler settings - previous_kernel_results = warm_sampler$kernel_results, - # what to trace (nothing) - trace_fn = function(current_state, kernel_results) { - # could compute badness here to save memory? - # is.finite(kernel_results$inner_results$inner_results$inner_results$log_accept_ratio) - kernel_results - } - ) - # return the parameter states and the kernel results - list( - all_states = results$all_states, - kernel_results = results$trace - ) - } - - # compile it into a concrete function and return - sample <- tf_function(sample_raw, - list( - as_tensorspec(warm_sampler$current_state), - tf$TensorSpec(shape = c(), - dtype = tf$int32) - )) - - sample - - }, - - run_warmup = function( - n_samples, - pb_update, - ideal_burst_size, - verbose - ) { - perform_warmup <- self$warmup > 0 - if (perform_warmup) { - # adapt and warm up - # self$kernel? - # self$init? - result <- self$warm_up_sampler(kernel, init) - } - - result - }, - - run_sampling = function( - n_samples, - pb_update, - ideal_burst_size, - trace_batch_size, - thin, - verbose - ) { - perform_sampling <- n_samples > 0 - if (perform_sampling) { - # on exiting during the main sampling period (even if killed by the - # user) trace the free state values - - on.exit(self$trace_values(trace_batch_size), add = TRUE) - - # main sampling - if (verbose) { - pb_sampling <- create_progress_bar( - phase = "sampling", - iter = c(self$warmup, n_samples), - pb_update = pb_update, - width = self$pb_width + # given a warmed up sampler object, return a compiled TF function + # that generates a new burst of samples from samples from it + make_sampler_function = function(warm_sampler) { + # make the uncompiled function (with curried arguments) + sample_raw <- function(current_state, n_samples) { + results <- tfp$mcmc$sample_chain( + # how many iterations + num_results = n_samples, + # where to start from + current_state = current_state, + # kernel + kernel = warm_sampler$kernel, + # tuned sampler settings + previous_kernel_results = warm_sampler$kernel_results, + # what to trace (nothing) + trace_fn = function(current_state, kernel_results) { + # could compute badness here to save memory? + # is.finite(kernel_results$inner_results$inner_results$inner_results$log_accept_ratio) + kernel_results + } ) - iterate_progress_bar( - pb = pb_sampling, - it = 0, - rejects = 0, - chains = self$n_chains, - file = self$pb_file + # return the parameter states and the kernel results + list( + all_states = results$all_states, + kernel_results = results$trace ) - } else { - pb_sampling <- NULL } - ### Adaptive start - print("Sampling parameters") - for (burst in seq_len(n_bursts)) { - burst_result <- sample( - current_state = current_state, - n_samples = burst_size + # compile it into a concrete function and return + sample <- tensorflow::tf_function( + sample_raw, + list( + as_tensorspec(warm_sampler$current_state), + tf$TensorSpec(shape = c(), dtype = tf$int32) ) + ) - # trace the MCMC results from this burst - burst_idx <- (burst - 1) * burst_size + seq_len(burst_size) - trace[burst_idx, , ] <- as.array(burst_result$all_states) - - # overwrite the current state - current_state <- get_last_state(burst_result$all_states) - - # accumulate and report on the badness - new_badness <- sum(bad_steps(burst_result$kernel_results)) - n_bad <- n_bad + new_badness - n_evaluations <- burst * burst_size * n_chains - perc_badness <- round(100 * n_bad / n_evaluations) - - # report on progress - print(sprintf("burst %i of %i (%i%s bad)", - burst, - n_bursts, - perc_badness, - "%")) + sample + }, + run_warmup = function( + n_samples, + pb_update, + ideal_burst_size, + verbose + ) { + perform_warmup <- self$warmup > 0 + if (perform_warmup) { + # adapt and warm up + param_vec <- unlist(self$sampler_parameter_values()) + sampler_kernel <- self$define_tf_kernel( + sampler_param_vec = param_vec + ) + init <- self$free_state + n_adapt <- as.integer(param_vec) + result <- self$warm_up_sampler( + sampler_kernel = sampler_kernel, + n_adapt = n_adapt, + free_state = init + ) } - ### Adaptive end - # split up warmup iterations into bursts of sampling - burst_lengths <- self$burst_lengths(n_samples, ideal_burst_size) - completed_iterations <- cumsum(burst_lengths) - - for (burst in seq_along(burst_lengths)) { - # so these bursts are R objects being passed through to python - # and how often to return them - # TF1/2 check todo - # replace with define_tf_draws - self$run_burst(n_samples = burst_lengths[burst], thin = thin) - # trace is it receiving the python - self$trace() + result + }, + run_sampling = function( + n_samples, + pb_update, + ideal_burst_size, + trace_batch_size, + thin, + verbose + ) { + perform_sampling <- n_samples > 0 + if (perform_sampling) { + # on exiting during the main sampling period (even if killed by the + # user) trace the free state values + + on.exit(self$trace_values(trace_batch_size), add = TRUE) + + # main sampling if (verbose) { - # update the progress bar/percentage log + pb_sampling <- create_progress_bar( + phase = "sampling", + iter = c(self$warmup, n_samples), + pb_update = pb_update, + width = self$pb_width + ) iterate_progress_bar( pb = pb_sampling, - it = completed_iterations[burst], - rejects = self$numerical_rejections, + it = 0, + rejects = 0, chains = self$n_chains, file = self$pb_file ) + } else { + pb_sampling <- NULL + } - self$write_percentage_log( - total = n_samples, - completed = completed_iterations[burst], - stage = "sampling" + ### Adaptive start + print("Sampling parameters") + for (burst in seq_len(n_bursts)) { + burst_result <- sample( + current_state = current_state, + n_samples = burst_size ) + + # trace the MCMC results from this burst + burst_idx <- (burst - 1) * burst_size + seq_len(burst_size) + trace[burst_idx, , ] <- as.array(burst_result$all_states) + + # overwrite the current state + current_state <- get_last_state(burst_result$all_states) + + # accumulate and report on the badness + new_badness <- sum(bad_steps(burst_result$kernel_results)) + n_bad <- n_bad + new_badness + n_evaluations <- burst * burst_size * n_chains + perc_badness <- round(100 * n_bad / n_evaluations) + + # report on progress + print(sprintf( + "burst %i of %i (%i%s bad)", + burst, + n_bursts, + perc_badness, + "%" + )) } - } - } # end sampling - }, + ### Adaptive end + + # split up warmup iterations into bursts of sampling + burst_lengths <- self$burst_lengths(n_samples, ideal_burst_size) + completed_iterations <- cumsum(burst_lengths) + + for (burst in seq_along(burst_lengths)) { + # so these bursts are R objects being passed through to python + # and how often to return them + # TF1/2 check todo + # replace with define_tf_draws + self$run_burst(n_samples = burst_lengths[burst], thin = thin) + # trace is it receiving the python + self$trace() + + if (verbose) { + # update the progress bar/percentage log + iterate_progress_bar( + pb = pb_sampling, + it = completed_iterations[burst], + rejects = self$numerical_rejections, + chains = self$n_chains, + file = self$pb_file + ) + + self$write_percentage_log( + total = n_samples, + completed = completed_iterations[burst], + stage = "sampling" + ) + } + } + } # end sampling + } + ) ) From 5500bb0614ade127d05754bd5f73df125c9f990b Mon Sep 17 00:00:00 2001 From: njtierney Date: Wed, 12 Mar 2025 14:49:06 +1100 Subject: [PATCH 07/55] sampling now nearly runs! Get this error at the end of iterations: `Error in trace_list_batches[[1]] : subscript out of bounds` - need to investigate further --- R/samplers.R | 55 +++++++++++++++++++++++++++++----------------------- R/utils.R | 4 ++-- 2 files changed, 33 insertions(+), 26 deletions(-) diff --git a/R/samplers.R b/R/samplers.R index 99a70974..5c8bfb85 100644 --- a/R/samplers.R +++ b/R/samplers.R @@ -490,6 +490,9 @@ adaptive_hmc_sampler <- R6Class( sample }, + # prepare a slot to put the warmed up results into + warm_results = NULL, + run_warmup = function( n_samples, pb_update, @@ -512,6 +515,7 @@ adaptive_hmc_sampler <- R6Class( ) } + self$warm_results <- result result }, @@ -550,8 +554,25 @@ adaptive_hmc_sampler <- R6Class( } ### Adaptive start - print("Sampling parameters") - for (burst in seq_len(n_bursts)) { + # split up warmup iterations into bursts of sampling + burst_lengths <- as.integer(self$burst_lengths( + n_samples, + ideal_burst_size + )) + completed_iterations <- cumsum(burst_lengths) + + # use this to compile the warmed version + sample <- make_sampler_function(self$warm_results) + + # TODO + # unsure how to get the first current state? I think this is right? + current_state <- self$warm_results$current_state + trace <- array(NA, dim = c(n_samples, dim(current_state))) + # track numerical rejections + n_bad <- 0 + + for (burst in seq_along(burst_lengths)) { + burst_size <- burst_lengths[burst] burst_result <- sample( current_state = current_state, n_samples = burst_size @@ -571,28 +592,13 @@ adaptive_hmc_sampler <- R6Class( perc_badness <- round(100 * n_bad / n_evaluations) # report on progress - print(sprintf( - "burst %i of %i (%i%s bad)", - burst, - n_bursts, - perc_badness, - "%" - )) - } - ### Adaptive end - - # split up warmup iterations into bursts of sampling - burst_lengths <- self$burst_lengths(n_samples, ideal_burst_size) - completed_iterations <- cumsum(burst_lengths) - - for (burst in seq_along(burst_lengths)) { - # so these bursts are R objects being passed through to python - # and how often to return them - # TF1/2 check todo - # replace with define_tf_draws - self$run_burst(n_samples = burst_lengths[burst], thin = thin) - # trace is it receiving the python - self$trace() + # print(sprintf( + # "burst %i of %i (%i%s bad)", + # burst, + # n_bursts, + # perc_badness, + # "%" + # )) if (verbose) { # update the progress bar/percentage log @@ -612,6 +618,7 @@ adaptive_hmc_sampler <- R6Class( } } } # end sampling + ### Adaptive end } ) ) diff --git a/R/utils.R b/R/utils.R index 4ff9601a..5284724b 100644 --- a/R/utils.R +++ b/R/utils.R @@ -1080,7 +1080,7 @@ n_warmup <- function(x) { x_info$warmup } -as_tensor_spec <- function(tensor) { +as_tensorspec <- function(tensor) { tf$TensorSpec(shape = tensor$shape, dtype = tensor$dtype) } @@ -1151,7 +1151,7 @@ make_sampler_function <- function(warm_sampler) { } # compile it into a concrete function and return - sample <- tf_function( + sample <- tensorflow::tf_function( sample_raw, list( as_tensorspec(warm_sampler$current_state), From 7f8360ee5aa65e2f29a09dcf8c9f8876ebc909c6 Mon Sep 17 00:00:00 2001 From: njtierney Date: Wed, 12 Mar 2025 15:33:40 +1100 Subject: [PATCH 08/55] as a first pass adaptive HMC now works!! Still a bit of cleaning up to do, but we are getting there! --- NAMESPACE | 1 + R/samplers.R | 43 ++++++++++++++++++++++++++++++------------- man/samplers.Rd | 19 +++++++++++++++++++ 3 files changed, 50 insertions(+), 13 deletions(-) diff --git a/NAMESPACE b/NAMESPACE index ede86f97..ce9a7e57 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -157,6 +157,7 @@ export(adagrad) export(adagrad_da) export(adam) export(adamax) +export(adaptive_hmc) export(apply) export(are_null) export(as.greta_model) diff --git a/R/samplers.R b/R/samplers.R index 5c8bfb85..1b868e3e 100644 --- a/R/samplers.R +++ b/R/samplers.R @@ -95,7 +95,6 @@ slice <- function(max_doublings = 5) { } #' @rdname samplers -#' @export #' #' @param epsilon leapfrog stepsize hyperparameter (positive, will be tuned) #' @param diag_sd estimate of the posterior marginal standard deviations @@ -108,6 +107,7 @@ slice <- function(max_doublings = 5) { #' not provided in this. The number of chains cannot be less than 2, due to #' how adaptive HMC works. `diag_sd` is used to rescale the parameter space to #' make it more uniform, and make sampling more efficient. +#' @export adaptive_hmc <- function( max_leapfrog_steps = 1000, epsilon = 0.1, @@ -564,32 +564,49 @@ adaptive_hmc_sampler <- R6Class( # use this to compile the warmed version sample <- make_sampler_function(self$warm_results) - # TODO - # unsure how to get the first current state? I think this is right? current_state <- self$warm_results$current_state - trace <- array(NA, dim = c(n_samples, dim(current_state))) + # trace <- array(NA, dim = c(n_samples, dim(current_state))) # track numerical rejections - n_bad <- 0 + # n_bad <- 0 for (burst in seq_along(burst_lengths)) { burst_size <- burst_lengths[burst] - burst_result <- sample( + batch_results <- sample( current_state = current_state, n_samples = burst_size ) + free_state_draws <- as.array(batch_results$all_states) + + # if there is one sample at a time, and it's rejected, conversion from + # python back to R can drop a dimension, so handle that here. Ugh. + if (n_dim(free_state_draws) != 3) { + dim(free_state_draws) <- c(1, dim(free_state_draws)) + } + + self$last_burst_free_states <- split_chains(free_state_draws) + + n_draws <- nrow(free_state_draws) + if (n_draws > 0) { + free_state <- free_state_draws[n_draws, , , drop = FALSE] + dim(free_state) <- dim(free_state)[-1] + self$free_state <- free_state + } + + self$trace() + # trace the MCMC results from this burst - burst_idx <- (burst - 1) * burst_size + seq_len(burst_size) - trace[burst_idx, , ] <- as.array(burst_result$all_states) + # burst_idx <- (burst - 1) * burst_size + seq_len(burst_size) + # trace[burst_idx, , ] <- as.array(batch_results$all_states) # overwrite the current state - current_state <- get_last_state(burst_result$all_states) + current_state <- get_last_state(batch_results$all_states) # accumulate and report on the badness - new_badness <- sum(bad_steps(burst_result$kernel_results)) - n_bad <- n_bad + new_badness - n_evaluations <- burst * burst_size * n_chains - perc_badness <- round(100 * n_bad / n_evaluations) + # new_badness <- sum(bad_steps(batch_results$kernel_results)) + # n_bad <- n_bad + new_badness + # n_evaluations <- burst * burst_size * self$n_chains + # perc_badness <- round(100 * n_bad / n_evaluations) # report on progress # print(sprintf( diff --git a/man/samplers.Rd b/man/samplers.Rd index 076eded3..e5c825f3 100644 --- a/man/samplers.Rd +++ b/man/samplers.Rd @@ -5,6 +5,7 @@ \alias{hmc} \alias{rwmh} \alias{slice} +\alias{adaptive_hmc} \title{MCMC samplers} \usage{ hmc(Lmin = 5, Lmax = 10, epsilon = 0.1, diag_sd = 1) @@ -12,6 +13,13 @@ hmc(Lmin = 5, Lmax = 10, epsilon = 0.1, diag_sd = 1) rwmh(proposal = c("normal", "uniform"), epsilon = 0.1, diag_sd = 1) slice(max_doublings = 5) + +adaptive_hmc( + max_leapfrog_steps = 1000, + epsilon = 0.1, + diag_sd = 1, + method = "SNAPER" +) } \arguments{ \item{Lmin}{minimum number of leapfrog steps (positive integer, Lmin > Lmax)} @@ -27,6 +35,12 @@ slice(max_doublings = 5) \item{max_doublings}{the maximum number of iterations of the 'doubling' algorithm used to adapt the size of the slice} + +\item{max_leapfrog_steps}{numeric. Default 1000. Maximum number of leapfrog +steps used. The algorithm will determine the optimal number less than this.} + +\item{method}{character length one. Currently can only be "SNAPER" but in +the future this may expand to other adaptive samplers.} } \value{ a \code{sampler} object that can be passed to \code{\link[=mcmc]{mcmc()}}. @@ -53,4 +67,9 @@ is implemented for uniform and normal proposals. \code{slice()} implements a multivariate slice sampling algorithm. The parameter \code{max_doublings} is not tuned during warmup. + +For \code{adaptive_hmc()}. The Lmin and Lmax parameters are learnt and so +not provided in this. The number of chains cannot be less than 2, due to +how adaptive HMC works. \code{diag_sd} is used to rescale the parameter space to +make it more uniform, and make sampling more efficient. } From 1c8f4949131b6b625633b0c7122d288297330a2a Mon Sep 17 00:00:00 2001 From: njtierney Date: Tue, 18 Mar 2025 10:54:52 +1100 Subject: [PATCH 09/55] Use `max_leapfrog_steps` not `Lmin` in a test of adaptive_hmc --- tests/testthat/test_inference.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/testthat/test_inference.R b/tests/testthat/test_inference.R index 5776572e..a76e8bab 100644 --- a/tests/testthat/test_inference.R +++ b/tests/testthat/test_inference.R @@ -423,7 +423,7 @@ test_that("samplers print informatively", { hmc(Lmin = 1) ) expect_snapshot( - adaptive_hmc(Lmin = 1) + adaptive_hmc(max_leapfrog_steps = 1L) ) # # check print sees changed parameters From 015e3d5aa4f93358aa0dda55e20bde32dff604e0 Mon Sep 17 00:00:00 2001 From: njtierney Date: Tue, 18 Mar 2025 11:03:40 +1100 Subject: [PATCH 10/55] update adaptive hmc tests to correctly capture errors - these will not pass currently but this test driven development will be useful. --- tests/testthat/_snaps/adaptive-hmc.md | 41 +++++++++++++++++++++++ tests/testthat/test-adaptive-hmc.R | 47 ++++++++++++++++++++++++--- 2 files changed, 83 insertions(+), 5 deletions(-) create mode 100644 tests/testthat/_snaps/adaptive-hmc.md diff --git a/tests/testthat/_snaps/adaptive-hmc.md b/tests/testthat/_snaps/adaptive-hmc.md new file mode 100644 index 00000000..92c9ab92 --- /dev/null +++ b/tests/testthat/_snaps/adaptive-hmc.md @@ -0,0 +1,41 @@ +# adaptive_hmc() errors when given 1 chain + + Code + draws <- mcmc(m, n_samples = 1, warmup = 1, chains = 1, sampler = adaptive_hmc()) + Condition + Error in `py_call_impl()`: + ! ValueError: SNAPERHMC requires at least 2 chains. Got: 1 + Run `reticulate::py_last_error()` for details. + +# adaptive_hmc() does not error when given 0 warmup + + Code + draws <- mcmc(m, n_samples = 1, warmup = 0, chains = 2, sampler = adaptive_hmc()) + Message + running 2 chains simultaneously on up to 8 CPU cores + Output + + Condition + Error in `self$run_warmup()`: + ! object 'result' not found + +# bad mcmc proposals are rejected + + Code + draws <- mcmc(m, chains = 1, n_samples = 2, warmup = 1, verbose = FALSE, + sampler = adaptive_hmc(), initial_values = initials(z = 1e+120)) + Condition + Error in `self$check_valid_parameters()`: + ! The log density could not be evaluated at these initial values + Try using these initials as the `values` argument in `calculate()` to see what values of subsequent s these initial values lead to. + +--- + + Code + mcmc(m, chains = 1, n_samples = 1, warmup = 1, sampler = adaptive_hmc(), + verbose = FALSE) + Condition + Error in `self$check_reasonable_starting_values()`: + ! Could not find reasonable starting values after 20 attempts. + Please specify initial values manually via the `initial_values` argument + diff --git a/tests/testthat/test-adaptive-hmc.R b/tests/testthat/test-adaptive-hmc.R index cac35242..26c94c75 100644 --- a/tests/testthat/test-adaptive-hmc.R +++ b/tests/testthat/test-adaptive-hmc.R @@ -1,5 +1,37 @@ set.seed(2025 - 02 - 13) +test_that("adaptive_hmc() errors when given 1 chain", { + skip_if_not(check_tf_version()) + x <- normal(0, 1) + m <- model(x) + expect_snapshot( + error = TRUE, + draws <- mcmc( + m, + n_samples = 1, + warmup = 1, + chains = 1, + sampler = adaptive_hmc() + ) + ) +}) + +test_that("adaptive_hmc() does not error when given 0 warmup", { + skip_if_not(check_tf_version()) + x <- normal(0, 1) + m <- model(x) + expect_snapshot( + error = TRUE, + draws <- mcmc( + m, + n_samples = 1, + warmup = 0, + chains = 2, + sampler = adaptive_hmc() + ) + ) +}) + test_that("bad mcmc proposals are rejected", { skip_if_not(check_tf_version()) @@ -11,15 +43,22 @@ test_that("bad mcmc proposals are rejected", { # # catch badness in the progress bar out <- get_output( - mcmc(m, n_samples = 10, warmup = 0, pb_update = 10) + mcmc( + m, + n_samples = 10, + warmup = 0, + pb_update = 10, + sampler = adaptive_hmc() + ) ) + expect_match(out, "100% bad") expect_snapshot( error = TRUE, draws <- mcmc( m, - chains = 1, + chains = 2, n_samples = 2, warmup = 0, verbose = FALSE, @@ -52,10 +91,8 @@ test_that("bad mcmc proposals are rejected", { m, adaptive_hmc( epsilon = 100, - # Lmin = 1, - # Lmax = 1 ), - chains = 1, + chains = 2, n_samples = 5, warmup = 0, verbose = FALSE From a6b13d6f4a6677ab11848c1222935a432b447706 Mon Sep 17 00:00:00 2001 From: njtierney Date: Tue, 18 Mar 2025 16:37:17 +1100 Subject: [PATCH 11/55] Tidying up adaptive_hmc some: - remove unused n_adapt parameter in warm_up_sampler, instead refer to parameter set by user in adaptive_hmc() - removed unrequired sampler_param_vec in define_tf_kernel() --- R/samplers.R | 60 ++++++++++++++++++++++++++++------------------------ 1 file changed, 32 insertions(+), 28 deletions(-) diff --git a/R/samplers.R b/R/samplers.R index 1b868e3e..579dd02f 100644 --- a/R/samplers.R +++ b/R/samplers.R @@ -177,8 +177,9 @@ hmc_sampler <- R6Class( free_state_size <- length(sampler_param_vec) - 2 - # TF1/2 check - # this will likely get replaced... + # TODO we will wrap these up into adaptive sampling to learn better + # similar to the adaptive_hmc - to put more things into TFP, + # rather than going back and forth between R and TF hmc_l <- sampler_param_vec[0] hmc_epsilon <- sampler_param_vec[1] @@ -198,6 +199,12 @@ hmc_sampler <- R6Class( # build the kernel # nolint start + # TODO HMC with momentum would get adaptation addition + # one of the is dual averaging step size adapation (see adaptive_hmc) + # the other is adapt_momentum + # diagonal mass + # https://www.tensorflow.org/probability/api_docs/python/tfp/experimental/mcmc/DiagonalMassMatrixAdaptation + # sampler_kernel <- tfp$mcmc$HamiltonianMonteCarlo( target_log_prob_fn = dag$tf_log_prob_function_adjusted, @@ -365,30 +372,26 @@ adaptive_hmc_sampler <- R6Class( ), accept_target = 0.651, + # TODO + # Eventually we will adapt the other samplers to have an interface + # more in line with this - where we do the adaptation in tensorflow define_tf_kernel = function(sampler_param_vec) { dag <- self$model$dag tfe <- dag$tf_environment - # free_state_size <- length(sampler_param_vec) - 2 - free_state_size <- length(sampler_param_vec) - adaptive_hmc_max_leapfrog_steps <- tf$cast( - x = sampler_param_vec[1], + x = self$parameters$max_leapfrog_steps, dtype = tf$int32 ) - # TODO pipe that in properly - n_warmup <- sampler_param_vec[1] - # adaptive_hmc_epsilon <- sampler_param_vec[1] - # adaptive_hmc_diag_sd <- sampler_param_vec[2:(1+free_state_size)] kernel_base <- tfp$experimental$mcmc$SNAPERHamiltonianMonteCarlo( target_log_prob_fn = dag$tf_log_prob_function_adjusted, step_size = 1, num_adaptation_steps = as.integer(self$warmup), - # TODO potentially remove this line max_leapfrog_steps = adaptive_hmc_max_leapfrog_steps ) + # learns the step size sampler_kernel <- tfp$mcmc$DualAveragingStepSizeAdaptation( inner_kernel = kernel_base, num_adaptation_steps = as.integer(self$warmup) @@ -403,13 +406,13 @@ adaptive_hmc_sampler <- R6Class( # `free_state`, adapt the kernel tuning parameters whilst simultaneously # burning-in the model parameter state. Return both finalised kernel # tuning parameters and the burned-in model parameter state - warm_up_sampler = function(sampler_kernel, n_adapt, free_state) { + warm_up_sampler = function(sampler_kernel, free_state) { # get the predetermined adaptation period of the kernel # make the uncompiled function (with curried arguments) warmup_raw <- function() { tfp$mcmc$sample_chain( - num_results = n_adapt, + num_results = as.integer(sampler_kernel$num_adaptation_steps), current_state = free_state, kernel = sampler_kernel, return_final_kernel_results = TRUE, @@ -499,21 +502,12 @@ adaptive_hmc_sampler <- R6Class( ideal_burst_size, verbose ) { - perform_warmup <- self$warmup > 0 - if (perform_warmup) { - # adapt and warm up - param_vec <- unlist(self$sampler_parameter_values()) - sampler_kernel <- self$define_tf_kernel( - sampler_param_vec = param_vec - ) - init <- self$free_state - n_adapt <- as.integer(param_vec) - result <- self$warm_up_sampler( - sampler_kernel = sampler_kernel, - n_adapt = n_adapt, - free_state = init - ) - } + sampler_kernel <- self$define_tf_kernel() + init <- self$free_state + result <- self$warm_up_sampler( + sampler_kernel = sampler_kernel, + free_state = init + ) self$warm_results <- result result @@ -527,6 +521,16 @@ adaptive_hmc_sampler <- R6Class( thin, verbose ) { + # There's a couple of possibilities + ## warmup already happened previously, so get already-warmup sampler + ## warmup never happened, so get an unwarmed up sampler + + # warmup_status <- self$warmup_status() + # self$warm_results + # self$warmup + # warmup_has_happened <- self$has_warmup_happened() + # warmup + perform_sampling <- n_samples > 0 if (perform_sampling) { # on exiting during the main sampling period (even if killed by the From df52a50971cc606196679f543f4efcb2ae109573 Mon Sep 17 00:00:00 2001 From: njtierney Date: Wed, 19 Mar 2025 13:27:39 +1100 Subject: [PATCH 12/55] Update tests for adaptive_hmc(): - Pass chains argument through to `check_geweke` as adaptive_hmc() requires at least 2 chains - use expect_no_error to test that adaptive_hmc() works with 0 warmup - update snapshots --- tests/testthat/_snaps/adaptive-hmc.md | 18 ++++-------------- tests/testthat/_snaps/inference.md | 8 ++++++++ tests/testthat/helpers.R | 11 +++++++---- tests/testthat/test-adaptive-hmc.R | 3 +-- tests/testthat/test_posteriors_geweke.R | 3 ++- 5 files changed, 22 insertions(+), 21 deletions(-) diff --git a/tests/testthat/_snaps/adaptive-hmc.md b/tests/testthat/_snaps/adaptive-hmc.md index 92c9ab92..f8a67ab7 100644 --- a/tests/testthat/_snaps/adaptive-hmc.md +++ b/tests/testthat/_snaps/adaptive-hmc.md @@ -7,23 +7,13 @@ ! ValueError: SNAPERHMC requires at least 2 chains. Got: 1 Run `reticulate::py_last_error()` for details. -# adaptive_hmc() does not error when given 0 warmup - - Code - draws <- mcmc(m, n_samples = 1, warmup = 0, chains = 2, sampler = adaptive_hmc()) - Message - running 2 chains simultaneously on up to 8 CPU cores - Output - - Condition - Error in `self$run_warmup()`: - ! object 'result' not found - # bad mcmc proposals are rejected Code - draws <- mcmc(m, chains = 1, n_samples = 2, warmup = 1, verbose = FALSE, + draws <- mcmc(m, chains = 2, n_samples = 2, warmup = 0, verbose = FALSE, sampler = adaptive_hmc(), initial_values = initials(z = 1e+120)) + Message + Only one set of initial values was provided, and was used for all chains Condition Error in `self$check_valid_parameters()`: ! The log density could not be evaluated at these initial values @@ -32,7 +22,7 @@ --- Code - mcmc(m, chains = 1, n_samples = 1, warmup = 1, sampler = adaptive_hmc(), + mcmc(m, chains = 1, n_samples = 1, warmup = 0, sampler = adaptive_hmc(), verbose = FALSE) Condition Error in `self$check_reasonable_starting_values()`: diff --git a/tests/testthat/_snaps/inference.md b/tests/testthat/_snaps/inference.md index 1fa4c7f7..fcf305e3 100644 --- a/tests/testthat/_snaps/inference.md +++ b/tests/testthat/_snaps/inference.md @@ -218,3 +218,11 @@ hmc sampler object with parameters: Lmin = 1, Lmax = 10, epsilon = 0.1, diag_sd = 1 +--- + + Code + adaptive_hmc(max_leapfrog_steps = 1L) + Output + adaptive_hmc_sampler object with parameters: + max_leapfrog_steps = 1, epsilon = 0.1, diag_sd = 1 + diff --git a/tests/testthat/helpers.R b/tests/testthat/helpers.R index 5273e3d8..b1ddd6e5 100644 --- a/tests/testthat/helpers.R +++ b/tests/testthat/helpers.R @@ -661,7 +661,8 @@ check_geweke <- function( p_x_bar_theta, niter = 2000, warmup = 1000, - thin = 1 + thin = 1, + chains = 1 ) { # sample independently target_theta <- p_theta(niter) @@ -674,7 +675,8 @@ check_geweke <- function( p_theta = p_theta, p_x_bar_theta = p_x_bar_theta, sampler = sampler, - warmup = warmup + warmup = warmup, + chains = chains ) geweke_checks <- list( @@ -716,7 +718,8 @@ p_theta_greta <- function( p_x_bar_theta, # TODO note that we might want to change this to adaptive_hmc() sampler = hmc(), - warmup = 1000 + warmup = 1000, + chains ) { # set up and initialize trace theta <- rep(NA, niter) @@ -727,7 +730,7 @@ p_theta_greta <- function( model, warmup = warmup, n_samples = 1, - chains = 1, + chains = chains, sampler = sampler, verbose = FALSE ) diff --git a/tests/testthat/test-adaptive-hmc.R b/tests/testthat/test-adaptive-hmc.R index 26c94c75..0ccc3d21 100644 --- a/tests/testthat/test-adaptive-hmc.R +++ b/tests/testthat/test-adaptive-hmc.R @@ -20,8 +20,7 @@ test_that("adaptive_hmc() does not error when given 0 warmup", { skip_if_not(check_tf_version()) x <- normal(0, 1) m <- model(x) - expect_snapshot( - error = TRUE, + expect_no_error( draws <- mcmc( m, n_samples = 1, diff --git a/tests/testthat/test_posteriors_geweke.R b/tests/testthat/test_posteriors_geweke.R index 2f2148ea..dc82a3cb 100644 --- a/tests/testthat/test_posteriors_geweke.R +++ b/tests/testthat/test_posteriors_geweke.R @@ -98,7 +98,8 @@ test_that("adaptive hmc sampler passes geweke tests", { model = model, data = x, p_theta = p_theta, - p_x_bar_theta = p_x_bar_theta + p_x_bar_theta = p_x_bar_theta, + chains = 2 ) geweke_qq(geweke_adaptive_hmc, title = "adaptive hmc sampler Geweke test") From 763a7de1ff4cb48841b0fe596e95039fead39737 Mon Sep 17 00:00:00 2001 From: njtierney Date: Wed, 19 Mar 2025 13:30:11 +1100 Subject: [PATCH 13/55] If there is no warmup, make sure the current_state gets passed through as the free state, otherwise use the last warmed up parameters to pass through to sampling. --- R/samplers.R | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/R/samplers.R b/R/samplers.R index 579dd02f..2e740fed 100644 --- a/R/samplers.R +++ b/R/samplers.R @@ -428,12 +428,21 @@ adaptive_hmc_sampler <- R6Class( # execute it result <- warmup() + # If no warmup, return the free state to use as initials in sampling, + # Otherwise, use the last chunk of all_states, which has been tuned + no_warmup <- as.integer(sampler_kernel$num_adaptation_steps) <= 0 + if (no_warmup) { + current_state <- tensorflow::as_tensor(free_state) + } else { + current_state <- get_last_state(result$all_states) + } + # return the last (burned-in) state of the model parameters and the final # (tuned) kernel parameters list( kernel = sampler_kernel, kernel_results = result$final_kernel_results, - current_state = get_last_state(result$all_states) + current_state = current_state ) }, From 2e9a2212fb6314651b1d153d26a60095190caa27 Mon Sep 17 00:00:00 2001 From: njtierney Date: Wed, 19 Mar 2025 16:02:37 +1100 Subject: [PATCH 14/55] remove test checking if adaptive_hmc errors when warmup = 0, since this is implicitly checked in other tests, e.g., extra_samples etc --- tests/testthat/test-adaptive-hmc.R | 15 --------------- 1 file changed, 15 deletions(-) diff --git a/tests/testthat/test-adaptive-hmc.R b/tests/testthat/test-adaptive-hmc.R index 0ccc3d21..59bdf389 100644 --- a/tests/testthat/test-adaptive-hmc.R +++ b/tests/testthat/test-adaptive-hmc.R @@ -16,21 +16,6 @@ test_that("adaptive_hmc() errors when given 1 chain", { ) }) -test_that("adaptive_hmc() does not error when given 0 warmup", { - skip_if_not(check_tf_version()) - x <- normal(0, 1) - m <- model(x) - expect_no_error( - draws <- mcmc( - m, - n_samples = 1, - warmup = 0, - chains = 2, - sampler = adaptive_hmc() - ) - ) -}) - test_that("bad mcmc proposals are rejected", { skip_if_not(check_tf_version()) From 76e1e95f5be117edbe62fed1b269841e1dfbc5f9 Mon Sep 17 00:00:00 2001 From: njtierney Date: Wed, 19 Mar 2025 16:03:19 +1100 Subject: [PATCH 15/55] Create functions to handle doing the update acceptance (for tuning etc), and then for rejection (for counting bad samples) --- R/sampler_class.R | 17 ++++++++++++++--- R/samplers.R | 26 +++++++++++--------------- 2 files changed, 25 insertions(+), 18 deletions(-) diff --git a/R/sampler_class.R b/R/sampler_class.R index 988228e3..4c5c867e 100644 --- a/R/sampler_class.R +++ b/R/sampler_class.R @@ -142,7 +142,6 @@ sampler <- R6Class( # how big would we like the bursts to be ideal_burst_size <- ifelse(one_by_one, 1L, pb_update) - self$run_warmup( n_samples = n_samples, pb_update = pb_update, @@ -565,14 +564,26 @@ sampler <- R6Class( self$free_state <- free_state } + self$update_acceptance(results = batch_results$trace) + self$update_rejection(results = batch_results$trace) + }, + + update_acceptance = function(results) { if (self$uses_metropolis) { # log acceptance probability - log_accept_stats <- as.array(batch_results$trace$log_accept_ratio) - is_accepted <- as.array(batch_results$trace$is_accepted) + log_accept_stats <- as.array(results$log_accept_ratio) + is_accepted <- as.array(results$is_accepted) + # related to tuning information self$accept_history <- rbind(self$accept_history, is_accepted) accept_stats_batch <- pmin(1, exp(log_accept_stats)) + # also related to tuning information self$mean_accept_stat <- mean(accept_stats_batch, na.rm = TRUE) + } + }, + update_rejection = function(results) { + if (self$uses_metropolis) { + log_accept_stats <- as.array(results$log_accept_ratio) # numerical rejections parameter sets bad <- sum(!is.finite(log_accept_stats)) self$numerical_rejections <- self$numerical_rejections + bad diff --git a/R/samplers.R b/R/samplers.R index 2e740fed..c76755d3 100644 --- a/R/samplers.R +++ b/R/samplers.R @@ -517,6 +517,7 @@ adaptive_hmc_sampler <- R6Class( sampler_kernel = sampler_kernel, free_state = init ) + self$update_rejection(result) self$warm_results <- result result @@ -606,6 +607,8 @@ adaptive_hmc_sampler <- R6Class( self$free_state <- free_state } + self$update_rejection(batch_results) + self$trace() # trace the MCMC results from this burst @@ -615,21 +618,6 @@ adaptive_hmc_sampler <- R6Class( # overwrite the current state current_state <- get_last_state(batch_results$all_states) - # accumulate and report on the badness - # new_badness <- sum(bad_steps(batch_results$kernel_results)) - # n_bad <- n_bad + new_badness - # n_evaluations <- burst * burst_size * self$n_chains - # perc_badness <- round(100 * n_bad / n_evaluations) - - # report on progress - # print(sprintf( - # "burst %i of %i (%i%s bad)", - # burst, - # n_bursts, - # perc_badness, - # "%" - # )) - if (verbose) { # update the progress bar/percentage log iterate_progress_bar( @@ -649,6 +637,14 @@ adaptive_hmc_sampler <- R6Class( } } # end sampling ### Adaptive end + }, + + update_rejection = function(results) { + if (self$uses_metropolis) { + # accumulate and report on the badness + bad <- sum(bad_steps(results$kernel_results)) + self$numerical_rejections <- self$numerical_rejections + bad + } } ) ) From 9b4617a6f3c183034b23915edfb16a73d79e6f26 Mon Sep 17 00:00:00 2001 From: njtierney Date: Thu, 20 Mar 2025 14:01:40 +1100 Subject: [PATCH 16/55] fix some of the examples in geweke test for adaptive_hmc --- tests/testthat/test_posteriors_geweke.R | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/tests/testthat/test_posteriors_geweke.R b/tests/testthat/test_posteriors_geweke.R index dc82a3cb..2db3d146 100644 --- a/tests/testthat/test_posteriors_geweke.R +++ b/tests/testthat/test_posteriors_geweke.R @@ -99,10 +99,15 @@ test_that("adaptive hmc sampler passes geweke tests", { data = x, p_theta = p_theta, p_x_bar_theta = p_x_bar_theta, - chains = 2 + chains = 2, + niter = 200, + warmup = 100, + thin = 5 ) geweke_qq(geweke_adaptive_hmc, title = "adaptive hmc sampler Geweke test") - testthat::expect_gte(geweke_adaptive_hmc$p.value, 0.005) + geweke_stat_adaptive_hmc <- geweke_ks(geweke_adaptive_hmc) + + testthat::expect_gte(geweke_stat_adaptive_hmc$p.value, 0.005) }) From 4fb8d9ad3f1aed5a382a6f876c3cd227423320fe Mon Sep 17 00:00:00 2001 From: njtierney Date: Thu, 20 Mar 2025 14:26:49 +1100 Subject: [PATCH 17/55] remove some inefficiencies inside of warm_up_sampler: - now returns as `self$warm_results` --- R/samplers.R | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/R/samplers.R b/R/samplers.R index c76755d3..367085aa 100644 --- a/R/samplers.R +++ b/R/samplers.R @@ -439,7 +439,7 @@ adaptive_hmc_sampler <- R6Class( # return the last (burned-in) state of the model parameters and the final # (tuned) kernel parameters - list( + self$warm_results <- list( kernel = sampler_kernel, kernel_results = result$final_kernel_results, current_state = current_state @@ -513,14 +513,10 @@ adaptive_hmc_sampler <- R6Class( ) { sampler_kernel <- self$define_tf_kernel() init <- self$free_state - result <- self$warm_up_sampler( + self$warm_up_sampler( sampler_kernel = sampler_kernel, free_state = init ) - self$update_rejection(result) - - self$warm_results <- result - result }, run_sampling = function( From eaeb1ddae8a8b301c708d638f3ef11759673a879 Mon Sep 17 00:00:00 2001 From: njtierney Date: Thu, 20 Mar 2025 15:55:35 +1100 Subject: [PATCH 18/55] Ensure adaptive_hmc warmup only happens when required, not every single time. - Remove extra `make_sampler_function()` - Add tests for `extra_samples` to ensure it works with `adaptive_hmc()` --- R/sampler_class.R | 15 ++++++---- R/samplers.R | 48 +++++++++++++++++++----------- R/utils.R | 41 ------------------------- tests/testthat/test-adaptive-hmc.R | 22 ++++++++++++++ 4 files changed, 61 insertions(+), 65 deletions(-) diff --git a/R/sampler_class.R b/R/sampler_class.R index 4c5c867e..5a4fda42 100644 --- a/R/sampler_class.R +++ b/R/sampler_class.R @@ -142,12 +142,15 @@ sampler <- R6Class( # how big would we like the bursts to be ideal_burst_size <- ifelse(one_by_one, 1L, pb_update) - self$run_warmup( - n_samples = n_samples, - pb_update = pb_update, - ideal_burst_size = ideal_burst_size, - verbose = verbose - ) + do_warmup <- self$warmup > 0 + if (do_warmup) { + self$run_warmup( + n_samples = n_samples, + pb_update = pb_update, + ideal_burst_size = ideal_burst_size, + verbose = verbose + ) + } self$run_sampling( n_samples = n_samples, diff --git a/R/samplers.R b/R/samplers.R index 367085aa..8ea35051 100644 --- a/R/samplers.R +++ b/R/samplers.R @@ -375,6 +375,8 @@ adaptive_hmc_sampler <- R6Class( # TODO # Eventually we will adapt the other samplers to have an interface # more in line with this - where we do the adaptation in tensorflow + sampler_kernel = NULL, + define_tf_kernel = function(sampler_param_vec) { dag <- self$model$dag tfe <- dag$tf_environment @@ -392,14 +394,10 @@ adaptive_hmc_sampler <- R6Class( ) # learns the step size - sampler_kernel <- tfp$mcmc$DualAveragingStepSizeAdaptation( + self$sampler_kernel <- tfp$mcmc$DualAveragingStepSizeAdaptation( inner_kernel = kernel_base, num_adaptation_steps = as.integer(self$warmup) ) - - return( - sampler_kernel - ) }, # given MCMC kernel `sampler_kernel` and initial model parameter state @@ -464,7 +462,16 @@ adaptive_hmc_sampler <- R6Class( # given a warmed up sampler object, return a compiled TF function # that generates a new burst of samples from samples from it - make_sampler_function = function(warm_sampler) { + sampler_function = NULL, + make_sampler_function = function(sampler) { + # if warmed up? + warmed_up <- !is.null(self$warm_results) + if (warmed_up) { + sampler <- self$warm_results + } else { + sampler <- self$sampler_kernel + } + # make the uncompiled function (with curried arguments) sample_raw <- function(current_state, n_samples) { results <- tfp$mcmc$sample_chain( @@ -473,9 +480,9 @@ adaptive_hmc_sampler <- R6Class( # where to start from current_state = current_state, # kernel - kernel = warm_sampler$kernel, + kernel = sampler$kernel, # tuned sampler settings - previous_kernel_results = warm_sampler$kernel_results, + previous_kernel_results = sampler$kernel_results, # what to trace (nothing) trace_fn = function(current_state, kernel_results) { # could compute badness here to save memory? @@ -491,15 +498,13 @@ adaptive_hmc_sampler <- R6Class( } # compile it into a concrete function and return - sample <- tensorflow::tf_function( + self$sampler_function <- tensorflow::tf_function( sample_raw, list( - as_tensorspec(warm_sampler$current_state), + as_tensorspec(sampler$current_state), tf$TensorSpec(shape = c(), dtype = tf$int32) ) ) - - sample }, # prepare a slot to put the warmed up results into @@ -511,11 +516,10 @@ adaptive_hmc_sampler <- R6Class( ideal_burst_size, verbose ) { - sampler_kernel <- self$define_tf_kernel() - init <- self$free_state + self$define_tf_kernel() self$warm_up_sampler( - sampler_kernel = sampler_kernel, - free_state = init + sampler_kernel = self$sampler_kernel, + free_state = self$free_state ) }, @@ -571,8 +575,16 @@ adaptive_hmc_sampler <- R6Class( )) completed_iterations <- cumsum(burst_lengths) + # maybe warm up a sampler + if (is.null(self$warm_results)) { + self$run_warmup() + } + # maybe compile a sampling function + if (is.null(self$sampler_function)) { + self$make_sampler_function() + } # use this to compile the warmed version - sample <- make_sampler_function(self$warm_results) + # sample <- self$make_sampler_function() current_state <- self$warm_results$current_state # trace <- array(NA, dim = c(n_samples, dim(current_state))) @@ -581,7 +593,7 @@ adaptive_hmc_sampler <- R6Class( for (burst in seq_along(burst_lengths)) { burst_size <- burst_lengths[burst] - batch_results <- sample( + batch_results <- self$sampler_function( current_state = current_state, n_samples = burst_size ) diff --git a/R/utils.R b/R/utils.R index 5284724b..2ab79af0 100644 --- a/R/utils.R +++ b/R/utils.R @@ -1120,44 +1120,3 @@ recursive_get_log_accept_ratio <- function(kernel_results) { } as.array(log_accept_ratios) } - - -# given a warmed up sampler object, return a compiled TF function that generates -# a new burst of samples from samples from it -make_sampler_function <- function(warm_sampler) { - # make the uncompiled function (with curried arguments) - sample_raw <- function(current_state, n_samples) { - results <- tfp$mcmc$sample_chain( - # how many iterations - num_results = n_samples, - # where to start from - current_state = current_state, - # kernel - kernel = warm_sampler$kernel, - # tuned sampler settings - previous_kernel_results = warm_sampler$kernel_results, - # what to trace (nothing) - trace_fn = function(current_state, kernel_results) { - # could compute badness here to save memory? - # is.finite(kernel_results$inner_results$inner_results$inner_results$log_accept_ratio) - kernel_results - } - ) - # return the parameter states and the kernel results - list( - all_states = results$all_states, - kernel_results = results$trace - ) - } - - # compile it into a concrete function and return - sample <- tensorflow::tf_function( - sample_raw, - list( - as_tensorspec(warm_sampler$current_state), - tf$TensorSpec(shape = c(), dtype = tf$int32) - ) - ) - - sample -} diff --git a/tests/testthat/test-adaptive-hmc.R b/tests/testthat/test-adaptive-hmc.R index 59bdf389..825db0bf 100644 --- a/tests/testthat/test-adaptive-hmc.R +++ b/tests/testthat/test-adaptive-hmc.R @@ -82,3 +82,25 @@ test_that("bad mcmc proposals are rejected", { verbose = FALSE )) }) + +test_that("extra_samples works with adaptive_hmc()", { + skip_if_not(check_tf_version()) + + # set up model + a <- normal(0, 1) + m <- model(a) + + draws <- mcmc( + m, + warmup = 10, + n_samples = 10, + verbose = FALSE, + sampler = adaptive_hmc() + ) + + more_draws <- extra_samples(draws, 20, verbose = FALSE) + + expect_true(inherits(more_draws, "greta_mcmc_list")) + expect_true(coda::niter(more_draws) == 30) + expect_true(coda::nchain(more_draws) == 4) +}) From b5329dc8ba6af1be2521451e09174d3ad3f2a98a Mon Sep 17 00:00:00 2001 From: njtierney Date: Fri, 21 Mar 2025 11:47:10 +1100 Subject: [PATCH 19/55] try using self$current_state where possible --- R/samplers.R | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/R/samplers.R b/R/samplers.R index 8ea35051..394248ba 100644 --- a/R/samplers.R +++ b/R/samplers.R @@ -400,6 +400,8 @@ adaptive_hmc_sampler <- R6Class( ) }, + current_state = NULL, + # given MCMC kernel `sampler_kernel` and initial model parameter state # `free_state`, adapt the kernel tuning parameters whilst simultaneously # burning-in the model parameter state. Return both finalised kernel @@ -430,9 +432,9 @@ adaptive_hmc_sampler <- R6Class( # Otherwise, use the last chunk of all_states, which has been tuned no_warmup <- as.integer(sampler_kernel$num_adaptation_steps) <= 0 if (no_warmup) { - current_state <- tensorflow::as_tensor(free_state) + self$current_state <- tensorflow::as_tensor(free_state) } else { - current_state <- get_last_state(result$all_states) + self$current_state <- get_last_state(result$all_states) } # return the last (burned-in) state of the model parameters and the final @@ -440,7 +442,7 @@ adaptive_hmc_sampler <- R6Class( self$warm_results <- list( kernel = sampler_kernel, kernel_results = result$final_kernel_results, - current_state = current_state + current_state = self$current_state ) }, @@ -501,7 +503,7 @@ adaptive_hmc_sampler <- R6Class( self$sampler_function <- tensorflow::tf_function( sample_raw, list( - as_tensorspec(sampler$current_state), + as_tensorspec(self$current_state), tf$TensorSpec(shape = c(), dtype = tf$int32) ) ) From 38d8b46b8c9d6d661fed5aeafb2aa7cfcc4159b2 Mon Sep 17 00:00:00 2001 From: njtierney Date: Fri, 21 Mar 2025 12:08:13 +1100 Subject: [PATCH 20/55] Make tweak to geweke helper function --- tests/testthat/helpers.R | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/tests/testthat/helpers.R b/tests/testthat/helpers.R index b1ddd6e5..246d6cde 100644 --- a/tests/testthat/helpers.R +++ b/tests/testthat/helpers.R @@ -756,8 +756,19 @@ p_theta_greta <- function( # samples, now using this value of x (slow, but necessary in eager mode) dag$tf_log_prob_function <- NULL dag$define_tf_log_prob_function() + sampler <- attr(draws, "model_info")$samplers[[1]] - sampler$define_tf_evaluate_sample_batch() + + # we need a condition inside the geweke test where if it's adaptive_hmc we + # recreate the log prob function, recreate the default kernel, replace that + # inside the warmed up results, and then recreate the sampler function + if (inherits(sampler, "adaptive_hmc_sampler")) { + sampler$define_tf_kernel() + sampler$warm_results$kernel <- sampler$sampler_kernel + sampler$make_sampler_function() + } else { + sampler$define_tf_evaluate_sample_batch() + } # take anoteher sample draws <- extra_samples( From 1a254540ce8486bf6f5fbeab4cae95775d2c55dd Mon Sep 17 00:00:00 2001 From: njtierney Date: Fri, 21 Mar 2025 14:48:39 +1100 Subject: [PATCH 21/55] Use the current state from the one the last time the sampler was run, not the current state from the warmed up sampler --- R/samplers.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/samplers.R b/R/samplers.R index 394248ba..9ef77aa3 100644 --- a/R/samplers.R +++ b/R/samplers.R @@ -588,7 +588,7 @@ adaptive_hmc_sampler <- R6Class( # use this to compile the warmed version # sample <- self$make_sampler_function() - current_state <- self$warm_results$current_state + current_state <- self$current_state # trace <- array(NA, dim = c(n_samples, dim(current_state))) # track numerical rejections # n_bad <- 0 From f8c03d1ef8921bc9d7b990136a762de63695925d Mon Sep 17 00:00:00 2001 From: njtierney Date: Mon, 24 Mar 2025 18:53:24 +1100 Subject: [PATCH 22/55] Some fixes for geweke tests for ahmc: - replace `current_state` with `free_state` - warm_results doesn't need to return the kernel and current state - we actually just only need the kernel_results - make_sampler_function() doesn't need `sampler` arg - make_sampler_function() also doesn't need to have the logic inside of it to dispatch different results if it is already warmed up, we can move this inside of `sample_raw` - create sampling_results object that gets the result of `batch_results$kernel_results`, for use in geweke test --- R/samplers.R | 41 +++++++++++++++++++--------------------- tests/testthat/helpers.R | 6 ++++-- 2 files changed, 23 insertions(+), 24 deletions(-) diff --git a/R/samplers.R b/R/samplers.R index 9ef77aa3..811d7b5e 100644 --- a/R/samplers.R +++ b/R/samplers.R @@ -400,7 +400,7 @@ adaptive_hmc_sampler <- R6Class( ) }, - current_state = NULL, + # current_state = NULL, # given MCMC kernel `sampler_kernel` and initial model parameter state # `free_state`, adapt the kernel tuning parameters whilst simultaneously @@ -432,17 +432,17 @@ adaptive_hmc_sampler <- R6Class( # Otherwise, use the last chunk of all_states, which has been tuned no_warmup <- as.integer(sampler_kernel$num_adaptation_steps) <= 0 if (no_warmup) { - self$current_state <- tensorflow::as_tensor(free_state) + self$free_state <- tensorflow::as_tensor(free_state) } else { - self$current_state <- get_last_state(result$all_states) + # when there is warmup, get the warmed up result + self$free_state <- get_last_state(result$all_states) } - # return the last (burned-in) state of the model parameters and the final # (tuned) kernel parameters self$warm_results <- list( - kernel = sampler_kernel, - kernel_results = result$final_kernel_results, - current_state = self$current_state + # kernel = sampler_kernel, + kernel_results = result$final_kernel_results + # current_state = self$current_state ) }, @@ -465,15 +465,7 @@ adaptive_hmc_sampler <- R6Class( # given a warmed up sampler object, return a compiled TF function # that generates a new burst of samples from samples from it sampler_function = NULL, - make_sampler_function = function(sampler) { - # if warmed up? - warmed_up <- !is.null(self$warm_results) - if (warmed_up) { - sampler <- self$warm_results - } else { - sampler <- self$sampler_kernel - } - + make_sampler_function = function() { # make the uncompiled function (with curried arguments) sample_raw <- function(current_state, n_samples) { results <- tfp$mcmc$sample_chain( @@ -482,9 +474,9 @@ adaptive_hmc_sampler <- R6Class( # where to start from current_state = current_state, # kernel - kernel = sampler$kernel, + kernel = self$sampler_kernel, # tuned sampler settings - previous_kernel_results = sampler$kernel_results, + previous_kernel_results = self$warm_results$kernel_results, # what to trace (nothing) trace_fn = function(current_state, kernel_results) { # could compute badness here to save memory? @@ -503,7 +495,7 @@ adaptive_hmc_sampler <- R6Class( self$sampler_function <- tensorflow::tf_function( sample_raw, list( - as_tensorspec(self$current_state), + as_tensorspec(tensorflow::as_tensor(self$free_state)), tf$TensorSpec(shape = c(), dtype = tf$int32) ) ) @@ -525,6 +517,8 @@ adaptive_hmc_sampler <- R6Class( ) }, + sampling_results = NULL, + run_sampling = function( n_samples, pb_update, @@ -588,7 +582,7 @@ adaptive_hmc_sampler <- R6Class( # use this to compile the warmed version # sample <- self$make_sampler_function() - current_state <- self$current_state + # current_state <- self$current_state # trace <- array(NA, dim = c(n_samples, dim(current_state))) # track numerical rejections # n_bad <- 0 @@ -596,7 +590,7 @@ adaptive_hmc_sampler <- R6Class( for (burst in seq_along(burst_lengths)) { burst_size <- burst_lengths[burst] batch_results <- self$sampler_function( - current_state = current_state, + current_state = self$free_state, n_samples = burst_size ) @@ -626,7 +620,7 @@ adaptive_hmc_sampler <- R6Class( # trace[burst_idx, , ] <- as.array(batch_results$all_states) # overwrite the current state - current_state <- get_last_state(batch_results$all_states) + # self$current_state <- get_last_state(batch_results$all_states) if (verbose) { # update the progress bar/percentage log @@ -644,6 +638,9 @@ adaptive_hmc_sampler <- R6Class( stage = "sampling" ) } + self$sampling_results <- list( + kernel_results = batch_results$kernel_results + ) } } # end sampling ### Adaptive end diff --git a/tests/testthat/helpers.R b/tests/testthat/helpers.R index 246d6cde..378f5ff5 100644 --- a/tests/testthat/helpers.R +++ b/tests/testthat/helpers.R @@ -764,13 +764,15 @@ p_theta_greta <- function( # inside the warmed up results, and then recreate the sampler function if (inherits(sampler, "adaptive_hmc_sampler")) { sampler$define_tf_kernel() - sampler$warm_results$kernel <- sampler$sampler_kernel + # sample_results <- sampler$sampling_results$kernel_results + # sampler$warm_results$kernel_results <- sample_results + # sampler$warm_results$kernel <- sampler$sampler_kernel sampler$make_sampler_function() } else { sampler$define_tf_evaluate_sample_batch() } - # take anoteher sample + # take another sample draws <- extra_samples( draws, n_samples = 1, From 15ab12f72a45d2e54095c324a2686eecc1df0ad5 Mon Sep 17 00:00:00 2001 From: njtierney Date: Mon, 24 Mar 2025 18:58:18 +1100 Subject: [PATCH 23/55] Move the sampling results into the warmed up kernel results --- tests/testthat/helpers.R | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/tests/testthat/helpers.R b/tests/testthat/helpers.R index 378f5ff5..34ff3c58 100644 --- a/tests/testthat/helpers.R +++ b/tests/testthat/helpers.R @@ -764,9 +764,8 @@ p_theta_greta <- function( # inside the warmed up results, and then recreate the sampler function if (inherits(sampler, "adaptive_hmc_sampler")) { sampler$define_tf_kernel() - # sample_results <- sampler$sampling_results$kernel_results - # sampler$warm_results$kernel_results <- sample_results - # sampler$warm_results$kernel <- sampler$sampler_kernel + sample_results <- sampler$sampling_results$kernel_results + sampler$warm_results$kernel_results <- sample_results sampler$make_sampler_function() } else { sampler$define_tf_evaluate_sample_batch() From ac5ad8a91a0722f07e8b589a88260b230dfd1e54 Mon Sep 17 00:00:00 2001 From: njtierney Date: Mon, 24 Mar 2025 19:33:11 +1100 Subject: [PATCH 24/55] breadcrumb for @goldingn --- R/samplers.R | 2 ++ 1 file changed, 2 insertions(+) diff --git a/R/samplers.R b/R/samplers.R index 811d7b5e..5a4aa2e4 100644 --- a/R/samplers.R +++ b/R/samplers.R @@ -608,6 +608,8 @@ adaptive_hmc_sampler <- R6Class( if (n_draws > 0) { free_state <- free_state_draws[n_draws, , , drop = FALSE] dim(free_state) <- dim(free_state)[-1] + ## NG - free_state and self$free_state are the same here only + ## when doing the check_geweke step, otherwise they are different self$free_state <- free_state } From 2c8a5296b9c275d57bdaf695cf1decd9e426ed58 Mon Sep 17 00:00:00 2001 From: njtierney Date: Mon, 24 Mar 2025 19:39:03 +1100 Subject: [PATCH 25/55] Another breadcrumb for @goldingn --- tests/testthat/helpers.R | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/testthat/helpers.R b/tests/testthat/helpers.R index 34ff3c58..d1f702e7 100644 --- a/tests/testthat/helpers.R +++ b/tests/testthat/helpers.R @@ -763,6 +763,7 @@ p_theta_greta <- function( # recreate the log prob function, recreate the default kernel, replace that # inside the warmed up results, and then recreate the sampler function if (inherits(sampler, "adaptive_hmc_sampler")) { + ## NG another place to place a debugger to explore things sampler$define_tf_kernel() sample_results <- sampler$sampling_results$kernel_results sampler$warm_results$kernel_results <- sample_results From ad010524ccd0e94e9fe620b94595b70d9aa5bcf5 Mon Sep 17 00:00:00 2001 From: Nick Golding Date: Tue, 25 Mar 2025 10:58:00 +0800 Subject: [PATCH 26/55] make the geweke tests less stochastic --- tests/testthat/test_posteriors_geweke.R | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/tests/testthat/test_posteriors_geweke.R b/tests/testthat/test_posteriors_geweke.R index 2db3d146..f8699d96 100644 --- a/tests/testthat/test_posteriors_geweke.R +++ b/tests/testthat/test_posteriors_geweke.R @@ -6,9 +6,10 @@ Sys.setenv("RELEASE_CANDIDATE" = "false") # for i in N n <- 10 -mu1 <- rnorm(1, 0, 3) -sd1 <- rlnorm(1) -sd2 <- rlnorm(1) + +mu1 <- 0 +sd1 <- 2 +sd2 <- 1 # prior (n draws) p_theta <- function(n) { @@ -99,9 +100,9 @@ test_that("adaptive hmc sampler passes geweke tests", { data = x, p_theta = p_theta, p_x_bar_theta = p_x_bar_theta, - chains = 2, + chains = 10, niter = 200, - warmup = 100, + warmup = 2000, thin = 5 ) From 9ca2f9b5a10c68627bf8e8156f8d5ca5c844fd5e Mon Sep 17 00:00:00 2001 From: Nick Golding Date: Tue, 25 Mar 2025 10:58:42 +0800 Subject: [PATCH 27/55] undo sampling_results edits --- R/samplers.R | 6 +----- tests/testthat/helpers.R | 2 -- 2 files changed, 1 insertion(+), 7 deletions(-) diff --git a/R/samplers.R b/R/samplers.R index 5a4aa2e4..7114280d 100644 --- a/R/samplers.R +++ b/R/samplers.R @@ -439,6 +439,7 @@ adaptive_hmc_sampler <- R6Class( } # return the last (burned-in) state of the model parameters and the final # (tuned) kernel parameters + self$warm_results <- list( # kernel = sampler_kernel, kernel_results = result$final_kernel_results @@ -517,8 +518,6 @@ adaptive_hmc_sampler <- R6Class( ) }, - sampling_results = NULL, - run_sampling = function( n_samples, pb_update, @@ -640,9 +639,6 @@ adaptive_hmc_sampler <- R6Class( stage = "sampling" ) } - self$sampling_results <- list( - kernel_results = batch_results$kernel_results - ) } } # end sampling ### Adaptive end diff --git a/tests/testthat/helpers.R b/tests/testthat/helpers.R index d1f702e7..30b03412 100644 --- a/tests/testthat/helpers.R +++ b/tests/testthat/helpers.R @@ -765,8 +765,6 @@ p_theta_greta <- function( if (inherits(sampler, "adaptive_hmc_sampler")) { ## NG another place to place a debugger to explore things sampler$define_tf_kernel() - sample_results <- sampler$sampling_results$kernel_results - sampler$warm_results$kernel_results <- sample_results sampler$make_sampler_function() } else { sampler$define_tf_evaluate_sample_batch() From 5337f8366ae4ead40638a8bfc9ab75396e272494 Mon Sep 17 00:00:00 2001 From: Nick Golding Date: Tue, 25 Mar 2025 11:15:20 +0800 Subject: [PATCH 28/55] increase adaptive HMC geweke warmup and sampling --- tests/testthat/test_posteriors_geweke.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/testthat/test_posteriors_geweke.R b/tests/testthat/test_posteriors_geweke.R index f8699d96..8df71b6d 100644 --- a/tests/testthat/test_posteriors_geweke.R +++ b/tests/testthat/test_posteriors_geweke.R @@ -101,8 +101,8 @@ test_that("adaptive hmc sampler passes geweke tests", { p_theta = p_theta, p_x_bar_theta = p_x_bar_theta, chains = 10, - niter = 200, - warmup = 2000, + niter = 1000, + warmup = 10000, thin = 5 ) From d52776be29d54d1eea0bda7a5cb61eb5a440d6cd Mon Sep 17 00:00:00 2001 From: njtierney Date: Wed, 26 Mar 2025 16:07:25 +1100 Subject: [PATCH 29/55] tidying up some "old lettuce" - commented code that doesn't need to be there. --- R/samplers.R | 30 ------------------------------ tests/testthat/helpers.R | 1 - 2 files changed, 31 deletions(-) diff --git a/R/samplers.R b/R/samplers.R index 7114280d..f4afecee 100644 --- a/R/samplers.R +++ b/R/samplers.R @@ -441,9 +441,7 @@ adaptive_hmc_sampler <- R6Class( # (tuned) kernel parameters self$warm_results <- list( - # kernel = sampler_kernel, kernel_results = result$final_kernel_results - # current_state = self$current_state ) }, @@ -457,8 +455,6 @@ adaptive_hmc_sampler <- R6Class( # return named list for replacing tensors list( adaptive_hmc_max_leapfrog_steps = max_leapfrog_steps, - # adaptive_hmc_epsilon = epsilon, - # adaptive_hmc_diag_sd = diag_sd, method = method ) }, @@ -526,16 +522,6 @@ adaptive_hmc_sampler <- R6Class( thin, verbose ) { - # There's a couple of possibilities - ## warmup already happened previously, so get already-warmup sampler - ## warmup never happened, so get an unwarmed up sampler - - # warmup_status <- self$warmup_status() - # self$warm_results - # self$warmup - # warmup_has_happened <- self$has_warmup_happened() - # warmup - perform_sampling <- n_samples > 0 if (perform_sampling) { # on exiting during the main sampling period (even if killed by the @@ -578,13 +564,6 @@ adaptive_hmc_sampler <- R6Class( if (is.null(self$sampler_function)) { self$make_sampler_function() } - # use this to compile the warmed version - # sample <- self$make_sampler_function() - - # current_state <- self$current_state - # trace <- array(NA, dim = c(n_samples, dim(current_state))) - # track numerical rejections - # n_bad <- 0 for (burst in seq_along(burst_lengths)) { burst_size <- burst_lengths[burst] @@ -607,8 +586,6 @@ adaptive_hmc_sampler <- R6Class( if (n_draws > 0) { free_state <- free_state_draws[n_draws, , , drop = FALSE] dim(free_state) <- dim(free_state)[-1] - ## NG - free_state and self$free_state are the same here only - ## when doing the check_geweke step, otherwise they are different self$free_state <- free_state } @@ -616,13 +593,6 @@ adaptive_hmc_sampler <- R6Class( self$trace() - # trace the MCMC results from this burst - # burst_idx <- (burst - 1) * burst_size + seq_len(burst_size) - # trace[burst_idx, , ] <- as.array(batch_results$all_states) - - # overwrite the current state - # self$current_state <- get_last_state(batch_results$all_states) - if (verbose) { # update the progress bar/percentage log iterate_progress_bar( diff --git a/tests/testthat/helpers.R b/tests/testthat/helpers.R index 30b03412..7d17c318 100644 --- a/tests/testthat/helpers.R +++ b/tests/testthat/helpers.R @@ -763,7 +763,6 @@ p_theta_greta <- function( # recreate the log prob function, recreate the default kernel, replace that # inside the warmed up results, and then recreate the sampler function if (inherits(sampler, "adaptive_hmc_sampler")) { - ## NG another place to place a debugger to explore things sampler$define_tf_kernel() sampler$make_sampler_function() } else { From 193763a65003722028183dcc16a3c84dd32b4ee8 Mon Sep 17 00:00:00 2001 From: njtierney Date: Thu, 27 Mar 2025 13:39:24 +1100 Subject: [PATCH 30/55] first pass at HMC sampler warmup message --- R/samplers.R | 2 ++ 1 file changed, 2 insertions(+) diff --git a/R/samplers.R b/R/samplers.R index f4afecee..c4fc3cdc 100644 --- a/R/samplers.R +++ b/R/samplers.R @@ -508,10 +508,12 @@ adaptive_hmc_sampler <- R6Class( verbose ) { self$define_tf_kernel() + cli::cli_inform("warming up HMC sampler") self$warm_up_sampler( sampler_kernel = self$sampler_kernel, free_state = self$free_state ) + cli::cli_inform("HMC sampler warmed up!") }, run_sampling = function( From 417849b6e8ae15c45b676b76e1b7fbbf97398a97 Mon Sep 17 00:00:00 2001 From: njtierney Date: Thu, 27 Mar 2025 15:08:00 +1100 Subject: [PATCH 31/55] add warmup info for adaptive hmc --- R/samplers.R | 37 +++++++++++++++++++++++++++++++++++-- 1 file changed, 35 insertions(+), 2 deletions(-) diff --git a/R/samplers.R b/R/samplers.R index c4fc3cdc..208385b2 100644 --- a/R/samplers.R +++ b/R/samplers.R @@ -507,13 +507,46 @@ adaptive_hmc_sampler <- R6Class( ideal_burst_size, verbose ) { + if (verbose) { + pb_warmup <- create_progress_bar( + phase = "warmup", + iter = c(self$warmup, n_samples), + pb_update = self$warmup, + width = self$pb_width + ) + iterate_progress_bar( + pb = pb_warmup, + it = 0, + rejects = 0, + chains = self$n_chains, + file = self$pb_file + ) + } else { + pb_sampling <- NULL + } self$define_tf_kernel() - cli::cli_inform("warming up HMC sampler") self$warm_up_sampler( sampler_kernel = self$sampler_kernel, free_state = self$free_state ) - cli::cli_inform("HMC sampler warmed up!") + + if (verbose) { + # update the progress bar/percentage log + iterate_progress_bar( + pb = pb_warmup, + it = self$warmup, + # TODO check if we can add rejection info + rejects = 0, + chains = self$n_chains, + file = self$pb_file + ) + + self$write_percentage_log( + total = self$warmup, + completed = self$warmup, + stage = "warmup" + ) + } # close verbose }, run_sampling = function( From adc048638a907c2cc70bc9975dbed10db91f3bfe Mon Sep 17 00:00:00 2001 From: njtierney Date: Thu, 24 Apr 2025 13:59:18 +1000 Subject: [PATCH 32/55] update WORDLIST --- inst/WORDLIST | 16 ++-------------- 1 file changed, 2 insertions(+), 14 deletions(-) diff --git a/inst/WORDLIST b/inst/WORDLIST index 731a7864..89d589d8 100644 --- a/inst/WORDLIST +++ b/inst/WORDLIST @@ -30,23 +30,22 @@ Nesterov NumPy ORCID OpenBUGS -Optimizers PSAT ProximalAdagradOptimizer ProximalGradientDescentOptimizer README RStudio RStudio's +SNAPER Scalable TF TFP Tensorflow -Vectorised WinBUGS XLA -acyclic al analysed +arg args automagically bayesplot @@ -106,7 +105,6 @@ optimisation optimisations optimiser optimisers -optimizers ouputs parallelisation parallelise @@ -114,19 +112,11 @@ parallelising pkgdown poisson polygamma -posteriori -priori pythonic realisation realisations -recoded regularisation reinstalls -repo -reproducibility -rescale -rescaled -rescaling reticulate scalable schoolers @@ -139,12 +129,10 @@ tf tfp tril tuh -unadjusted uncentered uncoached vectorising visualise -warmup winbuilder wishart zhang From dd46ff40aa70f1f40e384df289c346e701c02809 Mon Sep 17 00:00:00 2001 From: njtierney Date: Thu, 24 Apr 2025 14:02:51 +1000 Subject: [PATCH 33/55] Various tweaks: * Fix typos (proability and verbsoity) * add "as_tensor_spec" to globalVariables * update local snapshot tests (only works for me, but worthwhile checking --- R/greta-sitrep.R | 4 ++-- R/package.R | 1 + man/greta_sitrep.Rd | 4 ++-- man/samplers.Rd | 8 ++++---- tests/testthat/_snaps/adaptive-hmc.md | 3 +++ tests/testthat/_snaps/greta-sitrep.md | 8 ++++---- 6 files changed, 16 insertions(+), 12 deletions(-) diff --git a/R/greta-sitrep.R b/R/greta-sitrep.R index 46f63abd..4240cd5b 100644 --- a/R/greta-sitrep.R +++ b/R/greta-sitrep.R @@ -6,12 +6,12 @@ #' @param verbosity character. How verbose the output of the situation report. #' Possible options: "minimal" (default), "detailed", and "quiet". "Minimal" #' provides just information in python version, tensorflow version, -#' tensorflow proability, and whether greta conda environment is available. +#' tensorflow probability, and whether greta conda environment is available. #' "Quiet" presents no information, but prepares greta to be used. "Detailed" #' gives information on the version and path for R, greta, python, #' tensorflow, tensorflow probability, the greta conda environment, and a #' statement on greta usability. -#' @return Message on greta situation report. See "verbsoity" parameter details +#' @return Message on greta situation report. See "verbosity" parameter details #' above for more information. #' @export #' diff --git a/R/package.R b/R/package.R index 3b6df155..6ceec328 100644 --- a/R/package.R +++ b/R/package.R @@ -47,6 +47,7 @@ reticulate::py_version #' @importFrom utils globalVariables utils::globalVariables( c( + "as_tensor_spec", "N", "greta_deps_tf_tfp", "greta_logfile", diff --git a/man/greta_sitrep.Rd b/man/greta_sitrep.Rd index 8fd90b95..79716fad 100644 --- a/man/greta_sitrep.Rd +++ b/man/greta_sitrep.Rd @@ -10,14 +10,14 @@ greta_sitrep(verbosity = c("minimal", "detailed", "quiet")) \item{verbosity}{character. How verbose the output of the situation report. Possible options: "minimal" (default), "detailed", and "quiet". "Minimal" provides just information in python version, tensorflow version, -tensorflow proability, and whether greta conda environment is available. +tensorflow probability, and whether greta conda environment is available. "Quiet" presents no information, but prepares greta to be used. "Detailed" gives information on the version and path for R, greta, python, tensorflow, tensorflow probability, the greta conda environment, and a statement on greta usability.} } \value{ -Message on greta situation report. See "verbsoity" parameter details +Message on greta situation report. See "verbosity" parameter details above for more information. } \description{ diff --git a/man/samplers.Rd b/man/samplers.Rd index e5c825f3..9a7e9f64 100644 --- a/man/samplers.Rd +++ b/man/samplers.Rd @@ -68,8 +68,8 @@ is implemented for uniform and normal proposals. \code{slice()} implements a multivariate slice sampling algorithm. The parameter \code{max_doublings} is not tuned during warmup. -For \code{adaptive_hmc()}. The Lmin and Lmax parameters are learnt and so -not provided in this. The number of chains cannot be less than 2, due to -how adaptive HMC works. \code{diag_sd} is used to rescale the parameter space to -make it more uniform, and make sampling more efficient. +For \code{adaptive_hmc()}. The Lmin and Lmax parameters are learned and +so not provided in this. The number of chains cannot be less than 2, due +to how adaptive HMC works. \code{diag_sd} is used to rescale the parameter +space to make it more uniform, and make sampling more efficient. } diff --git a/tests/testthat/_snaps/adaptive-hmc.md b/tests/testthat/_snaps/adaptive-hmc.md index f8a67ab7..2b3d1601 100644 --- a/tests/testthat/_snaps/adaptive-hmc.md +++ b/tests/testthat/_snaps/adaptive-hmc.md @@ -2,6 +2,9 @@ Code draws <- mcmc(m, n_samples = 1, warmup = 1, chains = 1, sampler = adaptive_hmc()) + Message + + warmup 0/1 | eta: ?s Condition Error in `py_call_impl()`: ! ValueError: SNAPERHMC requires at least 2 chains. Got: 1 diff --git a/tests/testthat/_snaps/greta-sitrep.md b/tests/testthat/_snaps/greta-sitrep.md index a0ac6579..10031eb7 100644 --- a/tests/testthat/_snaps/greta-sitrep.md +++ b/tests/testthat/_snaps/greta-sitrep.md @@ -246,18 +246,18 @@ Message -- R --------------------------------------------------------------------------- - * version: 4.4.2 + * version: 4.5.0 * path: '/Library/Frameworks/R.framework/Resources' -- greta ----------------------------------------------------------------------- * version: 0.5.0.9000 - * path: '/Users/nick/github/greta-dev/greta' + * path: '/Users/nick_1/github/greta-dev/greta' -- python ---------------------------------------------------------------------- i checking if python available v python (v3.10) available - * path: '/Users/nick/Library/r-miniconda-arm64' + * path: '/Users/nick_1/Library/r-miniconda-arm64' -- greta conda environment ----------------------------------------------------- i checking if greta conda environment available @@ -269,7 +269,7 @@ i checking if TensorFlow available v TensorFlow (v2.15.0) available - * R path: '/Users/nick/Library/R/arm64/4.4/library/tensorflow' + * R path: '/Users/nick_1/Library/R/arm64/4.5/library/tensorflow' * Exists in conda env: TRUE -- TensorFlow Probability ------------------------------------------------------ From 46cd71cddb1c93257c4899896ea09a2f01a24948 Mon Sep 17 00:00:00 2001 From: njtierney Date: Thu, 24 Apr 2025 14:04:05 +1000 Subject: [PATCH 34/55] When adaptive HMC warmup is 0, do not print warmup progress bars. --- R/samplers.R | 20 +++++++++++++------- 1 file changed, 13 insertions(+), 7 deletions(-) diff --git a/R/samplers.R b/R/samplers.R index 208385b2..51f5da8b 100644 --- a/R/samplers.R +++ b/R/samplers.R @@ -103,10 +103,10 @@ slice <- function(max_doublings = 5) { #' steps used. The algorithm will determine the optimal number less than this. #' @param method character length one. Currently can only be "SNAPER" but in #' the future this may expand to other adaptive samplers. -#' @details For `adaptive_hmc()`. The Lmin and Lmax parameters are learnt and so -#' not provided in this. The number of chains cannot be less than 2, due to -#' how adaptive HMC works. `diag_sd` is used to rescale the parameter space to -#' make it more uniform, and make sampling more efficient. +#' @details For `adaptive_hmc()`. The Lmin and Lmax parameters are learned and +#' so not provided in this. The number of chains cannot be less than 2, due +#' to how adaptive HMC works. `diag_sd` is used to rescale the parameter +#' space to make it more uniform, and make sampling more efficient. #' @export adaptive_hmc <- function( max_leapfrog_steps = 1000, @@ -507,7 +507,9 @@ adaptive_hmc_sampler <- R6Class( ideal_burst_size, verbose ) { - if (verbose) { + # skip warmup progress bar if warmup is 0 + nonzero_warmup <- self$warmup > 0 + if (verbose && nonzero_warmup) { pb_warmup <- create_progress_bar( phase = "warmup", iter = c(self$warmup, n_samples), @@ -530,7 +532,8 @@ adaptive_hmc_sampler <- R6Class( free_state = self$free_state ) - if (verbose) { + # skip warmup progress bar if warmup is 0 + if (verbose && nonzero_warmup) { # update the progress bar/percentage log iterate_progress_bar( pb = pb_warmup, @@ -593,7 +596,10 @@ adaptive_hmc_sampler <- R6Class( # maybe warm up a sampler if (is.null(self$warm_results)) { - self$run_warmup() + self$run_warmup( + verbose = verbose, + n_samples = n_samples + ) } # maybe compile a sampling function if (is.null(self$sampler_function)) { From efdea40b8dff0afecf23f672fed43ce29b692362 Mon Sep 17 00:00:00 2001 From: njtierney Date: Thu, 24 Apr 2025 17:27:38 +1000 Subject: [PATCH 35/55] see if this helps resolve an obscure windows GH error --- .github/workflows/R-CMD-check.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml index d0cbcc13..ab2567c4 100644 --- a/.github/workflows/R-CMD-check.yaml +++ b/.github/workflows/R-CMD-check.yaml @@ -32,7 +32,7 @@ jobs: continue-on-error: ${{ matrix.tf == 'nightly' || contains(matrix.tf, 'rc') || matrix.r == 'devel' }} env: R_REMOTES_NO_ERRORS_FROM_WARNINGS: 'true' - R_COMPILE_AND_INSTALL_PACKAGES: 'never' + # R_COMPILE_AND_INSTALL_PACKAGES: 'never' GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} steps: From 5227b5267261e206874f555fc1612365b3a9d6c5 Mon Sep 17 00:00:00 2001 From: njtierney Date: Tue, 29 Apr 2025 10:20:06 +1000 Subject: [PATCH 36/55] tweak GH action --- .github/workflows/R-CMD-check.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml index ab2567c4..d0cbcc13 100644 --- a/.github/workflows/R-CMD-check.yaml +++ b/.github/workflows/R-CMD-check.yaml @@ -32,7 +32,7 @@ jobs: continue-on-error: ${{ matrix.tf == 'nightly' || contains(matrix.tf, 'rc') || matrix.r == 'devel' }} env: R_REMOTES_NO_ERRORS_FROM_WARNINGS: 'true' - # R_COMPILE_AND_INSTALL_PACKAGES: 'never' + R_COMPILE_AND_INSTALL_PACKAGES: 'never' GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} steps: From 6ce96e3aaa829408fe7ce0c8c288d1406249b2bd Mon Sep 17 00:00:00 2001 From: njtierney Date: Tue, 29 Apr 2025 12:05:37 +1000 Subject: [PATCH 37/55] try out the geweke test script --- .github/workflows/geweke.yaml | 167 ++++++++++++++++++++++++++++++++++ 1 file changed, 167 insertions(+) create mode 100644 .github/workflows/geweke.yaml diff --git a/.github/workflows/geweke.yaml b/.github/workflows/geweke.yaml new file mode 100644 index 00000000..21762f7f --- /dev/null +++ b/.github/workflows/geweke.yaml @@ -0,0 +1,167 @@ +name: Geweke Tests + +on: + push: + branches: + - main + - master + pull_request: + branches: + - main + - master + workflow_dispatch: + inputs: + script_path: + description: 'Path to Geweke test script' + required: true + default: 'inst/greta-geweke-ahmc.R' + +jobs: + run-geweke-tests: + name: ${{ matrix.os }}, tf-${{ matrix.tf }}, R-${{ matrix.r}} + timeout-minutes: 45 + strategy: + fail-fast: false + matrix: + include: + # Uncomment additional environments as needed + # - {os: 'ubuntu-latest' , tf: 'default', r: 'release'} + # - {os: 'windows-latest', tf: 'default', r: 'release'} + - {os: 'macOS-latest' , tf: 'default', r: 'release'} + + runs-on: ${{ matrix.os }} + continue-on-error: ${{ matrix.tf == 'nightly' || contains(matrix.tf, 'rc') || matrix.r == 'devel' }} + env: + R_REMOTES_NO_ERRORS_FROM_WARNINGS: 'true' + R_COMPILE_AND_INSTALL_PACKAGES: 'never' + GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} + + steps: + - uses: actions/checkout@v4 + + - uses: r-lib/actions/setup-pandoc@v2 + + - uses: r-lib/actions/setup-r@v2 + id: setup-r + with: + r-version: ${{ matrix.r }} + Ncpus: '4L' + use-public-rspm: true + + - name: Get Date + id: get-date + shell: bash + run: | + echo "year-week=$(date -u '+%Y-%U')" >> $GITHUB_OUTPUT + echo "date=$(date -u '+%F')" >> $GITHUB_OUTPUT + + - name: Restore R package cache + uses: actions/cache@v4 + id: r-package-cache + with: + path: ${{ env.R_LIBS_USER }} + key: ${{ matrix.os }}-${{ steps.setup-r.outputs.installed-r-version }}-${{ steps.get-date.outputs.year-week }}-1 + + - uses: r-lib/actions/setup-r-dependencies@v2 + with: + extra-packages: any::rcmdcheck any::devtools any::reticulate local::. + cache-version: 4 + upgrade: 'TRUE' + + - name: Install greta deps + shell: Rscript {0} + run: | + library(greta) + greta::install_greta_deps(timeout = 50) + + - name: Situation Report on greta install + shell: Rscript {0} + run: greta::greta_sitrep() + + - name: Session info + shell: Rscript {0} + run: | + options(width = 100) + pkgs <- installed.packages()[, "Package"] + sessioninfo::session_info(pkgs, include_base = TRUE) + + - name: Create plots directory + shell: bash + run: mkdir -p plots + + - name: Set up wrapper script + shell: bash + run: | + # Create the wrapper script in the tools directory + mkdir -p tools + cat > tools/run_geweke_test.R << 'EOF' +#!/usr/bin/env Rscript + +# Script that runs a Geweke test and captures plots and statistics +# Usage: Rscript run_geweke_test.R path/to/geweke/script.R + +# Get script path from command line args +args <- commandArgs(trailingOnly = TRUE) +script_path <- if (length(args) > 0) args[1] else "inst/greta-geweke-ahmc.R" + +cat("Running Geweke test script:", script_path, "\n") + +# Create directory for plots if it doesn't exist +if (!dir.exists("plots")) dir.create("plots", recursive = TRUE) + +# Load the package +devtools::load_all(".") + +# Set up graphics device to capture plots +script_basename <- basename(tools::file_path_sans_ext(script_path)) +plot_filename <- paste0("plots/geweke_", script_basename, ".png") +png(plot_filename, width = 1200, height = 800, res = 150) + +# Source the user's script +tryCatch({ + source(script_path) + cat("Successfully sourced script:", script_path, "\n") +}, error = function(e) { + cat("Error occurred while running the script:\n") + cat(as.character(e), "\n") + quit(status = 1) +}) + +# Close the device to save the plot +dev.off() +cat("Plot saved to", plot_filename, "\n") + +# Save test statistics to a text file +stats_filename <- paste0("plots/geweke_stats_", script_basename, ".txt") +if (exists("geweke_stat_adaptive_hmc")) { + capture.output(geweke_stat_adaptive_hmc, file = stats_filename) + cat("Geweke test statistics saved to", stats_filename, "\n") +} else { + cat("Warning: No geweke_stat_adaptive_hmc object found\n") +} + +# List all created files +cat("Files in plots directory:\n") +print(list.files("plots")) +EOF + chmod +x tools/run_geweke_test.R + + - name: Run Geweke test script + shell: bash + run: | + script_path="${{ github.event.inputs.script_path || 'inst/greta-geweke-ahmc.R' }}" + Rscript tools/run_geweke_test.R "$script_path" + + - name: Upload plots and statistics as artifacts + uses: actions/upload-artifact@v3 + with: + name: geweke-results-${{ matrix.os }} + path: plots/ + if-no-files-found: warn + + - name: Upload test script + uses: actions/upload-artifact@v3 + with: + name: geweke-script-${{ matrix.os }} + path: ${{ github.event.inputs.script_path || 'inst/greta-geweke-ahmc.R' }} + if-no-files-found: warn From 1703180449b7285f40fd90300df6ca26556dac2f Mon Sep 17 00:00:00 2001 From: njtierney Date: Tue, 29 Apr 2025 15:42:15 +1000 Subject: [PATCH 38/55] Try out getting the Geweke test to run --- .github/workflows/geweke.yaml | 141 +++++----------------------------- 1 file changed, 19 insertions(+), 122 deletions(-) diff --git a/.github/workflows/geweke.yaml b/.github/workflows/geweke.yaml index 21762f7f..c735806b 100644 --- a/.github/workflows/geweke.yaml +++ b/.github/workflows/geweke.yaml @@ -1,4 +1,4 @@ -name: Geweke Tests +name: Geweke-tests on: push: @@ -10,51 +10,25 @@ on: - main - master workflow_dispatch: - inputs: - script_path: - description: 'Path to Geweke test script' - required: true - default: 'inst/greta-geweke-ahmc.R' jobs: - run-geweke-tests: - name: ${{ matrix.os }}, tf-${{ matrix.tf }}, R-${{ matrix.r}} - timeout-minutes: 45 - strategy: - fail-fast: false - matrix: - include: - # Uncomment additional environments as needed - # - {os: 'ubuntu-latest' , tf: 'default', r: 'release'} - # - {os: 'windows-latest', tf: 'default', r: 'release'} - - {os: 'macOS-latest' , tf: 'default', r: 'release'} - - runs-on: ${{ matrix.os }} - continue-on-error: ${{ matrix.tf == 'nightly' || contains(matrix.tf, 'rc') || matrix.r == 'devel' }} - env: - R_REMOTES_NO_ERRORS_FROM_WARNINGS: 'true' - R_COMPILE_AND_INSTALL_PACKAGES: 'never' - GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} - + run: + runs-on: ubuntu-latest + permissions: + contents: write steps: - uses: actions/checkout@v4 - - - uses: r-lib/actions/setup-pandoc@v2 - - uses: r-lib/actions/setup-r@v2 id: setup-r with: r-version: ${{ matrix.r }} Ncpus: '4L' use-public-rspm: true - - - name: Get Date - id: get-date - shell: bash - run: | - echo "year-week=$(date -u '+%Y-%U')" >> $GITHUB_OUTPUT - echo "date=$(date -u '+%F')" >> $GITHUB_OUTPUT - + - uses: r-lib/actions/setup-r-dependencies@v2\ + with: + extra-packages: any::rcmdcheck any::devtools any::reticulate local::. + cache-version: 4 + upgrade: 'TRUE' - name: Restore R package cache uses: actions/cache@v4 id: r-package-cache @@ -62,12 +36,6 @@ jobs: path: ${{ env.R_LIBS_USER }} key: ${{ matrix.os }}-${{ steps.setup-r.outputs.installed-r-version }}-${{ steps.get-date.outputs.year-week }}-1 - - uses: r-lib/actions/setup-r-dependencies@v2 - with: - extra-packages: any::rcmdcheck any::devtools any::reticulate local::. - cache-version: 4 - upgrade: 'TRUE' - - name: Install greta deps shell: Rscript {0} run: | @@ -78,6 +46,9 @@ jobs: shell: Rscript {0} run: greta::greta_sitrep() + - name: Run Geweke tests + run: Rscript inst/greta-geweke-ahmc.R + - name: Session info shell: Rscript {0} run: | @@ -85,83 +56,9 @@ jobs: pkgs <- installed.packages()[, "Package"] sessioninfo::session_info(pkgs, include_base = TRUE) - - name: Create plots directory - shell: bash - run: mkdir -p plots - - - name: Set up wrapper script - shell: bash - run: | - # Create the wrapper script in the tools directory - mkdir -p tools - cat > tools/run_geweke_test.R << 'EOF' -#!/usr/bin/env Rscript - -# Script that runs a Geweke test and captures plots and statistics -# Usage: Rscript run_geweke_test.R path/to/geweke/script.R - -# Get script path from command line args -args <- commandArgs(trailingOnly = TRUE) -script_path <- if (length(args) > 0) args[1] else "inst/greta-geweke-ahmc.R" - -cat("Running Geweke test script:", script_path, "\n") - -# Create directory for plots if it doesn't exist -if (!dir.exists("plots")) dir.create("plots", recursive = TRUE) - -# Load the package -devtools::load_all(".") - -# Set up graphics device to capture plots -script_basename <- basename(tools::file_path_sans_ext(script_path)) -plot_filename <- paste0("plots/geweke_", script_basename, ".png") -png(plot_filename, width = 1200, height = 800, res = 150) - -# Source the user's script -tryCatch({ - source(script_path) - cat("Successfully sourced script:", script_path, "\n") -}, error = function(e) { - cat("Error occurred while running the script:\n") - cat(as.character(e), "\n") - quit(status = 1) -}) - -# Close the device to save the plot -dev.off() -cat("Plot saved to", plot_filename, "\n") - -# Save test statistics to a text file -stats_filename <- paste0("plots/geweke_stats_", script_basename, ".txt") -if (exists("geweke_stat_adaptive_hmc")) { - capture.output(geweke_stat_adaptive_hmc, file = stats_filename) - cat("Geweke test statistics saved to", stats_filename, "\n") -} else { - cat("Warning: No geweke_stat_adaptive_hmc object found\n") -} - -# List all created files -cat("Files in plots directory:\n") -print(list.files("plots")) -EOF - chmod +x tools/run_geweke_test.R - - - name: Run Geweke test script - shell: bash - run: | - script_path="${{ github.event.inputs.script_path || 'inst/greta-geweke-ahmc.R' }}" - Rscript tools/run_geweke_test.R "$script_path" - - - name: Upload plots and statistics as artifacts - uses: actions/upload-artifact@v3 - with: - name: geweke-results-${{ matrix.os }} - path: plots/ - if-no-files-found: warn - - - name: Upload test script - uses: actions/upload-artifact@v3 - with: - name: geweke-script-${{ matrix.os }} - path: ${{ github.event.inputs.script_path || 'inst/greta-geweke-ahmc.R' }} - if-no-files-found: warn + # - name: Upload plots and statistics as artifacts + # uses: actions/upload-artifact@v3 + # with: + # name: geweke-results-${{ matrix.os }} + # path: plots/ + # if-no-files-found: warn From 6d6b7d789158a3517ec537bf62ec3a420b78a3bb Mon Sep 17 00:00:00 2001 From: njtierney Date: Tue, 29 Apr 2025 15:45:08 +1000 Subject: [PATCH 39/55] fix typo --- .github/workflows/geweke.yaml | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/.github/workflows/geweke.yaml b/.github/workflows/geweke.yaml index c735806b..8e1cb437 100644 --- a/.github/workflows/geweke.yaml +++ b/.github/workflows/geweke.yaml @@ -16,15 +16,14 @@ jobs: runs-on: ubuntu-latest permissions: contents: write + env: + GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} steps: - uses: actions/checkout@v4 - uses: r-lib/actions/setup-r@v2 - id: setup-r with: - r-version: ${{ matrix.r }} - Ncpus: '4L' use-public-rspm: true - - uses: r-lib/actions/setup-r-dependencies@v2\ + - uses: r-lib/actions/setup-r-dependencies@v2 with: extra-packages: any::rcmdcheck any::devtools any::reticulate local::. cache-version: 4 From 681bd73a0412967038a4d3c052be6e2eeae72c38 Mon Sep 17 00:00:00 2001 From: njtierney Date: Tue, 29 Apr 2025 15:55:09 +1000 Subject: [PATCH 40/55] oops, push code for GHA to run --- .github/workflows/geweke.yaml | 12 +----- inst/greta-geweke-ahmc.R | 79 +++++++++++++++++++++++++++++++++++ 2 files changed, 81 insertions(+), 10 deletions(-) create mode 100644 inst/greta-geweke-ahmc.R diff --git a/.github/workflows/geweke.yaml b/.github/workflows/geweke.yaml index 8e1cb437..9a15da44 100644 --- a/.github/workflows/geweke.yaml +++ b/.github/workflows/geweke.yaml @@ -19,7 +19,8 @@ jobs: env: GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} steps: - - uses: actions/checkout@v4 + - name: Checkout + uses: actions/checkout@v4 - uses: r-lib/actions/setup-r@v2 with: use-public-rspm: true @@ -28,23 +29,14 @@ jobs: extra-packages: any::rcmdcheck any::devtools any::reticulate local::. cache-version: 4 upgrade: 'TRUE' - - name: Restore R package cache - uses: actions/cache@v4 - id: r-package-cache - with: - path: ${{ env.R_LIBS_USER }} - key: ${{ matrix.os }}-${{ steps.setup-r.outputs.installed-r-version }}-${{ steps.get-date.outputs.year-week }}-1 - - name: Install greta deps shell: Rscript {0} run: | library(greta) greta::install_greta_deps(timeout = 50) - - name: Situation Report on greta install shell: Rscript {0} run: greta::greta_sitrep() - - name: Run Geweke tests run: Rscript inst/greta-geweke-ahmc.R diff --git a/inst/greta-geweke-ahmc.R b/inst/greta-geweke-ahmc.R new file mode 100644 index 00000000..c4683c02 --- /dev/null +++ b/inst/greta-geweke-ahmc.R @@ -0,0 +1,79 @@ +devtools::load_all(".") + +# Set up graphics device to capture plots +# Create directory for plots if it doesn't exist +# if (!dir.exists("plots")) dir.create("plots", recursive = TRUE) +# +# script_basename <- basename(tools::file_path_sans_ext(script_path)) +# plot_filename <- paste0("plots/geweke_", script_basename, ".png") +# png(plot_filename, width = 1200, height = 800, res = 150) +# +# # Source the user's script +# tryCatch( +# { +# source(script_path) +# cat("Successfully sourced script:", script_path, "\n") +# }, +# error = function(e) { +# cat("Error occurred while running the script:\n") +# cat(as.character(e), "\n") +# quit(status = 1) +# } +# ) +# +# # Close the device to save the plot +# dev.off() +# cat("Plot saved to", plot_filename, "\n") +# +# # Save test statistics to a text file +# stats_filename <- paste0("plots/geweke_stats_", script_basename, ".txt") +# if (exists("geweke_stat_adaptive_hmc")) { +# capture.output(geweke_stat_adaptive_hmc, file = stats_filename) +# cat("Geweke test statistics saved to", stats_filename, "\n") +# } else { +# cat("Warning: No geweke_stat_adaptive_hmc object found\n") +# } +# +# # List all created files +# cat("Files in plots directory:\n") +# print(list.files("plots")) + +n <- 10 + +mu1 <- 0 +sd1 <- 2 +sd2 <- 1 + +# prior (n draws) +p_theta <- function(n) { + rnorm(n, mu1, sd1) +} + +# likelihood +p_x_bar_theta <- function(theta) { + rnorm(n, theta, sd2) +} + +# define the greta model (single precision for slice sampler) +x <- as_data(rep(0, n)) +greta_theta <- normal(mu1, sd1) +distribution(x) <- normal(greta_theta, sd2) +model <- model(greta_theta, precision = "single") + +geweke_adaptive_hmc <- check_geweke( + sampler = adaptive_hmc(), + model = model, + data = x, + p_theta = p_theta, + p_x_bar_theta = p_x_bar_theta, + chains = 10, + niter = 100, + warmup = 100, + thin = 5 +) + +# geweke_qq(geweke_adaptive_hmc, title = "adaptive hmc sampler Geweke test") +# +# geweke_stat_adaptive_hmc <- geweke_ks(geweke_adaptive_hmc) + +geweke_stat_adaptive_hmc From f6fe8bae5803b2471fa3e04b470a0de2e938fa24 Mon Sep 17 00:00:00 2001 From: njtierney Date: Tue, 29 Apr 2025 16:30:37 +1000 Subject: [PATCH 41/55] upload plot onto action --- .github/workflows/geweke.yaml | 7 +++++++ inst/greta-geweke-ahmc.R | 33 +++++++++++++++++++-------------- 2 files changed, 26 insertions(+), 14 deletions(-) diff --git a/.github/workflows/geweke.yaml b/.github/workflows/geweke.yaml index 9a15da44..9f29da11 100644 --- a/.github/workflows/geweke.yaml +++ b/.github/workflows/geweke.yaml @@ -40,6 +40,13 @@ jobs: - name: Run Geweke tests run: Rscript inst/greta-geweke-ahmc.R + - name: Upload plot as artifact + uses: actions/upload-artifact@v3 + with: + name: qq-plot + path: geweke-plot.png + retention-days: 90 # How long to keep the artifact + - name: Session info shell: Rscript {0} run: | diff --git a/inst/greta-geweke-ahmc.R b/inst/greta-geweke-ahmc.R index c4683c02..f113d236 100644 --- a/inst/greta-geweke-ahmc.R +++ b/inst/greta-geweke-ahmc.R @@ -60,20 +60,25 @@ greta_theta <- normal(mu1, sd1) distribution(x) <- normal(greta_theta, sd2) model <- model(greta_theta, precision = "single") -geweke_adaptive_hmc <- check_geweke( - sampler = adaptive_hmc(), - model = model, - data = x, - p_theta = p_theta, - p_x_bar_theta = p_x_bar_theta, - chains = 10, - niter = 100, - warmup = 100, - thin = 5 -) +# geweke_adaptive_hmc <- check_geweke( +# sampler = adaptive_hmc(), +# model = model, +# data = x, +# p_theta = p_theta, +# p_x_bar_theta = p_x_bar_theta, +# chains = 10, +# niter = 100, +# warmup = 100, +# thin = 5 +# ) +png("geweke-plot.png", width = 1200, height = 800, res = 150) # geweke_qq(geweke_adaptive_hmc, title = "adaptive hmc sampler Geweke test") -# -# geweke_stat_adaptive_hmc <- geweke_ks(geweke_adaptive_hmc) +sym <- c(1, 16) +pal <- c("darkorange", "purple", "cyan4") +plot(bill_dep ~ bill_len, data = penguins, pch = sym[sex], col = pal[species]) +dev.off() -geweke_stat_adaptive_hmc +# geweke_stat_adaptive_hmc <- geweke_ks(geweke_adaptive_hmc) +# +# geweke_stat_adaptive_hmc From 764393ecdcf6e2c2764761da35dbf244f7f87b45 Mon Sep 17 00:00:00 2001 From: njtierney Date: Tue, 29 Apr 2025 16:32:09 +1000 Subject: [PATCH 42/55] use v4 artifact --- .github/workflows/geweke.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/geweke.yaml b/.github/workflows/geweke.yaml index 9f29da11..a9f659bb 100644 --- a/.github/workflows/geweke.yaml +++ b/.github/workflows/geweke.yaml @@ -41,7 +41,7 @@ jobs: run: Rscript inst/greta-geweke-ahmc.R - name: Upload plot as artifact - uses: actions/upload-artifact@v3 + uses: actions/upload-artifact@v4 with: name: qq-plot path: geweke-plot.png From 89bd4bc938003c5f55026868da36ab937ff4c4d1 Mon Sep 17 00:00:00 2001 From: njtierney Date: Tue, 29 Apr 2025 16:38:03 +1000 Subject: [PATCH 43/55] update to upload the actual qq plot --- inst/greta-geweke-ahmc.R | 35 ++++++++++++++++------------------- 1 file changed, 16 insertions(+), 19 deletions(-) diff --git a/inst/greta-geweke-ahmc.R b/inst/greta-geweke-ahmc.R index f113d236..a6ff8dec 100644 --- a/inst/greta-geweke-ahmc.R +++ b/inst/greta-geweke-ahmc.R @@ -60,25 +60,22 @@ greta_theta <- normal(mu1, sd1) distribution(x) <- normal(greta_theta, sd2) model <- model(greta_theta, precision = "single") -# geweke_adaptive_hmc <- check_geweke( -# sampler = adaptive_hmc(), -# model = model, -# data = x, -# p_theta = p_theta, -# p_x_bar_theta = p_x_bar_theta, -# chains = 10, -# niter = 100, -# warmup = 100, -# thin = 5 -# ) +geweke_adaptive_hmc <- check_geweke( + sampler = adaptive_hmc(), + model = model, + data = x, + p_theta = p_theta, + p_x_bar_theta = p_x_bar_theta, + chains = 10, + niter = 100, + warmup = 100, + thin = 5 +) + +geweke_stat_adaptive_hmc <- geweke_ks(geweke_adaptive_hmc) + +geweke_stat_adaptive_hmc png("geweke-plot.png", width = 1200, height = 800, res = 150) -# geweke_qq(geweke_adaptive_hmc, title = "adaptive hmc sampler Geweke test") -sym <- c(1, 16) -pal <- c("darkorange", "purple", "cyan4") -plot(bill_dep ~ bill_len, data = penguins, pch = sym[sex], col = pal[species]) +geweke_qq(geweke_adaptive_hmc, title = "adaptive hmc sampler Geweke test") dev.off() - -# geweke_stat_adaptive_hmc <- geweke_ks(geweke_adaptive_hmc) -# -# geweke_stat_adaptive_hmc From e855010cb63003db43154335ea534522a040c2ac Mon Sep 17 00:00:00 2001 From: njtierney Date: Tue, 29 Apr 2025 18:59:19 +1000 Subject: [PATCH 44/55] up the number of samples a lot --- inst/greta-geweke-ahmc.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/inst/greta-geweke-ahmc.R b/inst/greta-geweke-ahmc.R index a6ff8dec..f21f52e1 100644 --- a/inst/greta-geweke-ahmc.R +++ b/inst/greta-geweke-ahmc.R @@ -67,8 +67,8 @@ geweke_adaptive_hmc <- check_geweke( p_theta = p_theta, p_x_bar_theta = p_x_bar_theta, chains = 10, - niter = 100, - warmup = 100, + niter = 10000, + warmup = 10000, thin = 5 ) From 03af6dc52b0ad04d258ab101e5ea6c84b62a1e2d Mon Sep 17 00:00:00 2001 From: njtierney Date: Wed, 30 Apr 2025 11:04:40 +1000 Subject: [PATCH 45/55] get geweke timeout to be 6 hours --- .github/workflows/geweke.yaml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/geweke.yaml b/.github/workflows/geweke.yaml index a9f659bb..e3282f2c 100644 --- a/.github/workflows/geweke.yaml +++ b/.github/workflows/geweke.yaml @@ -16,6 +16,7 @@ jobs: runs-on: ubuntu-latest permissions: contents: write + timeout-minutes: 360 env: GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} steps: From e1348b0c4224da16f7d0c01c06e4a09172829423 Mon Sep 17 00:00:00 2001 From: njtierney Date: Wed, 7 May 2025 14:50:31 +1000 Subject: [PATCH 46/55] try getting slice sampler to run on GHA --- .github/workflows/geweke-slice.yaml | 69 +++++++++++++++++++++++++ inst/greta-geweke-slice.R | 80 +++++++++++++++++++++++++++++ tests/testthat/helpers.R | 58 ++++++++++++++++++--- 3 files changed, 201 insertions(+), 6 deletions(-) create mode 100644 .github/workflows/geweke-slice.yaml create mode 100644 inst/greta-geweke-slice.R diff --git a/.github/workflows/geweke-slice.yaml b/.github/workflows/geweke-slice.yaml new file mode 100644 index 00000000..49830774 --- /dev/null +++ b/.github/workflows/geweke-slice.yaml @@ -0,0 +1,69 @@ +name: geweke-tests-slice + +on: + push: + branches: + - main + - master + pull_request: + branches: + - main + - master + workflow_dispatch: + +jobs: + run: + runs-on: ubuntu-latest + permissions: + contents: write + timeout-minutes: 360 + env: + GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} + steps: + - name: Checkout + uses: actions/checkout@v4 + - uses: r-lib/actions/setup-r@v2 + with: + use-public-rspm: true + - uses: r-lib/actions/setup-r-dependencies@v2 + with: + extra-packages: | + any::rcmdcheck + any::devtools + any::reticulate + any::ggplot2 + any::bayesplot + local::. + cache-version: 4 + upgrade: 'TRUE' + - name: Install greta deps + shell: Rscript {0} + run: | + library(greta) + greta::install_greta_deps(timeout = 50) + - name: Situation Report on greta install + shell: Rscript {0} + run: greta::greta_sitrep() + - name: Run Geweke tests + run: Rscript inst/greta-geweke-slice.R + + - name: Upload plot as artifact + uses: actions/upload-artifact@v4 + with: + name: qq-plots + path: geweke-*.png + retention-days: 90 # How long to keep the artifact + + - name: Session info + shell: Rscript {0} + run: | + options(width = 100) + pkgs <- installed.packages()[, "Package"] + sessioninfo::session_info(pkgs, include_base = TRUE) + + # - name: Upload plots and statistics as artifacts + # uses: actions/upload-artifact@v3 + # with: + # name: geweke-results-${{ matrix.os }} + # path: plots/ + # if-no-files-found: warn diff --git a/inst/greta-geweke-slice.R b/inst/greta-geweke-slice.R new file mode 100644 index 00000000..4ca67eac --- /dev/null +++ b/inst/greta-geweke-slice.R @@ -0,0 +1,80 @@ +devtools::load_all(".") + +n <- 10 + +mu1 <- 0 +sd1 <- 2 +sd2 <- 1 + +# prior (n draws) +p_theta <- function(n) { + rnorm(n, mu1, sd1) +} + +# likelihood +p_x_bar_theta <- function(theta) { + rnorm(n, theta, sd2) +} + +# define the greta model (single precision for slice sampler) +x <- as_data(rep(0, n)) +greta_theta <- normal(mu1, sd1) +distribution(x) <- normal(greta_theta, sd2) +model <- model(greta_theta, precision = "single") +n_iter <- 10 +n_warmup <- 10 +n_chains <- 2 +n_thin <- 1 +time_taken <- system.time({ + geweke_slice <- check_geweke( + sampler = slice(), + model = model, + data = x, + p_theta = p_theta, + p_x_bar_theta = p_x_bar_theta, + chains = n_chains, + niter = n_iter, + warmup = n_warmup, + thin = n_thin + ) +}) + +geweke_thin <- apply_thinning(geweke_slice, n_thin) + +geweke_stat_slice <- geweke_ks(geweke_thin) + +geweke_stat_slice + +the_time <- round(time_taken[3]) + +cat("this took", the_time, "seconds") + +png("geweke-plot-slice.png", width = 1200, height = 800, res = 150) + +qq_title <- build_qq_title( + "slice", + n_warmup, + n_iter, + n_chains, + n_thin, + geweke_stat_rwmh, + the_time +) +geweke_qq(geweke_thin, title = qq_title) +dev.off() + +png("geweke-mcmc-plot-slice.png", width = 1200, height = 800, res = 150) +plot(geweke_slice$draws) +dev.off() + +bplot_dens <- bayesplot::mcmc_dens_chains(geweke_slice$draws) + +ggplot2::ggsave( + filename = "geweke-gg-mcmc-plot-slice.png", + plot = bplot_dens, + width = 1200, + height = 800, + dpi = 150, + units = "px", + bg = "white" +) diff --git a/tests/testthat/helpers.R b/tests/testthat/helpers.R index 7d17c318..22dd7bfb 100644 --- a/tests/testthat/helpers.R +++ b/tests/testthat/helpers.R @@ -680,8 +680,9 @@ check_geweke <- function( ) geweke_checks <- list( - target_theta = do_thinning(target_theta, thin), - greta_theta = do_thinning(greta_theta, thin) + target_theta = target_theta, + greta_theta = greta_theta$theta, + draws = greta_theta$draws ) geweke_checks @@ -690,9 +691,13 @@ check_geweke <- function( geweke_qq <- function(geweke_checks, title) { # visualise correspondence quants <- (1:99) / 100 + # target q1 <- stats::quantile(geweke_checks$target_theta, quants) + # greta sampled q2 <- stats::quantile(geweke_checks$greta_theta, quants) - plot(q2, q1, main = title) + # q1 is target + # q2 is greta + plot(q2, q1, main = title, xlab = "Q2 greta", ylab = "Q1 target") graphics::abline(0, 1) } @@ -716,7 +721,7 @@ p_theta_greta <- function( data, p_theta, p_x_bar_theta, - # TODO note that we might want to change this to adaptive_hmc() + # this argument gets passed along based on sampler used sampler = hmc(), warmup = 1000, chains @@ -776,14 +781,17 @@ p_theta_greta <- function( verbose = FALSE ) - # trace the sample + # trace the sample - just one chain theta[i] <- tail(as.numeric(draws[[1]]), 1) } # kill the progress_bar cli::cli_progress_done() - theta + list( + theta = theta, + draws = draws + ) } # test mcmc for models with analytic posteriors @@ -929,6 +937,44 @@ do_thinning <- function(x, thinning = 1) { x[idx] } +apply_thinning <- function(geweke_result, n_thin) { + geweke_thin <- list( + target_theta = do_thinning(geweke_result$target_theta, n_thin), + greta_theta = do_thinning(geweke_result$greta_theta, n_thin) + ) + geweke_thin +} + + +build_qq_title <- function( + sampler_name, + n_warmup, + n_iter, + n_chains, + n_thin, + geweke_stat_rwmh, + the_time +) { + qq_title <- paste0( + sampler_name, + " sampler Geweke test \nwarmup = ", + n_warmup, + ", iter = ", + n_iter, + ", chains = ", + n_chains, + ", thin = ", + n_thin, + ", p-value = ", + round(geweke_stat_rwmh$p.value, 8), + ", time = ", + the_time, + "s" + ) + + qq_title +} + get_distribution_name <- function(x) { x_node <- get_node(x) From f36a26a79783d5a20bce7d71584e98ab1c9e54ab Mon Sep 17 00:00:00 2001 From: njtierney Date: Wed, 7 May 2025 14:59:13 +1000 Subject: [PATCH 47/55] use slice stat name --- inst/greta-geweke-slice.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/inst/greta-geweke-slice.R b/inst/greta-geweke-slice.R index 4ca67eac..b3d2e3b2 100644 --- a/inst/greta-geweke-slice.R +++ b/inst/greta-geweke-slice.R @@ -57,7 +57,7 @@ qq_title <- build_qq_title( n_iter, n_chains, n_thin, - geweke_stat_rwmh, + geweke_slice, the_time ) geweke_qq(geweke_thin, title = qq_title) From 968a347bcc65c377820515b776ce06c055a6dc0b Mon Sep 17 00:00:00 2001 From: njtierney Date: Wed, 7 May 2025 15:04:58 +1000 Subject: [PATCH 48/55] tweak for slice sampler --- inst/greta-geweke-slice.R | 2 +- tests/testthat/helpers.R | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/inst/greta-geweke-slice.R b/inst/greta-geweke-slice.R index b3d2e3b2..ece3150c 100644 --- a/inst/greta-geweke-slice.R +++ b/inst/greta-geweke-slice.R @@ -57,7 +57,7 @@ qq_title <- build_qq_title( n_iter, n_chains, n_thin, - geweke_slice, + geweke_stat_slice, the_time ) geweke_qq(geweke_thin, title = qq_title) diff --git a/tests/testthat/helpers.R b/tests/testthat/helpers.R index 22dd7bfb..5e6de8a1 100644 --- a/tests/testthat/helpers.R +++ b/tests/testthat/helpers.R @@ -952,7 +952,7 @@ build_qq_title <- function( n_iter, n_chains, n_thin, - geweke_stat_rwmh, + geweke_stat, the_time ) { qq_title <- paste0( @@ -966,7 +966,7 @@ build_qq_title <- function( ", thin = ", n_thin, ", p-value = ", - round(geweke_stat_rwmh$p.value, 8), + round(geweke_stat$p.value, 8), ", time = ", the_time, "s" From 811e6c0728b52fdd10787b1e4b8fefbec3621ec0 Mon Sep 17 00:00:00 2001 From: njtierney Date: Wed, 7 May 2025 15:16:53 +1000 Subject: [PATCH 49/55] add AHMC and Slice sampling to GHA --- .../{geweke.yaml => geweke-ahmc.yaml} | 18 ++- .github/workflows/geweke-slice.yaml | 4 +- inst/greta-geweke-ahmc.R | 110 +++++++++--------- inst/greta-geweke-slice.R | 13 ++- 4 files changed, 80 insertions(+), 65 deletions(-) rename .github/workflows/{geweke.yaml => geweke-ahmc.yaml} (80%) diff --git a/.github/workflows/geweke.yaml b/.github/workflows/geweke-ahmc.yaml similarity index 80% rename from .github/workflows/geweke.yaml rename to .github/workflows/geweke-ahmc.yaml index e3282f2c..461ae433 100644 --- a/.github/workflows/geweke.yaml +++ b/.github/workflows/geweke-ahmc.yaml @@ -1,4 +1,4 @@ -name: Geweke-tests +name: geweke-tests-ahmc on: push: @@ -27,7 +27,13 @@ jobs: use-public-rspm: true - uses: r-lib/actions/setup-r-dependencies@v2 with: - extra-packages: any::rcmdcheck any::devtools any::reticulate local::. + extra-packages: | + any::rcmdcheck + any::devtools + any::reticulate + any::ggplot2 + any::bayesplot + local::. cache-version: 4 upgrade: 'TRUE' - name: Install greta deps @@ -39,13 +45,13 @@ jobs: shell: Rscript {0} run: greta::greta_sitrep() - name: Run Geweke tests - run: Rscript inst/greta-geweke-ahmc.R + run: Rscript inst/greta-geweke-slice.R - - name: Upload plot as artifact + - name: Upload diagnostic plots uses: actions/upload-artifact@v4 with: - name: qq-plot - path: geweke-plot.png + name: diagnostic-plots + path: geweke-*.png retention-days: 90 # How long to keep the artifact - name: Session info diff --git a/.github/workflows/geweke-slice.yaml b/.github/workflows/geweke-slice.yaml index 49830774..0545598d 100644 --- a/.github/workflows/geweke-slice.yaml +++ b/.github/workflows/geweke-slice.yaml @@ -47,10 +47,10 @@ jobs: - name: Run Geweke tests run: Rscript inst/greta-geweke-slice.R - - name: Upload plot as artifact + - name: Upload diagnostic plots uses: actions/upload-artifact@v4 with: - name: qq-plots + name: diagnostic-plots path: geweke-*.png retention-days: 90 # How long to keep the artifact diff --git a/inst/greta-geweke-ahmc.R b/inst/greta-geweke-ahmc.R index f21f52e1..43e804e8 100644 --- a/inst/greta-geweke-ahmc.R +++ b/inst/greta-geweke-ahmc.R @@ -1,43 +1,5 @@ devtools::load_all(".") -# Set up graphics device to capture plots -# Create directory for plots if it doesn't exist -# if (!dir.exists("plots")) dir.create("plots", recursive = TRUE) -# -# script_basename <- basename(tools::file_path_sans_ext(script_path)) -# plot_filename <- paste0("plots/geweke_", script_basename, ".png") -# png(plot_filename, width = 1200, height = 800, res = 150) -# -# # Source the user's script -# tryCatch( -# { -# source(script_path) -# cat("Successfully sourced script:", script_path, "\n") -# }, -# error = function(e) { -# cat("Error occurred while running the script:\n") -# cat(as.character(e), "\n") -# quit(status = 1) -# } -# ) -# -# # Close the device to save the plot -# dev.off() -# cat("Plot saved to", plot_filename, "\n") -# -# # Save test statistics to a text file -# stats_filename <- paste0("plots/geweke_stats_", script_basename, ".txt") -# if (exists("geweke_stat_adaptive_hmc")) { -# capture.output(geweke_stat_adaptive_hmc, file = stats_filename) -# cat("Geweke test statistics saved to", stats_filename, "\n") -# } else { -# cat("Warning: No geweke_stat_adaptive_hmc object found\n") -# } -# -# # List all created files -# cat("Files in plots directory:\n") -# print(list.files("plots")) - n <- 10 mu1 <- 0 @@ -59,23 +21,65 @@ x <- as_data(rep(0, n)) greta_theta <- normal(mu1, sd1) distribution(x) <- normal(greta_theta, sd2) model <- model(greta_theta, precision = "single") +n_iter <- 10 +n_warmup <- 10 +n_chains <- 2 +n_thin <- 1 +time_taken <- system.time({ + geweke_ahmc <- check_geweke( + sampler = adaptive_hmc(), + model = model, + data = x, + p_theta = p_theta, + p_x_bar_theta = p_x_bar_theta, + chains = n_chains, + niter = n_iter, + warmup = n_warmup, + thin = n_thin + ) +}) -geweke_adaptive_hmc <- check_geweke( - sampler = adaptive_hmc(), - model = model, - data = x, - p_theta = p_theta, - p_x_bar_theta = p_x_bar_theta, - chains = 10, - niter = 10000, - warmup = 10000, - thin = 5 -) +geweke_thin <- apply_thinning(geweke_ahmc, n_thin) + +geweke_stat_ahmc <- geweke_ks(geweke_thin) + +geweke_stat_ahmc + +the_time <- round(time_taken[3]) + +cat("This took", the_time, "seconds") -geweke_stat_adaptive_hmc <- geweke_ks(geweke_adaptive_hmc) +png("geweke-qq-plot-sampler-ahmc.png", width = 1200, height = 800, res = 150) -geweke_stat_adaptive_hmc +qq_title <- build_qq_title( + "Adaptive hmc", + n_warmup, + n_iter, + n_chains, + n_thin, + geweke_stat_ahmc, + the_time +) +geweke_qq(geweke_thin, title = qq_title) +dev.off() -png("geweke-plot.png", width = 1200, height = 800, res = 150) -geweke_qq(geweke_adaptive_hmc, title = "adaptive hmc sampler Geweke test") +png( + "geweke-mcmc-coda-plot-sampler-ahmc.png", + width = 1200, + height = 800, + res = 150 +) +plot(geweke_ahmc$draws) dev.off() + +bplot_dens <- bayesplot::mcmc_dens_chains(geweke_ahmc$draws) + +ggplot2::ggsave( + filename = "geweke-dens-chains-plot-sampler-ahmc.png", + plot = bplot_dens, + width = 1200, + height = 800, + dpi = 150, + units = "px", + bg = "white" +) diff --git a/inst/greta-geweke-slice.R b/inst/greta-geweke-slice.R index ece3150c..2297e3ce 100644 --- a/inst/greta-geweke-slice.R +++ b/inst/greta-geweke-slice.R @@ -47,9 +47,9 @@ geweke_stat_slice the_time <- round(time_taken[3]) -cat("this took", the_time, "seconds") +cat("This took", the_time, "seconds") -png("geweke-plot-slice.png", width = 1200, height = 800, res = 150) +png("geweke-qq-plot-sampler-slice.png", width = 1200, height = 800, res = 150) qq_title <- build_qq_title( "slice", @@ -63,14 +63,19 @@ qq_title <- build_qq_title( geweke_qq(geweke_thin, title = qq_title) dev.off() -png("geweke-mcmc-plot-slice.png", width = 1200, height = 800, res = 150) +png( + "geweke-mcmc-coda-plot-sampler-slice.png", + width = 1200, + height = 800, + res = 150 +) plot(geweke_slice$draws) dev.off() bplot_dens <- bayesplot::mcmc_dens_chains(geweke_slice$draws) ggplot2::ggsave( - filename = "geweke-gg-mcmc-plot-slice.png", + filename = "geweke-dens-chains-plot-sampler-slice.png", plot = bplot_dens, width = 1200, height = 800, From 618215b9eb6360471380905417d30405913d2100 Mon Sep 17 00:00:00 2001 From: njtierney Date: Wed, 7 May 2025 19:09:27 +1000 Subject: [PATCH 50/55] update geweke tests to bake in the metadata into the filenames. try for all four samplers --- .github/workflows/geweke-ahmc.yaml | 1 + .github/workflows/geweke-hmc.yaml | 63 +++++++++++++++++++++ .github/workflows/geweke-rwmh.yaml | 63 +++++++++++++++++++++ .github/workflows/geweke-slice.yaml | 1 + inst/greta-geweke-ahmc.R | 45 ++++++++------- inst/greta-geweke-hmc.R | 88 +++++++++++++++++++++++++++++ inst/greta-geweke-rwmh.R | 88 +++++++++++++++++++++++++++++ inst/greta-geweke-slice.R | 53 +++++++++-------- tests/testthat/helpers.R | 61 +++++++++++++------- 9 files changed, 396 insertions(+), 67 deletions(-) create mode 100644 .github/workflows/geweke-hmc.yaml create mode 100644 .github/workflows/geweke-rwmh.yaml create mode 100644 inst/greta-geweke-hmc.R create mode 100644 inst/greta-geweke-rwmh.R diff --git a/.github/workflows/geweke-ahmc.yaml b/.github/workflows/geweke-ahmc.yaml index 461ae433..613adba3 100644 --- a/.github/workflows/geweke-ahmc.yaml +++ b/.github/workflows/geweke-ahmc.yaml @@ -33,6 +33,7 @@ jobs: any::reticulate any::ggplot2 any::bayesplot + any::glue local::. cache-version: 4 upgrade: 'TRUE' diff --git a/.github/workflows/geweke-hmc.yaml b/.github/workflows/geweke-hmc.yaml new file mode 100644 index 00000000..e3282f2c --- /dev/null +++ b/.github/workflows/geweke-hmc.yaml @@ -0,0 +1,63 @@ +name: Geweke-tests + +on: + push: + branches: + - main + - master + pull_request: + branches: + - main + - master + workflow_dispatch: + +jobs: + run: + runs-on: ubuntu-latest + permissions: + contents: write + timeout-minutes: 360 + env: + GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} + steps: + - name: Checkout + uses: actions/checkout@v4 + - uses: r-lib/actions/setup-r@v2 + with: + use-public-rspm: true + - uses: r-lib/actions/setup-r-dependencies@v2 + with: + extra-packages: any::rcmdcheck any::devtools any::reticulate local::. + cache-version: 4 + upgrade: 'TRUE' + - name: Install greta deps + shell: Rscript {0} + run: | + library(greta) + greta::install_greta_deps(timeout = 50) + - name: Situation Report on greta install + shell: Rscript {0} + run: greta::greta_sitrep() + - name: Run Geweke tests + run: Rscript inst/greta-geweke-ahmc.R + + - name: Upload plot as artifact + uses: actions/upload-artifact@v4 + with: + name: qq-plot + path: geweke-plot.png + retention-days: 90 # How long to keep the artifact + + - name: Session info + shell: Rscript {0} + run: | + options(width = 100) + pkgs <- installed.packages()[, "Package"] + sessioninfo::session_info(pkgs, include_base = TRUE) + + # - name: Upload plots and statistics as artifacts + # uses: actions/upload-artifact@v3 + # with: + # name: geweke-results-${{ matrix.os }} + # path: plots/ + # if-no-files-found: warn diff --git a/.github/workflows/geweke-rwmh.yaml b/.github/workflows/geweke-rwmh.yaml new file mode 100644 index 00000000..e3282f2c --- /dev/null +++ b/.github/workflows/geweke-rwmh.yaml @@ -0,0 +1,63 @@ +name: Geweke-tests + +on: + push: + branches: + - main + - master + pull_request: + branches: + - main + - master + workflow_dispatch: + +jobs: + run: + runs-on: ubuntu-latest + permissions: + contents: write + timeout-minutes: 360 + env: + GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} + steps: + - name: Checkout + uses: actions/checkout@v4 + - uses: r-lib/actions/setup-r@v2 + with: + use-public-rspm: true + - uses: r-lib/actions/setup-r-dependencies@v2 + with: + extra-packages: any::rcmdcheck any::devtools any::reticulate local::. + cache-version: 4 + upgrade: 'TRUE' + - name: Install greta deps + shell: Rscript {0} + run: | + library(greta) + greta::install_greta_deps(timeout = 50) + - name: Situation Report on greta install + shell: Rscript {0} + run: greta::greta_sitrep() + - name: Run Geweke tests + run: Rscript inst/greta-geweke-ahmc.R + + - name: Upload plot as artifact + uses: actions/upload-artifact@v4 + with: + name: qq-plot + path: geweke-plot.png + retention-days: 90 # How long to keep the artifact + + - name: Session info + shell: Rscript {0} + run: | + options(width = 100) + pkgs <- installed.packages()[, "Package"] + sessioninfo::session_info(pkgs, include_base = TRUE) + + # - name: Upload plots and statistics as artifacts + # uses: actions/upload-artifact@v3 + # with: + # name: geweke-results-${{ matrix.os }} + # path: plots/ + # if-no-files-found: warn diff --git a/.github/workflows/geweke-slice.yaml b/.github/workflows/geweke-slice.yaml index 0545598d..c7c96caa 100644 --- a/.github/workflows/geweke-slice.yaml +++ b/.github/workflows/geweke-slice.yaml @@ -33,6 +33,7 @@ jobs: any::reticulate any::ggplot2 any::bayesplot + any::glue local::. cache-version: 4 upgrade: 'TRUE' diff --git a/inst/greta-geweke-ahmc.R b/inst/greta-geweke-ahmc.R index 43e804e8..2e8338ce 100644 --- a/inst/greta-geweke-ahmc.R +++ b/inst/greta-geweke-ahmc.R @@ -21,13 +21,18 @@ x <- as_data(rep(0, n)) greta_theta <- normal(mu1, sd1) distribution(x) <- normal(greta_theta, sd2) model <- model(greta_theta, precision = "single") + +# ---- parameters ---- n_iter <- 10 n_warmup <- 10 n_chains <- 2 n_thin <- 1 +geweke_sampler <- adaptive_hmc() + +# ---- checking ---- time_taken <- system.time({ geweke_ahmc <- check_geweke( - sampler = adaptive_hmc(), + sampler = geweke_sampler, model = model, data = x, p_theta = p_theta, @@ -39,43 +44,41 @@ time_taken <- system.time({ ) }) +# ---- processing ---- geweke_thin <- apply_thinning(geweke_ahmc, n_thin) -geweke_stat_ahmc <- geweke_ks(geweke_thin) +geweke_stat <- geweke_ks(geweke_thin) -geweke_stat_ahmc +geweke_stat the_time <- round(time_taken[3]) cat("This took", the_time, "seconds") -png("geweke-qq-plot-sampler-ahmc.png", width = 1200, height = 800, res = 150) - -qq_title <- build_qq_title( - "Adaptive hmc", - n_warmup, - n_iter, - n_chains, - n_thin, - geweke_stat_ahmc, - the_time +# ---- plotting ---- +plot_metadata <- build_plot_metadata( + sampler = geweke_sampler, + time = the_time, + n_warmup = n_warmup, + n_iter = n_iter, + n_chains = n_chains, + n_thin = n_thin, + ks_stat = geweke_stat ) -geweke_qq(geweke_thin, title = qq_title) + +png(plot_metadata$qq_filename, width = 1200, height = 800, res = 150) + +geweke_qq(geweke_thin, title = plot_metadata$qq_title) dev.off() -png( - "geweke-mcmc-coda-plot-sampler-ahmc.png", - width = 1200, - height = 800, - res = 150 -) +png(plot_metadata$coda_filename, width = 1200, height = 800, res = 150) plot(geweke_ahmc$draws) dev.off() bplot_dens <- bayesplot::mcmc_dens_chains(geweke_ahmc$draws) ggplot2::ggsave( - filename = "geweke-dens-chains-plot-sampler-ahmc.png", + filename = plot_metadata$bayesplot_filename, plot = bplot_dens, width = 1200, height = 800, diff --git a/inst/greta-geweke-hmc.R b/inst/greta-geweke-hmc.R new file mode 100644 index 00000000..267d7d5a --- /dev/null +++ b/inst/greta-geweke-hmc.R @@ -0,0 +1,88 @@ +devtools::load_all(".") + +n <- 10 + +mu1 <- 0 +sd1 <- 2 +sd2 <- 1 + +# prior (n draws) +p_theta <- function(n) { + rnorm(n, mu1, sd1) +} + +# likelihood +p_x_bar_theta <- function(theta) { + rnorm(n, theta, sd2) +} + +# define the greta model (single precision for slice sampler) +x <- as_data(rep(0, n)) +greta_theta <- normal(mu1, sd1) +distribution(x) <- normal(greta_theta, sd2) +model <- model(greta_theta, precision = "single") + +# ---- parameters ---- +n_iter <- 10 +n_warmup <- 10 +n_chains <- 2 +n_thin <- 1 +geweke_sampler <- hmc() + +# ---- checking ---- +time_taken <- system.time({ + geweke_ahmc <- check_geweke( + sampler = geweke_sampler, + model = model, + data = x, + p_theta = p_theta, + p_x_bar_theta = p_x_bar_theta, + chains = n_chains, + niter = n_iter, + warmup = n_warmup, + thin = n_thin + ) +}) + +# ---- processing ---- +geweke_thin <- apply_thinning(geweke_ahmc, n_thin) + +geweke_stat <- geweke_ks(geweke_thin) + +geweke_stat + +the_time <- round(time_taken[3]) + +cat("This took", the_time, "seconds") + +# ---- plotting ---- +plot_metadata <- build_plot_metadata( + sampler = geweke_sampler, + time = the_time, + n_warmup = n_warmup, + n_iter = n_iter, + n_chains = n_chains, + n_thin = n_thin, + ks_stat = geweke_stat +) + +png(plot_metadata$qq_filename, width = 1200, height = 800, res = 150) + +geweke_qq(geweke_thin, title = plot_metadata$qq_title) +dev.off() + +png(plot_metadata$coda_filename, width = 1200, height = 800, res = 150) +plot(geweke_ahmc$draws) +dev.off() + +bplot_dens <- bayesplot::mcmc_dens_chains(geweke_ahmc$draws) + +ggplot2::ggsave( + filename = plot_metadata$bayesplot_filename, + plot = bplot_dens, + width = 1200, + height = 800, + dpi = 150, + units = "px", + bg = "white" +) diff --git a/inst/greta-geweke-rwmh.R b/inst/greta-geweke-rwmh.R new file mode 100644 index 00000000..7f5a2121 --- /dev/null +++ b/inst/greta-geweke-rwmh.R @@ -0,0 +1,88 @@ +devtools::load_all(".") + +n <- 10 + +mu1 <- 0 +sd1 <- 2 +sd2 <- 1 + +# prior (n draws) +p_theta <- function(n) { + rnorm(n, mu1, sd1) +} + +# likelihood +p_x_bar_theta <- function(theta) { + rnorm(n, theta, sd2) +} + +# define the greta model (single precision for slice sampler) +x <- as_data(rep(0, n)) +greta_theta <- normal(mu1, sd1) +distribution(x) <- normal(greta_theta, sd2) +model <- model(greta_theta, precision = "single") + +# ---- parameters ---- +n_iter <- 10 +n_warmup <- 10 +n_chains <- 2 +n_thin <- 1 +geweke_sampler <- rwmh() + +# ---- checking ---- +time_taken <- system.time({ + geweke_ahmc <- check_geweke( + sampler = geweke_sampler, + model = model, + data = x, + p_theta = p_theta, + p_x_bar_theta = p_x_bar_theta, + chains = n_chains, + niter = n_iter, + warmup = n_warmup, + thin = n_thin + ) +}) + +# ---- processing ---- +geweke_thin <- apply_thinning(geweke_ahmc, n_thin) + +geweke_stat <- geweke_ks(geweke_thin) + +geweke_stat + +the_time <- round(time_taken[3]) + +cat("This took", the_time, "seconds") + +# ---- plotting ---- +plot_metadata <- build_plot_metadata( + sampler = geweke_sampler, + time = the_time, + n_warmup = n_warmup, + n_iter = n_iter, + n_chains = n_chains, + n_thin = n_thin, + ks_stat = geweke_stat +) + +png(plot_metadata$qq_filename, width = 1200, height = 800, res = 150) + +geweke_qq(geweke_thin, title = plot_metadata$qq_title) +dev.off() + +png(plot_metadata$coda_filename, width = 1200, height = 800, res = 150) +plot(geweke_ahmc$draws) +dev.off() + +bplot_dens <- bayesplot::mcmc_dens_chains(geweke_ahmc$draws) + +ggplot2::ggsave( + filename = plot_metadata$bayesplot_filename, + plot = bplot_dens, + width = 1200, + height = 800, + dpi = 150, + units = "px", + bg = "white" +) diff --git a/inst/greta-geweke-slice.R b/inst/greta-geweke-slice.R index 2297e3ce..197a2019 100644 --- a/inst/greta-geweke-slice.R +++ b/inst/greta-geweke-slice.R @@ -21,13 +21,18 @@ x <- as_data(rep(0, n)) greta_theta <- normal(mu1, sd1) distribution(x) <- normal(greta_theta, sd2) model <- model(greta_theta, precision = "single") + +# ---- parameters ---- n_iter <- 10 n_warmup <- 10 n_chains <- 2 n_thin <- 1 +geweke_sampler <- slice() + +# ---- checking ---- time_taken <- system.time({ - geweke_slice <- check_geweke( - sampler = slice(), + geweke_ahmc <- check_geweke( + sampler = geweke_sampler, model = model, data = x, p_theta = p_theta, @@ -39,43 +44,41 @@ time_taken <- system.time({ ) }) -geweke_thin <- apply_thinning(geweke_slice, n_thin) +# ---- processing ---- +geweke_thin <- apply_thinning(geweke_ahmc, n_thin) -geweke_stat_slice <- geweke_ks(geweke_thin) +geweke_stat <- geweke_ks(geweke_thin) -geweke_stat_slice +geweke_stat the_time <- round(time_taken[3]) cat("This took", the_time, "seconds") -png("geweke-qq-plot-sampler-slice.png", width = 1200, height = 800, res = 150) - -qq_title <- build_qq_title( - "slice", - n_warmup, - n_iter, - n_chains, - n_thin, - geweke_stat_slice, - the_time +# ---- plotting ---- +plot_metadata <- build_plot_metadata( + sampler = geweke_sampler, + time = the_time, + n_warmup = n_warmup, + n_iter = n_iter, + n_chains = n_chains, + n_thin = n_thin, + ks_stat = geweke_stat ) -geweke_qq(geweke_thin, title = qq_title) + +png(plot_metadata$qq_filename, width = 1200, height = 800, res = 150) + +geweke_qq(geweke_thin, title = plot_metadata$qq_title) dev.off() -png( - "geweke-mcmc-coda-plot-sampler-slice.png", - width = 1200, - height = 800, - res = 150 -) -plot(geweke_slice$draws) +png(plot_metadata$coda_filename, width = 1200, height = 800, res = 150) +plot(geweke_ahmc$draws) dev.off() -bplot_dens <- bayesplot::mcmc_dens_chains(geweke_slice$draws) +bplot_dens <- bayesplot::mcmc_dens_chains(geweke_ahmc$draws) ggplot2::ggsave( - filename = "geweke-dens-chains-plot-sampler-slice.png", + filename = plot_metadata$bayesplot_filename, plot = bplot_dens, width = 1200, height = 800, diff --git a/tests/testthat/helpers.R b/tests/testthat/helpers.R index 5e6de8a1..cfc75fb6 100644 --- a/tests/testthat/helpers.R +++ b/tests/testthat/helpers.R @@ -945,34 +945,53 @@ apply_thinning <- function(geweke_result, n_thin) { geweke_thin } - -build_qq_title <- function( - sampler_name, +build_plot_metadata <- function( + sampler, + time, n_warmup, n_iter, n_chains, n_thin, - geweke_stat, - the_time + ks_stat ) { - qq_title <- paste0( - sampler_name, - " sampler Geweke test \nwarmup = ", - n_warmup, - ", iter = ", - n_iter, - ", chains = ", - n_chains, - ", thin = ", - n_thin, - ", p-value = ", - round(geweke_stat$p.value, 8), - ", time = ", - the_time, - "s" + sampler_class <- class(sampler) + + sampler <- switch( + sampler_class[1], + "adaptive_hmc_sampler" = "ahmc", + "slice sampler" = "slice", + "rwmh sampler" = "rwmh", + "hmc sampler" = "hmc" + ) + + qq_title <- glue::glue( + "Sampler {sampler} Geweke test (time taken = {time}s)\\ + \nwarmup = {n_warmup}, iter = {n_iter}, chains = {n_chains}, thin = {n_thin}\np-value = {round(ks_stat$p.value, 8)}" + ) + + param_stub <- glue::glue( + "warm-{n_warmup}-it-{n_iter}\\ + -chain-{n_chains}-thin-{n_thin}.png" + ) + + qq_filename <- glue::glue( + "geweke-{sampler}-qq-{param_stub}" + ) + + coda_filename <- glue::glue( + "geweke-{sampler}-coda-{param_stub}" ) - qq_title + bayesplot_filename <- glue::glue( + "geweke-{sampler}-dens-{param_stub}" + ) + + list( + qq_title = qq_title, + qq_filename = qq_filename, + coda_filename = coda_filename, + bayesplot_filename = bayesplot_filename + ) } From 3e6b7872edb83c2689c1be21d0b2056209c7782d Mon Sep 17 00:00:00 2001 From: njtierney Date: Wed, 7 May 2025 19:19:37 +1000 Subject: [PATCH 51/55] make sure GHA run the correct tests, and give the artifacts a better name --- .github/workflows/geweke-ahmc.yaml | 4 ++-- .github/workflows/geweke-hmc.yaml | 8 ++++---- .github/workflows/geweke-rwmh.yaml | 8 ++++---- .github/workflows/geweke-slice.yaml | 2 +- inst/greta-geweke-ahmc.R | 16 ++++++++-------- inst/greta-geweke-hmc.R | 16 ++++++++-------- inst/greta-geweke-rwmh.R | 16 ++++++++-------- inst/greta-geweke-slice.R | 16 ++++++++-------- 8 files changed, 43 insertions(+), 43 deletions(-) diff --git a/.github/workflows/geweke-ahmc.yaml b/.github/workflows/geweke-ahmc.yaml index 613adba3..c16b226a 100644 --- a/.github/workflows/geweke-ahmc.yaml +++ b/.github/workflows/geweke-ahmc.yaml @@ -46,12 +46,12 @@ jobs: shell: Rscript {0} run: greta::greta_sitrep() - name: Run Geweke tests - run: Rscript inst/greta-geweke-slice.R + run: Rscript inst/greta-geweke-ahmc.R - name: Upload diagnostic plots uses: actions/upload-artifact@v4 with: - name: diagnostic-plots + name: geweke-ahmc-diagnostic-plots path: geweke-*.png retention-days: 90 # How long to keep the artifact diff --git a/.github/workflows/geweke-hmc.yaml b/.github/workflows/geweke-hmc.yaml index e3282f2c..c6e62160 100644 --- a/.github/workflows/geweke-hmc.yaml +++ b/.github/workflows/geweke-hmc.yaml @@ -1,4 +1,4 @@ -name: Geweke-tests +name: geweke-tests-hmc on: push: @@ -39,13 +39,13 @@ jobs: shell: Rscript {0} run: greta::greta_sitrep() - name: Run Geweke tests - run: Rscript inst/greta-geweke-ahmc.R + run: Rscript inst/greta-geweke-hmc.R - name: Upload plot as artifact uses: actions/upload-artifact@v4 with: - name: qq-plot - path: geweke-plot.png + name: geweke-hmc-diagnostic-plots + path: geweke-*.png retention-days: 90 # How long to keep the artifact - name: Session info diff --git a/.github/workflows/geweke-rwmh.yaml b/.github/workflows/geweke-rwmh.yaml index e3282f2c..7d954725 100644 --- a/.github/workflows/geweke-rwmh.yaml +++ b/.github/workflows/geweke-rwmh.yaml @@ -1,4 +1,4 @@ -name: Geweke-tests +name: geweke-tests-rwmh on: push: @@ -39,13 +39,13 @@ jobs: shell: Rscript {0} run: greta::greta_sitrep() - name: Run Geweke tests - run: Rscript inst/greta-geweke-ahmc.R + run: Rscript inst/greta-geweke-rwmh.R - name: Upload plot as artifact uses: actions/upload-artifact@v4 with: - name: qq-plot - path: geweke-plot.png + name: geweke-rwmh-diagnostic-plots + path: geweke-*.png retention-days: 90 # How long to keep the artifact - name: Session info diff --git a/.github/workflows/geweke-slice.yaml b/.github/workflows/geweke-slice.yaml index c7c96caa..5800cb83 100644 --- a/.github/workflows/geweke-slice.yaml +++ b/.github/workflows/geweke-slice.yaml @@ -51,7 +51,7 @@ jobs: - name: Upload diagnostic plots uses: actions/upload-artifact@v4 with: - name: diagnostic-plots + name: geweke-slice-diagnostic-plots path: geweke-*.png retention-days: 90 # How long to keep the artifact diff --git a/inst/greta-geweke-ahmc.R b/inst/greta-geweke-ahmc.R index 2e8338ce..3ac0dc22 100644 --- a/inst/greta-geweke-ahmc.R +++ b/inst/greta-geweke-ahmc.R @@ -1,35 +1,35 @@ devtools::load_all(".") +# parameters ---- n <- 10 - mu1 <- 0 sd1 <- 2 sd2 <- 1 -# prior (n draws) +# prior (n draws) ---- p_theta <- function(n) { rnorm(n, mu1, sd1) } -# likelihood +# likelihood ---- p_x_bar_theta <- function(theta) { rnorm(n, theta, sd2) } -# define the greta model (single precision for slice sampler) +# define the greta model ---- x <- as_data(rep(0, n)) greta_theta <- normal(mu1, sd1) distribution(x) <- normal(greta_theta, sd2) -model <- model(greta_theta, precision = "single") +model <- model(greta_theta) -# ---- parameters ---- +# mcmc parameters ---- n_iter <- 10 n_warmup <- 10 n_chains <- 2 n_thin <- 1 geweke_sampler <- adaptive_hmc() -# ---- checking ---- +# checking ---- time_taken <- system.time({ geweke_ahmc <- check_geweke( sampler = geweke_sampler, @@ -44,7 +44,7 @@ time_taken <- system.time({ ) }) -# ---- processing ---- +# processing ---- geweke_thin <- apply_thinning(geweke_ahmc, n_thin) geweke_stat <- geweke_ks(geweke_thin) diff --git a/inst/greta-geweke-hmc.R b/inst/greta-geweke-hmc.R index 267d7d5a..1cc84c1d 100644 --- a/inst/greta-geweke-hmc.R +++ b/inst/greta-geweke-hmc.R @@ -1,35 +1,35 @@ devtools::load_all(".") +# parameters ---- n <- 10 - mu1 <- 0 sd1 <- 2 sd2 <- 1 -# prior (n draws) +# prior (n draws) ---- p_theta <- function(n) { rnorm(n, mu1, sd1) } -# likelihood +# likelihood ---- p_x_bar_theta <- function(theta) { rnorm(n, theta, sd2) } -# define the greta model (single precision for slice sampler) +# define the greta model ---- x <- as_data(rep(0, n)) greta_theta <- normal(mu1, sd1) distribution(x) <- normal(greta_theta, sd2) -model <- model(greta_theta, precision = "single") +model <- model(greta_theta) -# ---- parameters ---- +# mcmc parameters ---- n_iter <- 10 n_warmup <- 10 n_chains <- 2 n_thin <- 1 geweke_sampler <- hmc() -# ---- checking ---- +# checking ---- time_taken <- system.time({ geweke_ahmc <- check_geweke( sampler = geweke_sampler, @@ -44,7 +44,7 @@ time_taken <- system.time({ ) }) -# ---- processing ---- +# processing ---- geweke_thin <- apply_thinning(geweke_ahmc, n_thin) geweke_stat <- geweke_ks(geweke_thin) diff --git a/inst/greta-geweke-rwmh.R b/inst/greta-geweke-rwmh.R index 7f5a2121..fc5dc1dd 100644 --- a/inst/greta-geweke-rwmh.R +++ b/inst/greta-geweke-rwmh.R @@ -1,35 +1,35 @@ devtools::load_all(".") +# parameters ---- n <- 10 - mu1 <- 0 sd1 <- 2 sd2 <- 1 -# prior (n draws) +# prior (n draws) ---- p_theta <- function(n) { rnorm(n, mu1, sd1) } -# likelihood +# likelihood ---- p_x_bar_theta <- function(theta) { rnorm(n, theta, sd2) } -# define the greta model (single precision for slice sampler) +# define the greta model ---- x <- as_data(rep(0, n)) greta_theta <- normal(mu1, sd1) distribution(x) <- normal(greta_theta, sd2) -model <- model(greta_theta, precision = "single") +model <- model(greta_theta) -# ---- parameters ---- +# mcmc parameters ---- n_iter <- 10 n_warmup <- 10 n_chains <- 2 n_thin <- 1 geweke_sampler <- rwmh() -# ---- checking ---- +# checking ---- time_taken <- system.time({ geweke_ahmc <- check_geweke( sampler = geweke_sampler, @@ -44,7 +44,7 @@ time_taken <- system.time({ ) }) -# ---- processing ---- +# processing ---- geweke_thin <- apply_thinning(geweke_ahmc, n_thin) geweke_stat <- geweke_ks(geweke_thin) diff --git a/inst/greta-geweke-slice.R b/inst/greta-geweke-slice.R index 197a2019..f0556953 100644 --- a/inst/greta-geweke-slice.R +++ b/inst/greta-geweke-slice.R @@ -1,35 +1,35 @@ devtools::load_all(".") +# parameters ---- n <- 10 - mu1 <- 0 sd1 <- 2 sd2 <- 1 -# prior (n draws) +# prior (n draws) ---- p_theta <- function(n) { rnorm(n, mu1, sd1) } -# likelihood +# likelihood ---- p_x_bar_theta <- function(theta) { rnorm(n, theta, sd2) } -# define the greta model (single precision for slice sampler) +# define the greta model ---- x <- as_data(rep(0, n)) greta_theta <- normal(mu1, sd1) distribution(x) <- normal(greta_theta, sd2) -model <- model(greta_theta, precision = "single") +model <- model(greta_theta) -# ---- parameters ---- +# mcmc parameters ---- n_iter <- 10 n_warmup <- 10 n_chains <- 2 n_thin <- 1 geweke_sampler <- slice() -# ---- checking ---- +# checking ---- time_taken <- system.time({ geweke_ahmc <- check_geweke( sampler = geweke_sampler, @@ -44,7 +44,7 @@ time_taken <- system.time({ ) }) -# ---- processing ---- +# processing ---- geweke_thin <- apply_thinning(geweke_ahmc, n_thin) geweke_stat <- geweke_ks(geweke_thin) From bbf97f205608896a446f54009a4a861e731e2eeb Mon Sep 17 00:00:00 2001 From: njtierney Date: Wed, 7 May 2025 19:39:41 +1000 Subject: [PATCH 52/55] try running each of the samplers for a similar amount of time except for ahmc --- inst/greta-geweke-ahmc.R | 8 ++++---- inst/greta-geweke-hmc.R | 6 +++--- inst/greta-geweke-rwmh.R | 6 +++--- inst/greta-geweke-slice.R | 4 ++-- 4 files changed, 12 insertions(+), 12 deletions(-) diff --git a/inst/greta-geweke-ahmc.R b/inst/greta-geweke-ahmc.R index 3ac0dc22..c5f187e1 100644 --- a/inst/greta-geweke-ahmc.R +++ b/inst/greta-geweke-ahmc.R @@ -23,10 +23,10 @@ distribution(x) <- normal(greta_theta, sd2) model <- model(greta_theta) # mcmc parameters ---- -n_iter <- 10 -n_warmup <- 10 -n_chains <- 2 -n_thin <- 1 +n_iter <- 500 +n_warmup <- 10000 +n_chains <- 20 +n_thin <- 5 geweke_sampler <- adaptive_hmc() # checking ---- diff --git a/inst/greta-geweke-hmc.R b/inst/greta-geweke-hmc.R index 1cc84c1d..54d15578 100644 --- a/inst/greta-geweke-hmc.R +++ b/inst/greta-geweke-hmc.R @@ -23,10 +23,10 @@ distribution(x) <- normal(greta_theta, sd2) model <- model(greta_theta) # mcmc parameters ---- -n_iter <- 10 -n_warmup <- 10 +n_iter <- 2000 +n_warmup <- 2000 n_chains <- 2 -n_thin <- 1 +n_thin <- 5 geweke_sampler <- hmc() # checking ---- diff --git a/inst/greta-geweke-rwmh.R b/inst/greta-geweke-rwmh.R index fc5dc1dd..c4d41bee 100644 --- a/inst/greta-geweke-rwmh.R +++ b/inst/greta-geweke-rwmh.R @@ -23,10 +23,10 @@ distribution(x) <- normal(greta_theta, sd2) model <- model(greta_theta) # mcmc parameters ---- -n_iter <- 10 -n_warmup <- 10 +n_iter <- 2000 +n_warmup <- 2000 n_chains <- 2 -n_thin <- 1 +n_thin <- 5 geweke_sampler <- rwmh() # checking ---- diff --git a/inst/greta-geweke-slice.R b/inst/greta-geweke-slice.R index f0556953..d40caf87 100644 --- a/inst/greta-geweke-slice.R +++ b/inst/greta-geweke-slice.R @@ -23,8 +23,8 @@ distribution(x) <- normal(greta_theta, sd2) model <- model(greta_theta) # mcmc parameters ---- -n_iter <- 10 -n_warmup <- 10 +n_iter <- 2000 +n_warmup <- 2000 n_chains <- 2 n_thin <- 1 geweke_sampler <- slice() From 3cbc13b2eae19e1fa4931b82adfa9522b386ce4f Mon Sep 17 00:00:00 2001 From: njtierney Date: Wed, 7 May 2025 20:25:32 +1000 Subject: [PATCH 53/55] try with 250 - half as many iterations --- inst/greta-geweke-ahmc.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/inst/greta-geweke-ahmc.R b/inst/greta-geweke-ahmc.R index c5f187e1..a23ea184 100644 --- a/inst/greta-geweke-ahmc.R +++ b/inst/greta-geweke-ahmc.R @@ -23,7 +23,7 @@ distribution(x) <- normal(greta_theta, sd2) model <- model(greta_theta) # mcmc parameters ---- -n_iter <- 500 +n_iter <- 250 n_warmup <- 10000 n_chains <- 20 n_thin <- 5 From 925cbcf7ceb24364950abd0ce843f53032fae6fc Mon Sep 17 00:00:00 2001 From: njtierney Date: Wed, 7 May 2025 20:36:42 +1000 Subject: [PATCH 54/55] only run these runners when their relevant files change...which is now, but still --- .github/workflows/R-CMD-check.yaml | 9 +++++++++ .github/workflows/geweke-ahmc.yaml | 8 ++++++++ .github/workflows/geweke-hmc.yaml | 8 ++++++++ .github/workflows/geweke-rwmh.yaml | 8 ++++++++ .github/workflows/geweke-slice.yaml | 8 ++++++++ .github/workflows/pkgdown.yaml | 3 +++ .github/workflows/test-coverage.yaml | 16 ++++++++++++++-- 7 files changed, 58 insertions(+), 2 deletions(-) diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml index d0cbcc13..5d7ceecb 100644 --- a/.github/workflows/R-CMD-check.yaml +++ b/.github/workflows/R-CMD-check.yaml @@ -3,10 +3,19 @@ on: branches: - main - master + paths: + - '.github/workflows/R-CMD-check.yaml' # Run when the workflow file changes + - 'R/**' # Run when R code changes + - 'tests/**' # Run when test code changes pull_request: branches: - main - master + paths: + - '.github/workflows/R-CMD-check.yaml' # Run when the workflow file changes + - 'R/**' + - 'tests/**' + workflow_dispatch: schedule: - cron: '1 23 * * Sun' diff --git a/.github/workflows/geweke-ahmc.yaml b/.github/workflows/geweke-ahmc.yaml index c16b226a..1bef2ac2 100644 --- a/.github/workflows/geweke-ahmc.yaml +++ b/.github/workflows/geweke-ahmc.yaml @@ -5,10 +5,18 @@ on: branches: - main - master + paths: + - '.github/workflows/geweke-ahmc.yaml' # Run when the workflow file changes + - 'R/**' # Run when R code changes + - 'inst/greta-geweke-ahmc.R' # Run when the test script changes pull_request: branches: - main - master + paths: + - '.github/workflows/geweke-ahmc.yaml' + - 'R/**' + - 'inst/greta-geweke-ahmc.R' workflow_dispatch: jobs: diff --git a/.github/workflows/geweke-hmc.yaml b/.github/workflows/geweke-hmc.yaml index c6e62160..e08df617 100644 --- a/.github/workflows/geweke-hmc.yaml +++ b/.github/workflows/geweke-hmc.yaml @@ -5,10 +5,18 @@ on: branches: - main - master + paths: + - '.github/workflows/geweke-hmc.yaml' # Run when the workflow file changes + - 'R/**' # Run when R code changes + - 'inst/greta-geweke-hmc.R' # Run when the test script changes pull_request: branches: - main - master + paths: + - '.github/workflows/geweke-hmc.yaml' + - 'R/**' + - 'inst/greta-geweke-hmc.R' workflow_dispatch: jobs: diff --git a/.github/workflows/geweke-rwmh.yaml b/.github/workflows/geweke-rwmh.yaml index 7d954725..f6f8d5c0 100644 --- a/.github/workflows/geweke-rwmh.yaml +++ b/.github/workflows/geweke-rwmh.yaml @@ -5,10 +5,18 @@ on: branches: - main - master + paths: + - '.github/workflows/geweke-rwmh.yaml' # Run when the workflow file changes + - 'R/**' # Run when R code changes + - 'inst/greta-geweke-rwmh.R' # Run when the test script changes pull_request: branches: - main - master + paths: + - '.github/workflows/geweke-rwmh.yaml' + - 'R/**' + - 'inst/greta-geweke-rwmh.R' workflow_dispatch: jobs: diff --git a/.github/workflows/geweke-slice.yaml b/.github/workflows/geweke-slice.yaml index 5800cb83..2fd2b99c 100644 --- a/.github/workflows/geweke-slice.yaml +++ b/.github/workflows/geweke-slice.yaml @@ -5,10 +5,18 @@ on: branches: - main - master + paths: + - '.github/workflows/geweke-slice.yaml' # Run when the workflow file changes + - 'R/**' # Run when R code changes + - 'inst/greta-geweke-slice.R' # Run when the test script changes pull_request: branches: - main - master + paths: + - '.github/workflows/geweke-slice.yaml' + - 'R/**' + - 'inst/greta-geweke-slice.R' workflow_dispatch: jobs: diff --git a/.github/workflows/pkgdown.yaml b/.github/workflows/pkgdown.yaml index 2468d333..1c6253b4 100644 --- a/.github/workflows/pkgdown.yaml +++ b/.github/workflows/pkgdown.yaml @@ -3,6 +3,9 @@ on: branches: - main - master + paths: + - 'R/**' # Run when R code changes + - '.github/workflows/pkgdown.yaml' # Run when the workflow file changes name: pkgdown diff --git a/.github/workflows/test-coverage.yaml b/.github/workflows/test-coverage.yaml index 558fef72..3db183f3 100644 --- a/.github/workflows/test-coverage.yaml +++ b/.github/workflows/test-coverage.yaml @@ -2,9 +2,21 @@ # Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help on: push: - branches: [main, master] + branches: + - main + - master + paths: + - '.github/workflows/test-coverage.yaml' # Run when the workflow file changes + - 'R/**' # Run when R code changes + - 'tests/**' # Run when test code changes pull_request: - branches: [main, master] + branches: + - main + - master + paths: + - '.github/workflows/test-coverage.yaml' + - 'R/**' + - 'tests/**' name: test-coverage.yaml From 4cc989fe4b4eb0745a19fcd5bd42112fce99baa8 Mon Sep 17 00:00:00 2001 From: njtierney Date: Tue, 13 May 2025 11:57:44 +1000 Subject: [PATCH 55/55] use 10 chains for HMC geweke test --- inst/greta-geweke-hmc.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/inst/greta-geweke-hmc.R b/inst/greta-geweke-hmc.R index 54d15578..7a577b8c 100644 --- a/inst/greta-geweke-hmc.R +++ b/inst/greta-geweke-hmc.R @@ -25,7 +25,7 @@ model <- model(greta_theta) # mcmc parameters ---- n_iter <- 2000 n_warmup <- 2000 -n_chains <- 2 +n_chains <- 10 n_thin <- 5 geweke_sampler <- hmc()