Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
187 changes: 108 additions & 79 deletions R/args.R
Original file line number Diff line number Diff line change
Expand Up @@ -81,7 +81,7 @@ CmdStanArgs <- R6::R6Class(
init <- process_init(init, num_inits, model_variables)
self$init <- init
self$opencl_ids <- opencl_ids
self$num_threads = NULL
self$num_threads <- NULL
self$method_args$validate(num_procs = length(self$proc_ids))
if (is.logical(self$save_cmdstan_config)) {
self$save_cmdstan_config <- as.integer(self$save_cmdstan_config)
Expand Down Expand Up @@ -984,20 +984,20 @@ validate_pathfinder_args <- function(self) {
self$num_elbo_draws <- as.integer(self$num_elbo_draws)
}
if (!is.null(self$save_single_paths) && is.logical(self$save_single_paths)) {
self$save_single_paths = as.integer(self$save_single_paths)
self$save_single_paths <- as.integer(self$save_single_paths)
}
checkmate::assert_integerish(self$save_single_paths, null.ok = TRUE,
lower = 0, upper = 1, len = 1)
if (!is.null(self$save_single_paths)) {
self$save_single_paths <- 0
}
if (!is.null(self$psis_resample) && is.logical(self$psis_resample)) {
self$psis_resample = as.integer(self$psis_resample)
self$psis_resample <- as.integer(self$psis_resample)
}
checkmate::assert_integerish(self$psis_resample, null.ok = TRUE,
lower = 0, upper = 1, len = 1)
if (!is.null(self$calculate_lp) && is.logical(self$calculate_lp)) {
self$calculate_lp = as.integer(self$calculate_lp)
self$calculate_lp <- as.integer(self$calculate_lp)
}
checkmate::assert_integerish(self$calculate_lp, null.ok = TRUE,
lower = 0, upper = 1, len = 1)
Expand All @@ -1018,23 +1018,45 @@ validate_pathfinder_args <- function(self) {
}


# Validation helpers ------------------------------------------------------

#' Validate exe file exists
# Init helpers ------------------------------------------------------------

#' Build a model_variables structure from a draws object
#'
#' When a model has been created without a Stan file,
#' `model$variables()` is unavailable. This helper infers parameter
#' names and dimensions from `posterior::variables(as_draws_df(...))`.
#' In that representation containers are expanded (e.g. `beta[1]`,
#' `gamma[1,2]`) while scalars appear as bare names (e.g. `sigma`).
#' The number of dimensions is inferred from the index pattern:
#' `mu[1]` has 1, `mu[1,2]` has 2, etc.
#'
#' @noRd
#' @param exe_file Path to executable.
#' @return Either throws an error or returns `invisible(TRUE)`
validate_exe_file <- function(exe_file) {
if (!length(exe_file) ||
!nzchar(exe_file) ||
!file.exists(exe_file)) {
stop("Model not compiled. Try running the compile() method first.",
call. = FALSE)
}
invisible(TRUE)
#' @param draws A draws object (any format supported by posterior).
#' @return A list with a `parameters` element in the same format as
#' `model$variables()`.
model_variables_from_draws <- function(draws) {
df_vars <- posterior::variables(posterior::as_draws_df(draws))
df_vars <- df_vars[!grepl("__$", df_vars)]
has_bracket <- grepl("\\[", df_vars)
scalars <- df_vars[!has_bracket]
# For containers, extract base name and count dimensions from the
# index pattern of the first occurrence (e.g. "mu[1,2]" -> 2 dims)
container_vars <- df_vars[has_bracket]
container_names <- sub("\\[.*", "", container_vars)
container_indices <- sub("^[^\\[]*\\[(.*)\\]$", "\\1", container_vars)
parameters <- list()
for (var_name in scalars) {
parameters[[var_name]] <- list(type = "real", dimensions = 0L)
}
for (var_name in unique(container_names)) {
idx <- match(var_name, container_names)
ndims <- length(strsplit(container_indices[idx], ",")[[1]])
parameters[[var_name]] <- list(type = "real", dimensions = ndims)
}
list(parameters = parameters)
}


#' Generic for processing inits
#' @noRd
process_init <- function(init, ...) {
Expand Down Expand Up @@ -1080,12 +1102,11 @@ process_init.default <- function(init, ...) {
process_init.draws <- function(init, num_procs, model_variables = NULL,
warn_partial = getOption("cmdstanr_warn_inits", TRUE),
...) {
if (!is.null(model_variables)) {
variable_names = names(model_variables$parameters)
} else {
variable_names = colnames(draws)[!grepl("__", colnames(draws))]
}
draws <- posterior::as_draws_df(init)
if (is.null(model_variables)) {
model_variables <- model_variables_from_draws(draws)
}
variable_names <- names(model_variables$parameters)
# Since all other process_init functions return `num_proc` inits
# This will only happen if a raw draws object is passed
if (nrow(draws) < num_procs) {
Expand All @@ -1095,38 +1116,38 @@ process_init.draws <- function(init, num_procs, model_variables = NULL,
draws <- posterior::resample_draws(draws, ndraws = num_procs,
method ="simple_no_replace")
}
draws_rvar = posterior::as_draws_rvars(draws)
draws_rvar <- posterior::as_draws_rvars(draws)
variable_names <- variable_names[variable_names %in% names(draws_rvar)]
draws_rvar <- posterior::subset_draws(draws_rvar, variable = variable_names)
inits = lapply(1:num_procs, function(draw_iter) {
init_i = lapply(variable_names, function(var_name) {
x = .remove_leftmost_dim(posterior::draws_of(
inits <- lapply(1:num_procs, function(draw_iter) {
init_i <- lapply(variable_names, function(var_name) {
x <- .remove_leftmost_dim(posterior::draws_of(
posterior::subset_draws(draws_rvar[[var_name]], draw=draw_iter)))
if (model_variables$parameters[[var_name]]$dimensions == 0) {
return(as.double(x))
} else {
return(x)
}
})
bad_names = unlist(lapply(variable_names, function(var_name) {
x = drop(posterior::draws_of(drop(
bad_names <- unlist(lapply(variable_names, function(var_name) {
x <- drop(posterior::draws_of(drop(
posterior::subset_draws(draws_rvar[[var_name]], draw=draw_iter))))
if (any(is.infinite(x)) || any(is.na(x))) {
return(var_name)
}
return("")
}))
any_na_or_inf = bad_names != ""
any_na_or_inf <- bad_names != ""
if (any(any_na_or_inf)) {
err_msg = paste0(paste(bad_names[any_na_or_inf], collapse = ", "), " contains NA or Inf values!")
err_msg <- paste0(paste(bad_names[any_na_or_inf], collapse = ", "), " contains NA or Inf values!")
if (length(any_na_or_inf) > 1) {
err_msg = paste0("Variables: ", err_msg)
err_msg <- paste0("Variables: ", err_msg)
} else {
err_msg = paste0("Variable: ", err_msg)
err_msg <- paste0("Variable: ", err_msg)
}
stop(err_msg)
}
names(init_i) = variable_names
names(init_i) <- variable_names
return(init_i)
})
return(process_init(inits, num_procs, model_variables, warn_partial))
Expand Down Expand Up @@ -1241,8 +1262,7 @@ process_init.function <- function(init, num_procs, model_variables = NULL,

#' Validate a fit is a valid init
#' @noRd
validate_fit_init = function(init, model_variables) {
# Convert from data.table to data.frame
validate_fit_init <- function(init, model_variables) {
if (all(init$return_codes() == 1)) {
stop("We are unable to create initial values from a model with no samples. Please check the results of the model used for inits before continuing.")
} else if (!is.null(model_variables) &&!any(names(model_variables$parameters) %in% init$metadata()$stan_variables)) {
Expand All @@ -1265,13 +1285,10 @@ process_init.CmdStanMCMC <- function(init, num_procs, model_variables = NULL,
warn_partial = getOption("cmdstanr_warn_inits", TRUE),
...) {
validate_fit_init(init, model_variables)
draws_df = init$draws(format = "df")
if (is.null(model_variables)) {
model_variables = list(parameters = colnames(draws_df)[2:(length(colnames(draws_df)) - 3)])
}
init_draws_df = posterior::resample_draws(draws_df, ndraws = num_procs,
draws_df <- init$draws(format = "df")
init_draws_df <- posterior::resample_draws(draws_df, ndraws = num_procs,
method = "simple_no_replace")
init_draws_lst = process_init(init_draws_df,
init_draws_lst <- process_init(init_draws_df,
num_procs = num_procs, model_variables = model_variables)
return(init_draws_lst)
}
Expand All @@ -1292,47 +1309,45 @@ process_init_approx <- function(init, num_procs, model_variables = NULL,
...) {
validate_fit_init(init, model_variables)
# Convert from data.table to data.frame
draws_df = init$draws(format = "df")
if (is.null(model_variables)) {
model_variables = list(parameters = colnames(draws_df)[3:(length(colnames(draws_df)) - 3)])
}
draws_df$lw = draws_df$lp__ - draws_df$lp_approx__
draws_df <- init$draws(format = "df")
draws_df$lw <- draws_df$lp__ - draws_df$lp_approx__
# Replace NaN and Inf with -Inf
draws_df$lw[!is.finite(draws_df$lw)] <- -Inf
# Calculate unique draws based on 'lw' using base R functions
unique_draws = length(unique(draws_df$lw))
unique_draws <- length(unique(draws_df$lw))
if (num_procs > unique_draws) {
if (inherits(init, "CmdStanPathfinder")) {
algo_name = " Pathfinder "
extra_msg = " Try running Pathfinder with psis_resample=FALSE."
algo_name <- " Pathfinder "
extra_msg <- " Try running Pathfinder with psis_resample=FALSE."
} else if (inherits(init, "CmdStanVB")) {
algo_name = " CmdStanVB "
extra_msg = ""
algo_name <- " CmdStanVB "
extra_msg <- ""
} else if (inherits(init, "CmdStanLaplace")) {
algo_name = " CmdStanLaplace "
extra_msg = ""
algo_name <- " CmdStanLaplace "
extra_msg <- ""
} else {
algo_name = ""
extra_msg = ""
algo_name <- ""
extra_msg <- ""
}
stop(paste0("Not enough distinct draws (", num_procs, ") in", algo_name ,
"fit to create inits.", extra_msg))
}
if (unique_draws < (0.95 * nrow(draws_df))) {
temp_df = stats::aggregate(.draw ~ lw, data = draws_df, FUN = min)
draws_df = posterior::as_draws_df(merge(temp_df, draws_df, by = 'lw'))
draws_df$weight = exp(draws_df$lw - max(draws_df$lw))
temp_df <- stats::aggregate(.draw ~ lw, data = draws_df, FUN = min)
draws_df <- posterior::as_draws_df(merge(temp_df, draws_df, by = 'lw'))
draws_df$weight <- exp(draws_df$lw - max(draws_df$lw))
} else {
draws_df$weight = posterior::pareto_smooth(
draws_df$weight <- posterior::pareto_smooth(
exp(draws_df$lw - max(draws_df$lw)), tail = "right", r_eff=1, return_k=FALSE)
}
init_draws_df = posterior::resample_draws(draws_df, ndraws = num_procs,
weights = draws_df$weight, method = "simple_no_replace")
init_draws_lst = process_init(init_draws_df,
num_procs = num_procs, model_variables = model_variables, warn_partial)
return(init_draws_lst)
}

init_draws_df <- posterior::resample_draws(draws_df, ndraws = num_procs,
weights = draws_df$weight, method = "simple_no_replace")
init_draws_df <- posterior::subset_draws(init_draws_df,
variable = setdiff(posterior::variables(init_draws_df), c("lw", "weight")))
init_draws_lst <- process_init(init_draws_df,
num_procs = num_procs, model_variables = model_variables, warn_partial)
return(init_draws_lst)
}

#' Write initial values to files if provided as a `CmdStanPathfinder` class
#' @noRd
Expand All @@ -1351,14 +1366,13 @@ process_init.CmdStanPathfinder <- function(init, num_procs, model_variables = NU
if (!init$metadata()$calculate_lp) {
validate_fit_init(init, model_variables)
# Convert from data.table to data.frame
draws_df = init$draws(format = "df")
if (is.null(model_variables)) {
model_variables = list(parameters = colnames(draws_df)[3:(length(colnames(draws_df)) - 3)])
}
draws_df$weight = rep(1.0, nrow(draws_df))
init_draws_df = posterior::resample_draws(draws_df, ndraws = num_procs,
draws_df <- init$draws(format = "df")
draws_df$weight <- rep(1.0, nrow(draws_df))
init_draws_df <- posterior::resample_draws(draws_df, ndraws = num_procs,
weights = draws_df$weight, method = "simple_no_replace")
init_draws_lst = process_init(init_draws_df,
init_draws_df <- posterior::subset_draws(init_draws_df,
variable = setdiff(posterior::variables(init_draws_df), "weight"))
init_draws_lst <- process_init(init_draws_df,
num_procs = num_procs, model_variables = model_variables, warn_partial)
return(init_draws_lst)
} else {
Expand Down Expand Up @@ -1417,16 +1431,31 @@ process_init.CmdStanMLE <- function(init, num_procs, model_variables = NULL,
...) {
# Convert from data.table to data.frame
validate_fit_init(init, model_variables)
draws_df = init$draws(format = "df")
if (is.null(model_variables)) {
model_variables = list(parameters = colnames(draws_df)[2:(length(colnames(draws_df)) - 3)])
}
init_draws_df = draws_df[rep(1, num_procs),]
init_draws_lst_lst = process_init(init_draws_df,
draws_df <- init$draws(format = "df")
init_draws_df <- draws_df[rep(1, num_procs),]
init_draws_lst_lst <- process_init(init_draws_df,
num_procs = num_procs, model_variables = model_variables, warn_partial)
return(init_draws_lst_lst)
}


# Validation helpers ------------------------------------------------------

#' Validate exe file exists
#' @noRd
#' @param exe_file Path to executable.
#' @return Either throws an error or returns `invisible(TRUE)`
validate_exe_file <- function(exe_file) {
if (!length(exe_file) ||
!nzchar(exe_file) ||
!file.exists(exe_file)) {
stop("Model not compiled. Try running the compile() method first.",
call. = FALSE)
}
invisible(TRUE)
}


#' Validate initial values
#'
#' For CmdStan `init` must be `NULL`, a single real number >= 0, or paths to
Expand Down
62 changes: 47 additions & 15 deletions tests/testthat/test-model-init.R
Original file line number Diff line number Diff line change
Expand Up @@ -310,33 +310,65 @@ test_that("Initial values for single-element containers treated correctly", {
)
})

test_that("Pathfinder inits do not drop dimensions", {
test_that("Inits from fit/draws work for exe-only models with various parameter types", {
modcode <- "
data {
int N;
vector[N] y;
}

parameters {
matrix[N, 1] mu;
matrix[1, N] mu_2;
vector<lower=0>[N] sigma;
real mu_scalar;
matrix[1, 1] mu_mat;
array[1] real mu_arr1;
array[1, 1] real mu_arr2;
array[1, 1, 1] real mu_arr3;
matrix[3, 1] mu_matN1;
matrix[1, 3] mu_mat1N;
}

model {
target += normal_lupdf(y | mu[:, 1], sigma);
target += normal_lupdf(y | mu_2[1], sigma);
target += normal_lupdf(mu_scalar | 0, 1);
target += normal_lupdf(to_vector(mu_mat) | 0, 1);
target += normal_lupdf(mu_arr1 | 0, 1);
target += normal_lupdf(to_array_1d(mu_arr2) | 0, 1);
target += normal_lupdf(to_array_1d(mu_arr3) | 0, 1);
target += normal_lupdf(to_vector(mu_matN1) | 0, 1);
target += normal_lupdf(to_vector(mu_mat1N) | 0, 1);
}
"
# model with stan file
mod <- cmdstan_model(write_stan_file(modcode), force_recompile = TRUE)
data <- list(N = 100, y = rnorm(100))
# Pathfinder
utils::capture.output(
pf <- mod$pathfinder(data = data, psis_resample = FALSE)
pf <- mod$pathfinder(psis_resample = FALSE, calculate_lp = FALSE)
)

# Pathfinder inits with stan file (1 chain)
expect_no_error(
utils::capture.output(
fit <- mod$sample(init = pf, chains = 1,
iter_warmup = 100, iter_sampling = 100)
)
)
# Pathfinder inits with stan file (2 chains)
expect_no_error(
utils::capture.output(
fit <- mod$sample(data = data, init = pf, chains = 1,
fit <- mod$sample(init = pf, chains = 2,
iter_warmup = 100, iter_sampling = 100)
)
)

# exe-only model (no stan file)
tmp_exe <- tempfile(fileext = cmdstan_ext())
file.copy(mod$exe_file(), tmp_exe)
mod_nostan <- cmdstan_model(exe_file = tmp_exe)
# Pathfinder inits without stan file (1 chain)
expect_no_error(
utils::capture.output(
fit <- mod_nostan$sample(init = pf, chains = 1,
iter_warmup = 100, iter_sampling = 100)
)
)
# Pathfinder inits without stan file (2 chains)
expect_no_error(
utils::capture.output(
fit <- mod_nostan$sample(init = pf, chains = 2,
iter_warmup = 100, iter_sampling = 100)
)
)
})
Loading