From d1762743f72d324c2cd5d1c571b943117493ef86 Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Thu, 23 Oct 2025 10:14:35 -0400 Subject: [PATCH 1/3] CI Updates --- .github/workflows/main.yml | 33 +++++++++++++++++++++++++----- Makefile | 13 ++++++++++-- clients/typescript/.prettierignore | 2 +- 3 files changed, 40 insertions(+), 8 deletions(-) diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 09c7fab..25626fd 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -8,7 +8,7 @@ on: workflow_dispatch: {} env: - CACHE_VERSION: 3 + CACHE_VERSION: 4 # only run one copy per PR concurrency: @@ -16,6 +16,29 @@ concurrency: cancel-in-progress: true jobs: + check-formatting: + runs-on: ubuntu-latest + steps: + - name: Check out github + uses: actions/checkout@v5 + with: + submodules: recursive + + - name: Set up Julia + uses: julia-actions/install-juliaup@v2 + with: + channel: "lts" + + - name: Install Formatters + run: | + sudo apt install -y clang-format black isort + curl --proto '=https' --tlsv1.2 -LsSf https://github.com/posit-dev/air/releases/download/0.8.0/air-installer.sh | sh + julia --project=clients/julia -e 'import Pkg; Pkg.add("JuliaFormatter")' + + - name: Run formatting checks + run: | + make format-check || (echo "::error title=Formatting::Formatting check failed, run \`make format\` locally" && exit 1) + build: runs-on: ${{matrix.os}} strategy: @@ -56,7 +79,7 @@ jobs: strategy: matrix: os: [ubuntu-latest, macos-latest, windows-latest] - python-version: ["3.9", "3.12"] + python-version: ["3.9", "3"] fail-fast: false steps: - name: Check out github @@ -175,9 +198,9 @@ jobs: submodules: recursive - name: Set up Julia - uses: julia-actions/setup-julia@v2 + uses: julia-actions/install-juliaup@v2 with: - version: ${{ matrix.julia-version }} + channel: ${{ matrix.julia-version }} - name: Restore Stan uses: actions/cache@v4 @@ -201,7 +224,7 @@ jobs: - name: Run tests run: | - julia --project=clients/julia -t 2 -e "using Pkg; Pkg.test()" + julia --project=clients/julia -t 2 --color yes -e "using Pkg; Pkg.test()" env: TINYSTAN: ${{ github.workspace }} diff --git a/Makefile b/Makefile index 787d82a..0f9a344 100644 --- a/Makefile +++ b/Makefile @@ -124,13 +124,22 @@ TEST_MODEL_LIBS = $(join $(addprefix test_models/, $(TEST_MODEL_NAMES)), $(addsu test_models: $(TEST_MODEL_LIBS) -.PHONY: format +.PHONY: format format-check format: clang-format -i src/*.cpp src/*.hpp src/*.h || true isort clients/python || true black clients/python || true julia --project=clients/julia -e 'using JuliaFormatter; format("clients/julia/")' || true - Rscript -e 'formatR::tidy_dir("clients/R/", recursive=TRUE)' || true + air format clients/R || true + cd clients/typescript && npx prettier --write . || true + +format-check: + clang-format --dry-run --Werror src/*.cpp src/*.hpp src/*.h + isort clients/python --check + black clients/python --check + julia --project=clients/julia -e 'using JuliaFormatter; if !format("clients/julia/", overwrite=false) exit(1) end' + air format clients/R --check + cd clients/typescript && npx prettier --check . .PHONY: clean clean: diff --git a/clients/typescript/.prettierignore b/clients/typescript/.prettierignore index a7e5fab..c301472 100644 --- a/clients/typescript/.prettierignore +++ b/clients/typescript/.prettierignore @@ -1,3 +1,3 @@ dist/ coverage/ - +doc/*.hbs From 6d95892448682fbfe028b384377ff81b1424abb8 Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Thu, 23 Oct 2025 12:13:58 -0400 Subject: [PATCH 2/3] workaround: use windows-2022 for Julia tests --- .github/workflows/main.yml | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 25626fd..bb22ac6 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -8,7 +8,7 @@ on: workflow_dispatch: {} env: - CACHE_VERSION: 4 + CACHE_VERSION: 5 # only run one copy per PR concurrency: @@ -44,7 +44,9 @@ jobs: strategy: fail-fast: false matrix: - os: [windows-latest, ubuntu-latest, macos-latest] + # windows-2022 is a workaround for a julia issue as of 10/23/2025 + # https://github.com/JuliaLang/julia/issues/59931 + os: [windows-2022, ubuntu-latest, macos-latest] env: GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} steps: @@ -78,7 +80,7 @@ jobs: runs-on: ${{matrix.os}} strategy: matrix: - os: [ubuntu-latest, macos-latest, windows-latest] + os: [ubuntu-latest, macos-latest, windows-2022] python-version: ["3.9", "3"] fail-fast: false steps: @@ -133,7 +135,7 @@ jobs: runs-on: ${{matrix.os}} strategy: matrix: - os: [ubuntu-latest, macos-latest, windows-latest] + os: [ubuntu-latest, macos-latest, windows-2022] fail-fast: false steps: - name: Check out github @@ -167,7 +169,7 @@ jobs: key: ${{ hashFiles('**/*.stan', 'src/*', 'stan/src/stan/version.hpp', 'Makefile') }}-${{ matrix.os }}-v${{ env.CACHE_VERSION }} - name: Run tests - if: matrix.os != 'windows-latest' + if: matrix.os != 'windows-2022' run: | cd clients/R Rscript -e "devtools::test(reporter = c(\"summary\", \"fail\"))" @@ -175,7 +177,7 @@ jobs: TINYSTAN: ${{ github.workspace }} - name: Run tests (windows) - if: matrix.os == 'windows-latest' + if: matrix.os == 'windows-2022' run: | cd clients/R Rscript -e 'devtools::test(reporter = c("summary", "fail"))' @@ -188,7 +190,7 @@ jobs: runs-on: ${{matrix.os}} strategy: matrix: - os: [ubuntu-latest, macos-latest, windows-latest] + os: [ubuntu-latest, macos-latest, windows-2022] julia-version: ["1"] fail-fast: false steps: From 90ae37640c5b92ee1653d4af75e1d5c9843554d0 Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Thu, 23 Oct 2025 10:14:41 -0400 Subject: [PATCH 3/3] Format --- clients/R/R/compile.R | 156 ++-- clients/R/R/download.R | 82 +- clients/R/R/output.R | 25 +- clients/R/R/tinystan.R | 867 +++++++++++++-------- clients/R/R/zzz.R | 11 +- clients/R/convert_docs.R | 17 +- clients/R/tests/testthat.R | 1 - clients/R/tests/testthat/setup.R | 32 +- clients/R/tests/testthat/test_compile.R | 39 +- clients/R/tests/testthat/test_laplace.R | 201 +++-- clients/R/tests/testthat/test_model.R | 15 +- clients/R/tests/testthat/test_optimize.R | 266 ++++--- clients/R/tests/testthat/test_pathfinder.R | 324 +++++--- clients/R/tests/testthat/test_sample.R | 607 ++++++++++----- clients/julia/src/model.jl | 5 +- clients/julia/test/test_laplace.jl | 6 +- clients/julia/test/test_pathfinder.jl | 4 +- clients/julia/test/test_sample.jl | 8 +- clients/python/tinystan/model.py | 4 +- src/buffer.hpp | 10 +- src/errors.hpp | 4 +- 21 files changed, 1713 insertions(+), 971 deletions(-) diff --git a/clients/R/R/compile.R b/clients/R/R/compile.R index 4743c4d..3a4461c 100644 --- a/clients/R/R/compile.R +++ b/clients/R/R/compile.R @@ -3,17 +3,27 @@ MAKE <- Sys.getenv("MAKE", "make") verify_tinystan_path <- function(path) { - suppressWarnings({ - folder <- normalizePath(path) - }) - if (!dir.exists(folder)) { - stop(paste0("TinyStan folder '", folder, "' does not exist!\n", "If you need to set a different location, call 'set_tinystan_path()'")) - } - makefile <- file.path(folder, "Makefile") - if (!file.exists(makefile)) { - stop(paste0("TinyStan folder '", folder, "' does not contain file 'Makefile',", - " please ensure it is built properly!\n", "If you need to set a different location, call 'set_tinystan_path()'")) - } + suppressWarnings({ + folder <- normalizePath(path) + }) + if (!dir.exists(folder)) { + stop(paste0( + "TinyStan folder '", + folder, + "' does not exist!\n", + "If you need to set a different location, call 'set_tinystan_path()'" + )) + } + makefile <- file.path(folder, "Makefile") + if (!file.exists(makefile)) { + stop(paste0( + "TinyStan folder '", + folder, + "' does not contain file 'Makefile',", + " please ensure it is built properly!\n", + "If you need to set a different location, call 'set_tinystan_path()'" + )) + } } #' @title Function `set_tinystan_path()` @@ -21,8 +31,8 @@ verify_tinystan_path <- function(path) { #' @details This should point to the top-level folder of the repository. #' @export set_tinystan_path <- function(path) { - verify_tinystan_path(path) - Sys.setenv(TINYSTAN = normalizePath(path)) + verify_tinystan_path(path) + Sys.setenv(TINYSTAN = normalizePath(path)) } #' Get the path to TinyStan. @@ -36,21 +46,28 @@ set_tinystan_path <- function(path) { #' #' @seealso [set_tinystan_path] get_tinystan_path <- function() { - # try to get from environment - path <- Sys.getenv("TINYSTAN", unset = "") - if (path == "") { - path <- CURRENT_TINYSTAN - tryCatch({ - verify_tinystan_path(path) - }, error = function(e) { - print(paste0("TinyStan not found at location specified by $TINYSTAN ", - "environment variable, downloading version ", packageVersion("tinystan"), - " to ", path)) - get_tinystan_src() - }) - } + # try to get from environment + path <- Sys.getenv("TINYSTAN", unset = "") + if (path == "") { + path <- CURRENT_TINYSTAN + tryCatch( + { + verify_tinystan_path(path) + }, + error = function(e) { + print(paste0( + "TinyStan not found at location specified by $TINYSTAN ", + "environment variable, downloading version ", + packageVersion("tinystan"), + " to ", + path + )) + get_tinystan_src() + } + ) + } - return(path) + return(path) } @@ -73,52 +90,67 @@ get_tinystan_path <- function() { #' @seealso [tinystan::set_tinystan_path()] #' @export compile_model <- function(stan_file, stanc_args = NULL, make_args = NULL) { - verify_tinystan_path(get_tinystan_path()) - suppressWarnings({ - file_path <- normalizePath(stan_file) - }) - if (tools::file_ext(file_path) != "stan") { - stop(paste0("File '", file_path, "' does not end with '.stan'")) - } - if (!file.exists(file_path)) { - stop(paste0("File '", file_path, "' does not exist!")) - } + verify_tinystan_path(get_tinystan_path()) + suppressWarnings({ + file_path <- normalizePath(stan_file) + }) + if (tools::file_ext(file_path) != "stan") { + stop(paste0("File '", file_path, "' does not end with '.stan'")) + } + if (!file.exists(file_path)) { + stop(paste0("File '", file_path, "' does not exist!")) + } - output <- paste0(tools::file_path_sans_ext(file_path), "_model.so") - stancflags <- paste("--include-paths=.", paste(stanc_args, collapse = " ")) + output <- paste0(tools::file_path_sans_ext(file_path), "_model.so") + stancflags <- paste("--include-paths=.", paste(stanc_args, collapse = " ")) - flags <- c(paste("-C", get_tinystan_path()), make_args, paste0("STANCFLAGS=\"", - stancflags, "\""), output) + flags <- c( + paste("-C", get_tinystan_path()), + make_args, + paste0("STANCFLAGS=\"", stancflags, "\""), + output + ) - suppressWarnings({ - res <- system2(MAKE, args = flags, stdout = TRUE, stderr = TRUE) - }) - res_attrs <- attributes(res) - if ("status" %in% names(res_attrs) && res_attrs$status != 0) { - stop(paste0("Compilation failed with error code ", res_attrs$status, "\noutput:\n", - paste(res, collapse = "\n"))) - } + suppressWarnings({ + res <- system2(MAKE, args = flags, stdout = TRUE, stderr = TRUE) + }) + res_attrs <- attributes(res) + if ("status" %in% names(res_attrs) && res_attrs$status != 0) { + stop(paste0( + "Compilation failed with error code ", + res_attrs$status, + "\noutput:\n", + paste(res, collapse = "\n") + )) + } - return(output) + return(output) } tbb_found <- function() { - suppressWarnings(out <- system2("where.exe", "tbb.dll", stdout = NULL, stderr = NULL)) - return(out == 0) + suppressWarnings( + out <- system2("where.exe", "tbb.dll", stdout = NULL, stderr = NULL) + ) + return(out == 0) } WINDOWS_PATH_SET <- FALSE windows_dll_path_setup <- function() { - if (.Platform$OS.type == "windows" && !WINDOWS_PATH_SET) { - - if (tbb_found()) { - assign("WINDOWS_PATH_SET", TRUE, envir = .GlobalEnv) - } else { - tbb_path <- file.path(get_tinystan_path(), "stan", "lib", "stan_math", - "lib", "tbb") - Sys.setenv(PATH = paste(tbb_path, Sys.getenv("PATH"), sep = ";")) - assign("WINDOWS_PATH_SET", tbb_found(), envir = .GlobalEnv) - } + if (.Platform$OS.type == "windows" && !WINDOWS_PATH_SET) { + if (tbb_found()) { + assign("WINDOWS_PATH_SET", TRUE, envir = .GlobalEnv) + } else { + tbb_path <- file.path( + get_tinystan_path(), + "stan", + "lib", + "stan_math", + "lib", + "tbb" + ) + Sys.setenv(PATH = paste(tbb_path, Sys.getenv("PATH"), sep = ";")) + assign("WINDOWS_PATH_SET", tbb_found(), envir = .GlobalEnv) } + } } diff --git a/clients/R/R/download.R b/clients/R/R/download.R index edaf56d..b5448b6 100644 --- a/clients/R/R/download.R +++ b/clients/R/R/download.R @@ -1,37 +1,65 @@ current_version <- packageVersion("tinystan") -current_version_list <- list(major = current_version$major, minor = current_version$minor, - patch = current_version$patch) +current_version_list <- list( + major = current_version$major, + minor = current_version$minor, + patch = current_version$patch +) HOME_TINYSTAN <- path.expand(file.path("~", ".tinystan")) -CURRENT_TINYSTAN <- file.path(HOME_TINYSTAN, paste0("tinystan-", current_version)) +CURRENT_TINYSTAN <- file.path( + HOME_TINYSTAN, + paste0("tinystan-", current_version) +) RETRIES <- 5 get_tinystan_src <- function() { - url <- paste0("https://github.com/WardBrian/tinystan/releases/download/", "v", - current_version, "/tinystan-", current_version, ".tar.gz") + url <- paste0( + "https://github.com/WardBrian/tinystan/releases/download/", + "v", + current_version, + "/tinystan-", + current_version, + ".tar.gz" + ) - dir.create(HOME_TINYSTAN, showWarnings = FALSE, recursive = TRUE) - temp <- tempfile() - err_text <- paste("Failed to download TinyStan", current_version, "from github.com.") - for (i in 1:RETRIES) { - tryCatch({ - download.file(url, destfile = temp, mode = "wb", quiet = TRUE, method = "auto") - }, error = function(e) { - cat(err_text, "\n") - if (i == RETRIES) { - stop(err_text, call. = FALSE) - } else { - cat("Retrying (", i + 1, "/", RETRIES, ")...\n", sep = "") - Sys.sleep(1) - } - }) - } + dir.create(HOME_TINYSTAN, showWarnings = FALSE, recursive = TRUE) + temp <- tempfile() + err_text <- paste( + "Failed to download TinyStan", + current_version, + "from github.com." + ) + for (i in 1:RETRIES) { + tryCatch( + { + download.file( + url, + destfile = temp, + mode = "wb", + quiet = TRUE, + method = "auto" + ) + }, + error = function(e) { + cat(err_text, "\n") + if (i == RETRIES) { + stop(err_text, call. = FALSE) + } else { + cat("Retrying (", i + 1, "/", RETRIES, ")...\n", sep = "") + Sys.sleep(1) + } + } + ) + } - tryCatch({ - untar(temp, exdir = HOME_TINYSTAN) - }, error = function(e) { - stop(paste("Failed to unpack", url, "during installation"), call. = FALSE) - }) + tryCatch( + { + untar(temp, exdir = HOME_TINYSTAN) + }, + error = function(e) { + stop(paste("Failed to unpack", url, "during installation"), call. = FALSE) + } + ) - unlink(temp) + unlink(temp) } diff --git a/clients/R/R/output.R b/clients/R/R/output.R index 25f7fd0..55c81cd 100644 --- a/clients/R/R/output.R +++ b/clients/R/R/output.R @@ -1,23 +1,22 @@ - # copied from cmdstanr, definitely doesn't handle tuples, but then neither does # posterior repair_variable_names <- function(names) { - names <- sub("\\.", "[", names) - names <- gsub("\\.", ",", names) - names[grep("\\[", names)] <- paste0(names[grep("\\[", names)], "]") - names + names <- sub("\\.", "[", names) + names <- gsub("\\.", ",", names) + names[grep("\\[", names)] <- paste0(names[grep("\\[", names)], "]") + names } output_as_rvars <- function(names, num_draws, num_chains, draws) { - names <- repair_variable_names(names) - num_params <- length(names) - dims <- c(num_params, num_draws, num_chains) + names <- repair_variable_names(names) + num_params <- length(names) + dims <- c(num_params, num_draws, num_chains) - # all our outputs are row-major - draws <- array(draws, dim = dims, dimnames = list(names, NULL, NULL)) - # so we need to rearrange. posterior likes draws x chains x params - draws <- aperm(draws, c(2, 3, 1)) + # all our outputs are row-major + draws <- array(draws, dim = dims, dimnames = list(names, NULL, NULL)) + # so we need to rearrange. posterior likes draws x chains x params + draws <- aperm(draws, c(2, 3, 1)) - posterior::as_draws_rvars(draws) + posterior::as_draws_rvars(draws) } diff --git a/clients/R/R/tinystan.R b/clients/R/R/tinystan.R index f651d65..3dc0124 100644 --- a/clients/R/R/tinystan.R +++ b/clients/R/R/tinystan.R @@ -7,133 +7,197 @@ #' fit = sampler(model = mod, data = data_file) #' fit #' -tinystan_model = function(lib, stanc_args = NULL, make_args = NULL, warn = TRUE) { - if (tools::file_ext(lib) == "stan") { - # FIXME: This is a hack to get the code from the stan file - stan_code <- lib - lib <- compile_model(lib, stanc_args, make_args) - code <- paste0(readLines(stan_code), collapse = "\n") - built_with_so <- FALSE - } else { - build_with_so <- TRUE - code <- "Built with .so object. No code available for printing" - } - if (.Platform$OS.type == "windows") { - lib_old <- lib - lib <- paste0(tools::file_path_sans_ext(lib), ".dll") - file.copy(from = lib_old, to = lib) - windows_dll_path_setup() - } - lib <- tools::file_path_as_absolute(lib) - lib_name <- tools::file_path_sans_ext(basename(lib)) - if (warn && is.loaded("tinystan_create_model_R", PACKAGE = lib_name)) { - warning(paste0("Loading a shared object '", lib, "' which is already loaded.\n", - "If the file has changed since the last time it was loaded, this load may not update the library!")) - } - dyn.load(lib, PACKAGE = lib_name) - sep <- .C("tinystan_separator_char_R", sep = raw(1), PACKAGE = lib_name)$sep - sep <- rawToChar(sep) - - api_ver = .C("tinystan_api_version", major = integer(1), minor = integer(1), - patch = integer(1), PACKAGE = lib_name) - - if (api_ver$major != current_version_list$major) { - msg = paste0("Incompatible TinyStan API version. Expected ", paste(current_version_list, - collapse = "."), " but found ", paste(api_ver, collapse = "."), ".\nYou need to re-compile your model.") - stop(msg) - } else if (api_ver$minor != current_version_list$minor || api_ver$patch != current_version_list$patch) { - msg = paste0("TinyStan API version mismatch. Expected ", paste(current_version_list, - collapse = "."), " but found ", paste(api_ver, collapse = "."), ".\nYou may need to re-compile your model.") - warning(msg) - } - - ret <- list(lib = lib, lib_name = lib_name, sep = sep, code = code, built_with_so = built_with_so) - class(ret) <- c("tinystan_model", class(ret)) - return(ret) +tinystan_model = function( + lib, + stanc_args = NULL, + make_args = NULL, + warn = TRUE +) { + if (tools::file_ext(lib) == "stan") { + # FIXME: This is a hack to get the code from the stan file + stan_code <- lib + lib <- compile_model(lib, stanc_args, make_args) + code <- paste0(readLines(stan_code), collapse = "\n") + built_with_so <- FALSE + } else { + build_with_so <- TRUE + code <- "Built with .so object. No code available for printing" + } + if (.Platform$OS.type == "windows") { + lib_old <- lib + lib <- paste0(tools::file_path_sans_ext(lib), ".dll") + file.copy(from = lib_old, to = lib) + windows_dll_path_setup() + } + lib <- tools::file_path_as_absolute(lib) + lib_name <- tools::file_path_sans_ext(basename(lib)) + if (warn && is.loaded("tinystan_create_model_R", PACKAGE = lib_name)) { + warning(paste0( + "Loading a shared object '", + lib, + "' which is already loaded.\n", + "If the file has changed since the last time it was loaded, this load may not update the library!" + )) + } + dyn.load(lib, PACKAGE = lib_name) + sep <- .C("tinystan_separator_char_R", sep = raw(1), PACKAGE = lib_name)$sep + sep <- rawToChar(sep) + + api_ver = .C( + "tinystan_api_version", + major = integer(1), + minor = integer(1), + patch = integer(1), + PACKAGE = lib_name + ) + + if (api_ver$major != current_version_list$major) { + msg = paste0( + "Incompatible TinyStan API version. Expected ", + paste(current_version_list, collapse = "."), + " but found ", + paste(api_ver, collapse = "."), + ".\nYou need to re-compile your model." + ) + stop(msg) + } else if ( + api_ver$minor != current_version_list$minor || + api_ver$patch != current_version_list$patch + ) { + msg = paste0( + "TinyStan API version mismatch. Expected ", + paste(current_version_list, collapse = "."), + " but found ", + paste(api_ver, collapse = "."), + ".\nYou may need to re-compile your model." + ) + warning(msg) + } + + ret <- list( + lib = lib, + lib_name = lib_name, + sep = sep, + code = code, + built_with_so = built_with_so + ) + class(ret) <- c("tinystan_model", class(ret)) + return(ret) } #' @export print.tinystan_model <- function(model, ...) { - cat(model$code, ...) - if (model$built_with_so) { - cat("Library: ", model$lib, "\n") - } + cat(model$code, ...) + if (model$built_with_so) { + cat("Library: ", model$lib, "\n") + } } #'@export api_version = function(stan_model) { - .C("tinystan_api_version", major = integer(1), minor = integer(1), patch = integer(1), - PACKAGE = stan_model$lib_name) + .C( + "tinystan_api_version", + major = integer(1), + minor = integer(1), + patch = integer(1), + PACKAGE = stan_model$lib_name + ) } #'@export stan_version = function(stan_model) { - .C("tinystan_stan_version", major = integer(1), minor = integer(1), patch = integer(1), - PACKAGE = stan_model$lib_name) + .C( + "tinystan_stan_version", + major = integer(1), + minor = integer(1), + patch = integer(1), + PACKAGE = stan_model$lib_name + ) } #' @noRd with_model = function(model, data, seed, block) { - ffi_ret <- .C("tinystan_create_model_R", model_ptr = raw(8), as.character(data), - as.integer(seed), err = raw(8), NAOK = TRUE, PACKAGE = model$lib_name) - handle_error(all(ffi_ret$model_ptr == 0), model$lib_name, ffi_ret$err) - tryCatch({ - # this is the equivalent of base R's `with` function - eval(substitute(block), ffi_ret, enclos = parent.frame()) - }, finally = { - .C("tinystan_destroy_model_R", as.raw(ffi_ret$model_ptr), PACKAGE = model$lib_name) - }) + ffi_ret <- .C( + "tinystan_create_model_R", + model_ptr = raw(8), + as.character(data), + as.integer(seed), + err = raw(8), + NAOK = TRUE, + PACKAGE = model$lib_name + ) + handle_error(all(ffi_ret$model_ptr == 0), model$lib_name, ffi_ret$err) + tryCatch( + { + # this is the equivalent of base R's `with` function + eval(substitute(block), ffi_ret, enclos = parent.frame()) + }, + finally = { + .C( + "tinystan_destroy_model_R", + as.raw(ffi_ret$model_ptr), + PACKAGE = model$lib_name + ) + } + ) } #' @noRd get_parameter_names <- function(...) { - UseMethod("get_parameter_names") + UseMethod("get_parameter_names") } #' @noRd get_parameter_names.tinystan_model = function(model, model_ptr) { - param_names_raw <- .C("tinystan_model_param_names_R", as.raw(model_ptr), names = as.character(""), - PACKAGE = model$lib_name)$names - if (param_names_raw == "") { - return(c()) - } - strsplit(param_names_raw, ",")[[1]] - + param_names_raw <- .C( + "tinystan_model_param_names_R", + as.raw(model_ptr), + names = as.character(""), + PACKAGE = model$lib_name + )$names + if (param_names_raw == "") { + return(c()) + } + strsplit(param_names_raw, ",")[[1]] } #' @noRd get_free_params <- function(...) { - UseMethod("get_free_params") + UseMethod("get_free_params") } #' @noRd get_free_params.tinystan_model = function(model, model_ptr) { - .C("tinystan_model_num_free_params_R", as.raw(model_ptr), params = as.integer(0), - PACKAGE = model$lib_name)$params + .C( + "tinystan_model_num_free_params_R", + as.raw(model_ptr), + params = as.integer(0), + PACKAGE = model$lib_name + )$params } #' @noRd encode_inits = function(model, inits) { - if (is.null(inits)) { - return(as.character("")) - } - if (is.character(inits)) { - if (length(inits) == 1) { - return(as.character(inits)) - } else { - return(as.character(paste0(inits, collapse = model$sep))) - } - } - if (is.list(inits)) { - return(as.character(paste0(inits, collapse = model$sep))) + if (is.null(inits)) { + return(as.character("")) + } + if (is.character(inits)) { + if (length(inits) == 1) { + return(as.character(inits)) + } else { + return(as.character(paste0(inits, collapse = model$sep))) } - stop("inits must be a character vector or a list") + } + if (is.list(inits)) { + return(as.character(paste0(inits, collapse = model$sep))) + } + stop("inits must be a character vector or a list") } #' @export sampler <- function(...) { - UseMethod("sampler") + UseMethod("sampler") } #' @title Generic function `sampler` @@ -144,278 +208,459 @@ sampler <- function(...) { #' mod <- tinystan_model(system.file('bernoulli.stan', package = 'tinystan')) #' fit = sampler(model = mod, data = data_file) #' fit -sampler.tinystan_model = function(model, data = "", num_chains = 4, inits = NULL, - seed = NULL, id = 1, init_radius = 2, num_warmup = 1000, num_samples = 1000, - metric = HMCMetric$DIAGONAL, init_inv_metric = NULL, save_inv_metric = FALSE, - adapt = TRUE, delta = 0.8, gamma = 0.05, kappa = 0.75, t0 = 10, init_buffer = 75, - term_buffer = 50, window = 25, save_warmup = FALSE, stepsize = 1, stepsize_jitter = 0, - max_depth = 10, refresh = 0, num_threads = -1) { - - if (num_chains < 1) { - stop("num_chains must be at least 1") - } - if (save_warmup && num_warmup < 0) { - stop("num_warmup must be non-negative") - } - if (num_samples < 1) { - stop("num_samples must be at least 1") +sampler.tinystan_model = function( + model, + data = "", + num_chains = 4, + inits = NULL, + seed = NULL, + id = 1, + init_radius = 2, + num_warmup = 1000, + num_samples = 1000, + metric = HMCMetric$DIAGONAL, + init_inv_metric = NULL, + save_inv_metric = FALSE, + adapt = TRUE, + delta = 0.8, + gamma = 0.05, + kappa = 0.75, + t0 = 10, + init_buffer = 75, + term_buffer = 50, + window = 25, + save_warmup = FALSE, + stepsize = 1, + stepsize_jitter = 0, + max_depth = 10, + refresh = 0, + num_threads = -1 +) { + if (num_chains < 1) { + stop("num_chains must be at least 1") + } + if (save_warmup && num_warmup < 0) { + stop("num_warmup must be non-negative") + } + if (num_samples < 1) { + stop("num_samples must be at least 1") + } + + if (is.null(seed)) { + seed <- as.integer(runif(1, min = 0, max = (2^31))) + } + + with_model(model, data, seed, { + free_params <- get_free_params(model, model_ptr) + + params <- c(HMC_SAMPLER_VARIABLES, get_parameter_names(model, model_ptr)) + num_params <- length(params) + num_draws <- as.integer(save_warmup) * num_warmup + num_samples + output_size <- num_params * num_chains * num_draws + + if (metric == HMCMetric$DENSE) { + metric_shape <- rep(free_params, 2) + } else { + metric_shape <- free_params } - if (is.null(seed)) { - seed <- as.integer(runif(1, min = 0, max = (2^31))) + if (is.null(init_inv_metric)) { + metric_has_init <- FALSE + inv_metric_init <- 0 + } else { + metric_has_init <- TRUE + metric_dims <- dim(init_inv_metric) + + if ( + length(metric_dims) == length(metric_shape) && + all(metric_dims == metric_shape) + ) { + inv_metric_init <- replicate(num_chains, init_inv_metric) + } else if ( + length(metric_dims) == (length(metric_shape) + 1) && + all(metric_dims == c(metric_shape, num_chains)) + ) { + inv_metric_init <- init_inv_metric + } else { + stop( + "Invalid initial metric size. Expected a ", + paste(metric_shape, collapse = " x "), + " or ", + paste(c(metric_shape, num_chains), collapse = " x "), + " matrix" + ) + } + } + metric_size <- 1 + stepsizes_size <- 1 + if (adapt) { + stepsizes_size <- num_chains + if (save_inv_metric) { + metric_size <- num_chains * prod(metric_shape) + } } - with_model(model, data, seed, { - free_params <- get_free_params(model, model_ptr) - - params <- c(HMC_SAMPLER_VARIABLES, get_parameter_names(model, model_ptr)) - num_params <- length(params) - num_draws <- as.integer(save_warmup) * num_warmup + num_samples - output_size <- num_params * num_chains * num_draws - + vars <- .C( + "tinystan_sample_R", + return_code = as.integer(0), + as.raw(model_ptr), + as.integer(num_chains), + encode_inits(model, inits), + as.integer(seed), + as.integer(id), + as.double(init_radius), + as.integer(num_warmup), + as.integer(num_samples), + as.integer(metric), + as.logical(metric_has_init), + as.double(inv_metric_init), + as.logical(adapt), + as.double(delta), + as.double(gamma), + as.double(kappa), + as.double(t0), + as.integer(init_buffer), + as.integer(term_buffer), + as.integer(window), + as.logical(save_warmup), + as.double(stepsize), + as.double(stepsize_jitter), + as.integer(max_depth), + as.integer(refresh), + as.integer(num_threads), + out = double(output_size), + as.integer(output_size), + save_stepsizes = as.logical(adapt), + stepsize_out = double(stepsizes_size), + save_inv_metric = as.logical( + adapt && + save_inv_metric + ), + inv_metric = double(metric_size), + err = raw(8), + PACKAGE = model$lib_name + ) + handle_error(vars$return_code, model$lib_name, vars$err) + # reshape the output matrix + out <- output_as_rvars(params, num_draws, num_chains, vars$out) + + if (adapt) { + if (save_inv_metric) { if (metric == HMCMetric$DENSE) { - metric_shape <- rep(free_params, 2) + inv_metric <- aperm( + array( + vars$inv_metric, + dim = c(free_params, free_params, num_chains) + ), + c(3, 2, 1) + ) } else { - metric_shape <- free_params + inv_metric <- aperm( + array(vars$inv_metric, dim = c(free_params, num_chains)), + c(2, 1) + ) } - - if (is.null(init_inv_metric)) { - metric_has_init <- FALSE - inv_metric_init <- 0 - } else { - metric_has_init <- TRUE - metric_dims <- dim(init_inv_metric) - - if (length(metric_dims) == length(metric_shape) && all(metric_dims == - metric_shape)) { - inv_metric_init <- replicate(num_chains, init_inv_metric) - } else if (length(metric_dims) == (length(metric_shape) + 1) && all(metric_dims == - c(metric_shape, num_chains))) { - inv_metric_init <- init_inv_metric - } else { - stop("Invalid initial metric size. Expected a ", paste(metric_shape, - collapse = " x "), " or ", paste(c(metric_shape, num_chains), collapse = " x "), - " matrix") - } - } - metric_size <- 1 - stepsizes_size <- 1 - if (adapt) { - stepsizes_size <- num_chains - if (save_inv_metric) { - metric_size <- num_chains * prod(metric_shape) - } - } - - vars <- .C("tinystan_sample_R", return_code = as.integer(0), as.raw(model_ptr), - as.integer(num_chains), encode_inits(model, inits), as.integer(seed), - as.integer(id), as.double(init_radius), as.integer(num_warmup), as.integer(num_samples), - as.integer(metric), as.logical(metric_has_init), as.double(inv_metric_init), - as.logical(adapt), as.double(delta), as.double(gamma), as.double(kappa), - as.double(t0), as.integer(init_buffer), as.integer(term_buffer), as.integer(window), - as.logical(save_warmup), as.double(stepsize), as.double(stepsize_jitter), - as.integer(max_depth), as.integer(refresh), as.integer(num_threads), - out = double(output_size), as.integer(output_size), save_stepsizes = as.logical(adapt), - stepsize_out = double(stepsizes_size), save_inv_metric = as.logical(adapt && - save_inv_metric), inv_metric = double(metric_size), err = raw(8), - PACKAGE = model$lib_name) - handle_error(vars$return_code, model$lib_name, vars$err) - # reshape the output matrix - out <- output_as_rvars(params, num_draws, num_chains, vars$out) - - if (adapt) { - if (save_inv_metric) { - if (metric == HMCMetric$DENSE) { - inv_metric <- aperm(array(vars$inv_metric, dim = c(free_params, - free_params, num_chains)), c(3, 2, 1)) - } else { - inv_metric <- aperm(array(vars$inv_metric, dim = c(free_params, - num_chains)), c(2, 1)) - } - return(list(draws = out, stepsize = vars$stepsize_out, inv_metric = inv_metric)) - } - return(list(draws = out, stepsize = vars$stepsize_out)) - } - list(draws = out) - }) + return(list( + draws = out, + stepsize = vars$stepsize_out, + inv_metric = inv_metric + )) + } + return(list(draws = out, stepsize = vars$stepsize_out)) + } + list(draws = out) + }) } #' @export pathfinder <- function(...) { - UseMethod("pathfinder") + UseMethod("pathfinder") } #' @title Generic function `pathfinder` #' @description Run Stan's pathfinder algorithm #' @export -pathfinder.tinystan_model = function(model, data = "", num_paths = 4, inits = NULL, - seed = NULL, id = 1, init_radius = 2, num_draws = 1000, max_history_size = 5, - init_alpha = 0.001, tol_obj = 1e-12, tol_rel_obj = 10000, tol_grad = 1e-08, tol_rel_grad = 1e+07, - tol_param = 1e-08, num_iterations = 1000, num_elbo_draws = 25, num_multi_draws = 1000, - calculate_lp = TRUE, psis_resample = TRUE, refresh = 0, num_threads = -1) { - if (num_draws < 1) { - stop("num_draws must be at least 1") +pathfinder.tinystan_model = function( + model, + data = "", + num_paths = 4, + inits = NULL, + seed = NULL, + id = 1, + init_radius = 2, + num_draws = 1000, + max_history_size = 5, + init_alpha = 0.001, + tol_obj = 1e-12, + tol_rel_obj = 10000, + tol_grad = 1e-08, + tol_rel_grad = 1e+07, + tol_param = 1e-08, + num_iterations = 1000, + num_elbo_draws = 25, + num_multi_draws = 1000, + calculate_lp = TRUE, + psis_resample = TRUE, + refresh = 0, + num_threads = -1 +) { + if (num_draws < 1) { + stop("num_draws must be at least 1") + } + if (num_paths < 1) { + stop("num_paths must be at least 1") + } + if (num_multi_draws < 1) { + stop("num_multi_draws must be at least 1") + } + + if (calculate_lp && psis_resample) { + num_output <- num_multi_draws + } else { + num_output <- num_draws * num_paths + } + if (is.null(seed)) { + seed <- as.integer(runif(1, min = 0, max = (2^31))) + } + + with_model(model, data, seed, { + free_params <- get_free_params(model, model_ptr) + if (free_params == 0) { + stop("Model has no parameters") } - if (num_paths < 1) { - stop("num_paths must be at least 1") - } - if (num_multi_draws < 1) { - stop("num_multi_draws must be at least 1") - } - - if (calculate_lp && psis_resample) { - num_output <- num_multi_draws - } else { - num_output <- num_draws * num_paths - } - if (is.null(seed)) { - seed <- as.integer(runif(1, min = 0, max = (2^31))) - } - - with_model(model, data, seed, { - free_params <- get_free_params(model, model_ptr) - if (free_params == 0) { - stop("Model has no parameters") - } - params <- c(PATHFINDER_VARIABLES, get_parameter_names(model, model_ptr)) - num_params <- length(params) - output_size <- num_params * num_output - - vars <- .C("tinystan_pathfinder_R", return_code = as.integer(0), as.raw(model_ptr), - as.integer(num_paths), encode_inits(model, inits), as.integer(seed), - as.integer(id), as.double(init_radius), as.integer(num_draws), as.integer(max_history_size), - as.double(init_alpha), as.double(tol_obj), as.double(tol_rel_obj), as.double(tol_grad), - as.double(tol_rel_grad), as.double(tol_param), as.integer(num_iterations), - as.integer(num_elbo_draws), as.integer(num_multi_draws), as.integer(calculate_lp), - as.integer(psis_resample), as.integer(refresh), as.integer(num_threads), - out = double(output_size), as.integer(output_size), err = raw(8), PACKAGE = model$lib_name) - handle_error(vars$return_code, model$lib_name, vars$err) - - list(draws = output_as_rvars(params, num_output, 1, vars$out)) - }) + params <- c(PATHFINDER_VARIABLES, get_parameter_names(model, model_ptr)) + num_params <- length(params) + output_size <- num_params * num_output + + vars <- .C( + "tinystan_pathfinder_R", + return_code = as.integer(0), + as.raw(model_ptr), + as.integer(num_paths), + encode_inits(model, inits), + as.integer(seed), + as.integer(id), + as.double(init_radius), + as.integer(num_draws), + as.integer(max_history_size), + as.double(init_alpha), + as.double(tol_obj), + as.double(tol_rel_obj), + as.double(tol_grad), + as.double(tol_rel_grad), + as.double(tol_param), + as.integer(num_iterations), + as.integer(num_elbo_draws), + as.integer(num_multi_draws), + as.integer(calculate_lp), + as.integer(psis_resample), + as.integer(refresh), + as.integer(num_threads), + out = double(output_size), + as.integer(output_size), + err = raw(8), + PACKAGE = model$lib_name + ) + handle_error(vars$return_code, model$lib_name, vars$err) + + list(draws = output_as_rvars(params, num_output, 1, vars$out)) + }) } #' @export optimizer <- function(...) { - UseMethod("optimizer") + UseMethod("optimizer") } #' @title Generic function `optimizer` #' @description Run Stan's Optimization algorithms #' @export -optimizer.tinystan_model = function(model, data = "", init = NULL, seed = NULL, id = 1, - init_radius = 2, algorithm = OptimizationAlgorithm$LBFGS, jacobian = FALSE, num_iterations = 2000, - max_history_size = 5, init_alpha = 0.001, tol_obj = 1e-12, tol_rel_obj = 10000, - tol_grad = 1e-08, tol_rel_grad = 1e+07, tol_param = 1e-08, refresh = 0, num_threads = -1) { - - if (is.null(seed)) { - seed <- as.integer(runif(1, min = 0, max = (2^31))) - } - - with_model(model, data, seed, { - params <- c(OPTIMIZATION_VARIABLES, get_parameter_names(model, model_ptr)) - num_params <- length(params) - output_size <- num_params - - vars <- .C("tinystan_optimize_R", return_code = as.integer(0), as.raw(model_ptr), - encode_inits(model, init), as.integer(seed), as.integer(id), as.double(init_radius), - as.integer(algorithm), as.integer(num_iterations), as.logical(jacobian), - as.integer(max_history_size), as.double(init_alpha), as.double(tol_obj), - as.double(tol_rel_obj), as.double(tol_grad), as.double(tol_rel_grad), - as.double(tol_param), as.integer(refresh), as.integer(num_threads), out = double(output_size), - as.integer(output_size), err = raw(8), PACKAGE = model$lib_name) - handle_error(vars$return_code, model$lib_name, vars$err) - - list(draws = output_as_rvars(params, 1, 1, vars$out)) - }) +optimizer.tinystan_model = function( + model, + data = "", + init = NULL, + seed = NULL, + id = 1, + init_radius = 2, + algorithm = OptimizationAlgorithm$LBFGS, + jacobian = FALSE, + num_iterations = 2000, + max_history_size = 5, + init_alpha = 0.001, + tol_obj = 1e-12, + tol_rel_obj = 10000, + tol_grad = 1e-08, + tol_rel_grad = 1e+07, + tol_param = 1e-08, + refresh = 0, + num_threads = -1 +) { + if (is.null(seed)) { + seed <- as.integer(runif(1, min = 0, max = (2^31))) + } + + with_model(model, data, seed, { + params <- c(OPTIMIZATION_VARIABLES, get_parameter_names(model, model_ptr)) + num_params <- length(params) + output_size <- num_params + + vars <- .C( + "tinystan_optimize_R", + return_code = as.integer(0), + as.raw(model_ptr), + encode_inits(model, init), + as.integer(seed), + as.integer(id), + as.double(init_radius), + as.integer(algorithm), + as.integer(num_iterations), + as.logical(jacobian), + as.integer(max_history_size), + as.double(init_alpha), + as.double(tol_obj), + as.double(tol_rel_obj), + as.double(tol_grad), + as.double(tol_rel_grad), + as.double(tol_param), + as.integer(refresh), + as.integer(num_threads), + out = double(output_size), + as.integer(output_size), + err = raw(8), + PACKAGE = model$lib_name + ) + handle_error(vars$return_code, model$lib_name, vars$err) + + list(draws = output_as_rvars(params, 1, 1, vars$out)) + }) } #' @export laplace_sampler <- function(...) { - UseMethod("laplace_sampler") + UseMethod("laplace_sampler") } #' @title Generic function `laplace_sampler` #' @description Run Stan's Laplace approximation algorithm #' @export -laplace_sampler.tinystan_model = function(model, mode, data = "", num_draws = 1000, - jacobian = TRUE, calculate_lp = TRUE, save_hessian = FALSE, seed = NULL, refresh = 0, - num_threads = -1) { - - if (num_draws < 1) { - stop("num_draws must be at least 1") - } - if (is.null(seed)) { - seed <- as.integer(runif(1, min = 0, max = (2^31))) +laplace_sampler.tinystan_model = function( + model, + mode, + data = "", + num_draws = 1000, + jacobian = TRUE, + calculate_lp = TRUE, + save_hessian = FALSE, + seed = NULL, + refresh = 0, + num_threads = -1 +) { + if (num_draws < 1) { + stop("num_draws must be at least 1") + } + if (is.null(seed)) { + seed <- as.integer(runif(1, min = 0, max = (2^31))) + } + + with_model(model, data, seed, { + req_params <- .C( + "tinystan_model_num_constrained_params_for_unconstraining_R", + as.raw(model_ptr), + params = as.integer(0), + PACKAGE = model$lib_name + )$params + + if (is.numeric(mode)) { + if (length(mode) < req_params) { + stop(paste0( + "Mode array has incorrect length.", + " Expected at least ", + req_params, + " elements." + )) + } + mode_array <- as.double(mode) + mode_json <- as.character("") + use_array <- TRUE + } else { + mode_array <- as.double(0) + mode_json <- as.character(mode) + use_array <- FALSE } - with_model(model, data, seed, { - req_params <- .C("tinystan_model_num_constrained_params_for_unconstraining_R", as.raw(model_ptr), params = as.integer(0), - PACKAGE = model$lib_name)$params - - if (is.numeric(mode)) { - if (length(mode) < req_params) { - stop(paste0("Mode array has incorrect length.", " Expected at least ", req_params, " elements.")) - } - mode_array <- as.double(mode) - mode_json <- as.character("") - use_array <- TRUE - } else { - mode_array <- as.double(0) - mode_json <- as.character(mode) - use_array <- FALSE - } + params <- c(LAPLACE_VARIABLES, get_parameter_names(model, model_ptr)) + num_params <- length(params) + free_params <- get_free_params(model, model_ptr) - params <- c(LAPLACE_VARIABLES, get_parameter_names(model, model_ptr)) - num_params <- length(params) - free_params <- get_free_params(model, model_ptr) - - if (save_hessian) { - hessian_size <- free_params * free_params - } else { - hessian_size <- 1 - } - - vars <- .C("tinystan_laplace_sample_R", return_code = as.integer(0), as.raw(model_ptr), - as.logical(use_array), mode_array, mode_json, as.integer(seed), as.integer(num_draws), - as.logical(jacobian), as.logical(calculate_lp), as.integer(refresh), - as.integer(num_threads), out = double(num_params * num_draws), as.integer(num_params * - num_draws), as.logical(save_hessian), hessian = double(hessian_size), - err = raw(8), PACKAGE = model$lib_name) - handle_error(vars$return_code, model$lib_name, vars$err) - - out <- output_as_rvars(params, num_draws, 1, vars$out) - if (save_hessian) { - hessian <- array(vars$hessian, dim = c(free_params, free_params)) - return(list(draws = out, hessian = hessian)) - } - list(draws = out) - }) + if (save_hessian) { + hessian_size <- free_params * free_params + } else { + hessian_size <- 1 + } + vars <- .C( + "tinystan_laplace_sample_R", + return_code = as.integer(0), + as.raw(model_ptr), + as.logical(use_array), + mode_array, + mode_json, + as.integer(seed), + as.integer(num_draws), + as.logical(jacobian), + as.logical(calculate_lp), + as.integer(refresh), + as.integer(num_threads), + out = double(num_params * num_draws), + as.integer( + num_params * + num_draws + ), + as.logical(save_hessian), + hessian = double(hessian_size), + err = raw(8), + PACKAGE = model$lib_name + ) + handle_error(vars$return_code, model$lib_name, vars$err) + + out <- output_as_rvars(params, num_draws, 1, vars$out) + if (save_hessian) { + hessian <- array(vars$hessian, dim = c(free_params, free_params)) + return(list(draws = out, hessian = hessian)) + } + list(draws = out) + }) } - #' Get and free the error message stored at the C++ pointer #' @noRd #' @keywords internal handle_error <- function(return_code, lib_name, err_ptr) { - if (return_code != 0) { - if (all(err_ptr == 0)) { - stop(paste("Unknown error, function returned code", return_code)) - } - msg <- .C("tinystan_get_error_message_R", as.raw(err_ptr), err_msg = as.character(""), - PACKAGE = lib_name)$err_msg - type <- .C("tinystan_get_error_type_R", as.raw(err_ptr), err_type = as.integer(0), - PACKAGE = lib_name)$err_type - .C("tinystan_destroy_error_R", as.raw(err_ptr), PACKAGE = lib_name) - if (type == 3) { - if (requireNamespace("rlang", quietly = TRUE)) { - rlang::interrupt() - } - msg <- "User interrupt" - } - stop(msg) + if (return_code != 0) { + if (all(err_ptr == 0)) { + stop(paste("Unknown error, function returned code", return_code)) + } + msg <- .C( + "tinystan_get_error_message_R", + as.raw(err_ptr), + err_msg = as.character(""), + PACKAGE = lib_name + )$err_msg + type <- .C( + "tinystan_get_error_type_R", + as.raw(err_ptr), + err_type = as.integer(0), + PACKAGE = lib_name + )$err_type + .C("tinystan_destroy_error_R", as.raw(err_ptr), PACKAGE = lib_name) + if (type == 3) { + if (requireNamespace("rlang", quietly = TRUE)) { + rlang::interrupt() + } + msg <- "User interrupt" } + stop(msg) + } } diff --git a/clients/R/R/zzz.R b/clients/R/R/zzz.R index ad7ec0e..90b7400 100644 --- a/clients/R/R/zzz.R +++ b/clients/R/R/zzz.R @@ -1,5 +1,12 @@ -HMC_SAMPLER_VARIABLES = c("lp__", "accept_stat__", "stepsize__", "treedepth__", "n_leapfrog__", - "divergent__", "energy__") +HMC_SAMPLER_VARIABLES = c( + "lp__", + "accept_stat__", + "stepsize__", + "treedepth__", + "n_leapfrog__", + "divergent__", + "energy__" +) PATHFINDER_VARIABLES = c("lp_approx__", "lp__", "path__") diff --git a/clients/R/convert_docs.R b/clients/R/convert_docs.R index feb533d..6fe4e43 100644 --- a/clients/R/convert_docs.R +++ b/clients/R/convert_docs.R @@ -9,14 +9,15 @@ roxygen2::roxygenize() files <- list.files("man", pattern = "*.Rd") for (f in files) { - name <- substr(f, 1, nchar(f) - 3) + name <- substr(f, 1, nchar(f) - 3) - # read .Rd file and convert to markdown - rd <- rd2markdown::get_rd(file = file.path(".", "man", f)) - md <- rd2markdown::rd2markdown(rd, fragments = c(), level = 3) - - # write it to the docs folder - writeLines(md, file.path("..", "..", "docs", "languages", "_r", paste0(name, - ".md"))) + # read .Rd file and convert to markdown + rd <- rd2markdown::get_rd(file = file.path(".", "man", f)) + md <- rd2markdown::rd2markdown(rd, fragments = c(), level = 3) + # write it to the docs folder + writeLines( + md, + file.path("..", "..", "docs", "languages", "_r", paste0(name, ".md")) + ) } diff --git a/clients/R/tests/testthat.R b/clients/R/tests/testthat.R index 4f39674..a12d757 100644 --- a/clients/R/tests/testthat.R +++ b/clients/R/tests/testthat.R @@ -1,4 +1,3 @@ - library(testthat) library(tinystan) diff --git a/clients/R/tests/testthat/setup.R b/clients/R/tests/testthat/setup.R index 474a6c0..d04f0a4 100644 --- a/clients/R/tests/testthat/setup.R +++ b/clients/R/tests/testthat/setup.R @@ -7,22 +7,38 @@ top_level_directory <- file.path("../../../../") set_tinystan_path(top_level_directory) stan_folder <- file.path(top_level_directory, "test_models") -bernoulli_model <- tinystan_model(file.path(stan_folder, "bernoulli", "bernoulli.stan")) +bernoulli_model <- tinystan_model(file.path( + stan_folder, + "bernoulli", + "bernoulli.stan" +)) BERNOULLI_DATA <- "{\"N\": 10, \"y\": [0, 1, 0, 0, 0, 0, 0, 0, 0, 1]}" -gaussian_model <- tinystan_model(file.path(stan_folder, "gaussian", "gaussian.stan")) +gaussian_model <- tinystan_model(file.path( + stan_folder, + "gaussian", + "gaussian.stan" +)) empty_model <- tinystan_model(file.path(stan_folder, "empty", "empty.stan")) -multimodal_model <- tinystan_model(file.path(stan_folder, "multimodal", "multimodal.stan")) -simple_jacobian_model <- tinystan_model(file.path(stan_folder, "simple_jacobian", - "simple_jacobian.stan")) +multimodal_model <- tinystan_model(file.path( + stan_folder, + "multimodal", + "multimodal.stan" +)) +simple_jacobian_model <- tinystan_model(file.path( + stan_folder, + "simple_jacobian", + "simple_jacobian.stan" +)) # hack around the fact that comparisons to NULL result in logical(0) and # all(logical(0)) is TRUE, for some reason. .builtin_all <- all all <- function(l) { - if (length(l) == 0) - return(isTRUE(l)) - .builtin_all(l) + if (length(l) == 0) { + return(isTRUE(l)) + } + .builtin_all(l) } options(warn = 2) diff --git a/clients/R/tests/testthat/test_compile.R b/clients/R/tests/testthat/test_compile.R index 07e8633..fd6eab5 100644 --- a/clients/R/tests/testthat/test_compile.R +++ b/clients/R/tests/testthat/test_compile.R @@ -1,35 +1,42 @@ test_that("compilation works", { - name <- "gaussian" - file <- file.path(stan_folder, name, paste0(name, ".stan")) + name <- "gaussian" + file <- file.path(stan_folder, name, paste0(name, ".stan")) - lib <- file.path(stan_folder, name, paste0(name, "_model.so")) - unlink(lib, force = TRUE) + lib <- file.path(stan_folder, name, paste0(name, "_model.so")) + unlink(lib, force = TRUE) - out <- compile_model(file, stanc_args = c("--O1")) + out <- compile_model(file, stanc_args = c("--O1")) - expect_true(file.exists(lib)) - expect_equal(normalizePath(lib), normalizePath(out)) + expect_true(file.exists(lib)) + expect_equal(normalizePath(lib), normalizePath(out)) - unlink(lib, force = TRUE) + unlink(lib, force = TRUE) - out <- compile_model(file) + out <- compile_model(file) }) test_that("compilation fails on non-stan file", { - expect_error(compile_model(file.path(stan_folder, "bernoulli", "bernoulli.data.json")), - "does not end with '.stan'") + expect_error( + compile_model(file.path(stan_folder, "bernoulli", "bernoulli.data.json")), + "does not end with '.stan'" + ) }) test_that("compilation fails on missing file", { - expect_error(compile_model("badpath.stan"), "does not exist!") + expect_error(compile_model("badpath.stan"), "does not exist!") }) test_that("compilation fails on bad syntax", { - expect_error(compile_model(file.path(stan_folder, "syntax_error", "syntax_error.stan")), - "Compilation failed") + expect_error( + compile_model(file.path(stan_folder, "syntax_error", "syntax_error.stan")), + "Compilation failed" + ) }) test_that("bad paths fail", { - expect_error(set_tinystan_path("badpath"), "does not exist!") - expect_error(set_tinystan_path(file.path(stan_folder)), "does not contain file 'Makefile'") + expect_error(set_tinystan_path("badpath"), "does not exist!") + expect_error( + set_tinystan_path(file.path(stan_folder)), + "does not contain file 'Makefile'" + ) }) diff --git a/clients/R/tests/testthat/test_laplace.R b/clients/R/tests/testthat/test_laplace.R index 3cded20..6bcf335 100644 --- a/clients/R/tests/testthat/test_laplace.R +++ b/clients/R/tests/testthat/test_laplace.R @@ -1,108 +1,167 @@ BERNOULLI_MODE <- "{\"theta\": 0.25}" test_that("data arguments work", { + out1 <- laplace_sampler(bernoulli_model, BERNOULLI_MODE, BERNOULLI_DATA) + expect_true(mean(out1$draws$theta) > 0.22 && mean(out1$draws$theta) < 0.28) - out1 <- laplace_sampler(bernoulli_model, BERNOULLI_MODE, BERNOULLI_DATA) - expect_true(mean(out1$draws$theta) > 0.22 && mean(out1$draws$theta) < 0.28) - - data_file <- file.path(stan_folder, "bernoulli", "bernoulli.data.json") - out2 <- laplace_sampler(bernoulli_model, BERNOULLI_MODE, data = data_file) - expect_true(mean(out2$draws$theta) > 0.22 && mean(out2$draws$theta) < 0.28) + data_file <- file.path(stan_folder, "bernoulli", "bernoulli.data.json") + out2 <- laplace_sampler(bernoulli_model, BERNOULLI_MODE, data = data_file) + expect_true(mean(out2$draws$theta) > 0.22 && mean(out2$draws$theta) < 0.28) }) test_that("output sizes are correct", { - out1 <- laplace_sampler(bernoulli_model, BERNOULLI_MODE, BERNOULLI_DATA, num_draws = 324) - expect_equal(posterior::ndraws(out1$draws), 324) + out1 <- laplace_sampler( + bernoulli_model, + BERNOULLI_MODE, + BERNOULLI_DATA, + num_draws = 324 + ) + expect_equal(posterior::ndraws(out1$draws), 324) }) test_that("calculate_lp works", { - out <- laplace_sampler(bernoulli_model, BERNOULLI_MODE, BERNOULLI_DATA, num_draws = 500, - calculate_lp = TRUE) - expect_equal(sum(is.nan(posterior::draws_of(out$draws$log_p__))), 0) - - out2 <- laplace_sampler(bernoulli_model, BERNOULLI_MODE, BERNOULLI_DATA, num_draws = 500, - calculate_lp = FALSE) - expect_equal(sum(is.nan(posterior::draws_of(out2$draws$log_p__))), 500) + out <- laplace_sampler( + bernoulli_model, + BERNOULLI_MODE, + BERNOULLI_DATA, + num_draws = 500, + calculate_lp = TRUE + ) + expect_equal(sum(is.nan(posterior::draws_of(out$draws$log_p__))), 0) + + out2 <- laplace_sampler( + bernoulli_model, + BERNOULLI_MODE, + BERNOULLI_DATA, + num_draws = 500, + calculate_lp = FALSE + ) + expect_equal(sum(is.nan(posterior::draws_of(out2$draws$log_p__))), 500) }) test_that("jacobian arg works", { - for (jacobian in c(TRUE, FALSE)) { - - out <- optimizer(simple_jacobian_model, jacobian = jacobian, seed = 12345) - sigma <- posterior::extract_variable(out$draws, "sigma") - - draws <- laplace_sampler(simple_jacobian_model, c(sigma), jacobian = jacobian, - seed = 12345) - sigma <- mean(posterior::extract_variable(draws$draws, "sigma")) - if (jacobian) { - expect_equal(sigma, 3.3, tolerance = 0.2, ignore_attr = TRUE) - } else { - expect_equal(sigma, 3, tolerance = 0.2, ignore_attr = TRUE) - } - + for (jacobian in c(TRUE, FALSE)) { + out <- optimizer(simple_jacobian_model, jacobian = jacobian, seed = 12345) + sigma <- posterior::extract_variable(out$draws, "sigma") + + draws <- laplace_sampler( + simple_jacobian_model, + c(sigma), + jacobian = jacobian, + seed = 12345 + ) + sigma <- mean(posterior::extract_variable(draws$draws, "sigma")) + if (jacobian) { + expect_equal(sigma, 3.3, tolerance = 0.2, ignore_attr = TRUE) + } else { + expect_equal(sigma, 3, tolerance = 0.2, ignore_attr = TRUE) } + } }) test_that("save_hessian works", { - data <- "{\"N\": 3}" - mode <- "{\"alpha\": [0.1,0.2,0.3]}" + data <- "{\"N\": 3}" + mode <- "{\"alpha\": [0.1,0.2,0.3]}" - out <- laplace_sampler(gaussian_model, mode, data, save_hessian = TRUE) - expect_true("hessian" %in% names(out)) - expect_equal(dim(out$hessian), c(3, 3)) - expect_equal(out$hessian, matrix(c(-1, 0, 0, 0, -1, 0, 0, 0, -1), nrow = 3)) + out <- laplace_sampler(gaussian_model, mode, data, save_hessian = TRUE) + expect_true("hessian" %in% names(out)) + expect_equal(dim(out$hessian), c(3, 3)) + expect_equal(out$hessian, matrix(c(-1, 0, 0, 0, -1, 0, 0, 0, -1), nrow = 3)) }) test_that("seed works", { - - out1 <- laplace_sampler(bernoulli_model, BERNOULLI_MODE, BERNOULLI_DATA, seed = 123) - out2 <- laplace_sampler(bernoulli_model, BERNOULLI_MODE, BERNOULLI_DATA, seed = 123) - - expect_equal(out1$draws$theta, out2$draws$theta) - - out3 <- laplace_sampler(bernoulli_model, BERNOULLI_MODE, BERNOULLI_DATA, seed = 456) - expect_error(expect_equal(out1$draws$theta, out3$draws$theta)) - + out1 <- laplace_sampler( + bernoulli_model, + BERNOULLI_MODE, + BERNOULLI_DATA, + seed = 123 + ) + out2 <- laplace_sampler( + bernoulli_model, + BERNOULLI_MODE, + BERNOULLI_DATA, + seed = 123 + ) + + expect_equal(out1$draws$theta, out2$draws$theta) + + out3 <- laplace_sampler( + bernoulli_model, + BERNOULLI_MODE, + BERNOULLI_DATA, + seed = 456 + ) + expect_error(expect_equal(out1$draws$theta, out3$draws$theta)) }) test_that("bad data handled properly", { - - data <- "{\"N\": -1}" - expect_error(laplace_sampler(bernoulli_model, BERNOULLI_MODE, data), "greater than or equal to 0") - - data <- "{\"N\": 1, \"y\": [1,2]}" - expect_error(laplace_sampler(bernoulli_model, BERNOULLI_MODE, data), "mismatch in dimension") - - expect_error(laplace_sampler(bernoulli_model, BERNOULLI_MODE, "{\"bad\"}"), "Error in JSON parsing") - - expect_error(laplace_sampler(bernoulli_model, BERNOULLI_MODE, "not/real/path.json"), - "Could not open data file") - + data <- "{\"N\": -1}" + expect_error( + laplace_sampler(bernoulli_model, BERNOULLI_MODE, data), + "greater than or equal to 0" + ) + + data <- "{\"N\": 1, \"y\": [1,2]}" + expect_error( + laplace_sampler(bernoulli_model, BERNOULLI_MODE, data), + "mismatch in dimension" + ) + + expect_error( + laplace_sampler(bernoulli_model, BERNOULLI_MODE, "{\"bad\"}"), + "Error in JSON parsing" + ) + + expect_error( + laplace_sampler(bernoulli_model, BERNOULLI_MODE, "not/real/path.json"), + "Could not open data file" + ) }) test_that("bad mode array handled properly", { - mode1 = c(2) - expect_error(laplace_sampler(bernoulli_model, mode1, BERNOULLI_DATA), "Bounded variable is 2") - - mode2 = c(1.0) - expect_error(laplace_sampler(gaussian_model, mode2, "{\"N\": 4 }"), "incorrect length") + mode1 = c(2) + expect_error( + laplace_sampler(bernoulli_model, mode1, BERNOULLI_DATA), + "Bounded variable is 2" + ) + + mode2 = c(1.0) + expect_error( + laplace_sampler(gaussian_model, mode2, "{\"N\": 4 }"), + "incorrect length" + ) }) test_that("bad mode json handled properly", { - mode <- "{\"theta\": 2}" - expect_error(laplace_sampler(bernoulli_model, mode, BERNOULLI_DATA), "Bounded variable is 2") - - mode <- "{\"theta\": [0.5, 0.5]}" - expect_error(laplace_sampler(bernoulli_model, mode, BERNOULLI_DATA), "mismatch in number") - - expect_error(laplace_sampler(bernoulli_model, "bad/path.json", BERNOULLI_DATA), - "Could not open data file") + mode <- "{\"theta\": 2}" + expect_error( + laplace_sampler(bernoulli_model, mode, BERNOULLI_DATA), + "Bounded variable is 2" + ) + + mode <- "{\"theta\": [0.5, 0.5]}" + expect_error( + laplace_sampler(bernoulli_model, mode, BERNOULLI_DATA), + "mismatch in number" + ) + + expect_error( + laplace_sampler(bernoulli_model, "bad/path.json", BERNOULLI_DATA), + "Could not open data file" + ) }) test_that("bad num_draws raises errors", { - expect_error(laplace_sampler(bernoulli_model, BERNOULLI_MODE, BERNOULLI_DATA, - num_draws = 0), "at least 1") + expect_error( + laplace_sampler( + bernoulli_model, + BERNOULLI_MODE, + BERNOULLI_DATA, + num_draws = 0 + ), + "at least 1" + ) }) diff --git a/clients/R/tests/testthat/test_model.R b/clients/R/tests/testthat/test_model.R index c870980..e49f9e2 100644 --- a/clients/R/tests/testthat/test_model.R +++ b/clients/R/tests/testthat/test_model.R @@ -1,17 +1,16 @@ - bernoulli_file <- file.path(stan_folder, "bernoulli", "bernoulli.stan") test_that("model loads", { - suppressWarnings(model <- tinystan_model(bernoulli_file)) - expect_true(!is.null(model)) + suppressWarnings(model <- tinystan_model(bernoulli_file)) + expect_true(!is.null(model)) }) test_that("api_version is correct", { - expect_equal(api_version(bernoulli_model), current_version_list) + expect_equal(api_version(bernoulli_model), current_version_list) }) test_that("stan version is valid", { - stan_version <- stan_version(bernoulli_model) - expect_equal(stan_version$major, 2) - expect_gte(stan_version$minor, 34) - expect_gte(stan_version$patch, 0) + stan_version <- stan_version(bernoulli_model) + expect_equal(stan_version$major, 2) + expect_gte(stan_version$minor, 34) + expect_gte(stan_version$patch, 0) }) diff --git a/clients/R/tests/testthat/test_optimize.R b/clients/R/tests/testthat/test_optimize.R index f075dc3..273caf9 100644 --- a/clients/R/tests/testthat/test_optimize.R +++ b/clients/R/tests/testthat/test_optimize.R @@ -1,141 +1,199 @@ -ALGORITHMS <- c(tinystan::OptimizationAlgorithm$NEWTON, tinystan::OptimizationAlgorithm$BFGS, - tinystan::OptimizationAlgorithm$LBFGS) +ALGORITHMS <- c( + tinystan::OptimizationAlgorithm$NEWTON, + tinystan::OptimizationAlgorithm$BFGS, + tinystan::OptimizationAlgorithm$LBFGS +) test_that("data args work", { + out1 <- optimizer(bernoulli_model, BERNOULLI_DATA) + expect_true(mean(out1$draws$theta) > 0.19 && mean(out1$draws$theta) < 0.21) - out1 <- optimizer(bernoulli_model, BERNOULLI_DATA) - expect_true(mean(out1$draws$theta) > 0.19 && mean(out1$draws$theta) < 0.21) - - data_file <- file.path(stan_folder, "bernoulli", "bernoulli.data.json") - out2 <- optimizer(bernoulli_model, data = data_file) - expect_true(mean(out2$draws$theta) > 0.19 && mean(out2$draws$theta) < 0.21) - + data_file <- file.path(stan_folder, "bernoulli", "bernoulli.data.json") + out2 <- optimizer(bernoulli_model, data = data_file) + expect_true(mean(out2$draws$theta) > 0.19 && mean(out2$draws$theta) < 0.21) }) test_that("algorithm and jacobian args work", { - - for (algorithm in ALGORITHMS) { - for (jacobian in c(TRUE, FALSE)) { - - out <- optimizer(simple_jacobian_model, algorithm = algorithm, jacobian = jacobian, - seed = 1234) - - sigma <- posterior::extract_variable(out$draws, "sigma") - - if (jacobian) { - expect_equal(sigma, 3.3, tolerance = 0.01, ignore_attr = TRUE) - } else { - expect_equal(sigma, 3, tolerance = 0.01, ignore_attr = TRUE) - } - - } + for (algorithm in ALGORITHMS) { + for (jacobian in c(TRUE, FALSE)) { + out <- optimizer( + simple_jacobian_model, + algorithm = algorithm, + jacobian = jacobian, + seed = 1234 + ) + + sigma <- posterior::extract_variable(out$draws, "sigma") + + if (jacobian) { + expect_equal(sigma, 3.3, tolerance = 0.01, ignore_attr = TRUE) + } else { + expect_equal(sigma, 3, tolerance = 0.01, ignore_attr = TRUE) + } } - + } }) test_that("seed works", { + out1 <- optimizer(bernoulli_model, BERNOULLI_DATA, seed = 123) + out2 <- optimizer(bernoulli_model, BERNOULLI_DATA, seed = 123) - out1 <- optimizer(bernoulli_model, BERNOULLI_DATA, seed = 123) - out2 <- optimizer(bernoulli_model, BERNOULLI_DATA, seed = 123) - - expect_equal(out1$draws, out2$draws) - - out3 <- optimizer(bernoulli_model, BERNOULLI_DATA, seed = 456) - expect_error(expect_equal(out1$draws, out3$draws)) + expect_equal(out1$draws, out2$draws) + out3 <- optimizer(bernoulli_model, BERNOULLI_DATA, seed = 456) + expect_error(expect_equal(out1$draws, out3$draws)) }) test_that("inits work", { + init <- "{\"mu\": -100}" + out1 <- optimizer(multimodal_model, init = init) + expect_true(all(out1$draws$mu < 0)) - init <- "{\"mu\": -100}" - out1 <- optimizer(multimodal_model, init = init) - expect_true(all(out1$draws$mu < 0)) - - init <- "{\"mu\": 100}" - temp_file <- tempfile(fileext = ".json") - write(init, temp_file) - - out2 <- optimizer(multimodal_model, init = temp_file) - expect_true(all(out2$draws$mu > 0)) + init <- "{\"mu\": 100}" + temp_file <- tempfile(fileext = ".json") + write(init, temp_file) + out2 <- optimizer(multimodal_model, init = temp_file) + expect_true(all(out2$draws$mu > 0)) }) test_that("bad data handled properly", { + data <- "{\"N\": -1}" + expect_error(optimizer(bernoulli_model, data), "greater than or equal to 0") - data <- "{\"N\": -1}" - expect_error(optimizer(bernoulli_model, data), "greater than or equal to 0") + data <- "{\"N\": 1, \"y\": [1,2]}" + expect_error(optimizer(bernoulli_model, data), "mismatch in dimension") - data <- "{\"N\": 1, \"y\": [1,2]}" - expect_error(optimizer(bernoulli_model, data), "mismatch in dimension") - - expect_error(optimizer(bernoulli_model, "{\"bad\"}"), "Error in JSON parsing") - - expect_error(optimizer(bernoulli_model, "not/real/path.json"), "Could not open data file") + expect_error(optimizer(bernoulli_model, "{\"bad\"}"), "Error in JSON parsing") + expect_error( + optimizer(bernoulli_model, "not/real/path.json"), + "Could not open data file" + ) }) test_that("bad init handled properly", { - - init <- "{\"theta\": 2}" - expect_error(optimizer(bernoulli_model, BERNOULLI_DATA, init = init), "Initialization failed") - - expect_error(optimizer(bernoulli_model, BERNOULLI_DATA, init = "bad/path.json"), - "Could not open data file") - + init <- "{\"theta\": 2}" + expect_error( + optimizer(bernoulli_model, BERNOULLI_DATA, init = init), + "Initialization failed" + ) + + expect_error( + optimizer(bernoulli_model, BERNOULLI_DATA, init = "bad/path.json"), + "Could not open data file" + ) }) test_that("empty model ok", { - - expect_no_error(optimizer(empty_model)) - + expect_no_error(optimizer(empty_model)) }) test_that("bad args raise errors", { - - for (algorithm in ALGORITHMS) { - - expect_error(optimizer(bernoulli_model, BERNOULLI_DATA, algorithm = algorithm, - id = 0), "positive") - expect_error(optimizer(bernoulli_model, BERNOULLI_DATA, algorithm = algorithm, - num_iterations = 0), "positive") - expect_error(optimizer(bernoulli_model, BERNOULLI_DATA, algorithm = algorithm, - init_radius = -0.1), "non-negative") - - if (algorithm != tinystan::OptimizationAlgorithm$NEWTON) { - expected <- expect_error - } else { - expected <- function(e, ...) { - expect_no_error(e) - } - } - - expected(optimizer(bernoulli_model, BERNOULLI_DATA, algorithm = algorithm, - init_alpha = 0), "positive") - expected(optimizer(bernoulli_model, BERNOULLI_DATA, algorithm = algorithm, - tol_obj = 0), "positive") - expected(optimizer(bernoulli_model, BERNOULLI_DATA, algorithm = algorithm, - tol_rel_obj = 0), "positive") - expected(optimizer(bernoulli_model, BERNOULLI_DATA, algorithm = algorithm, - tol_grad = 0), "positive") - expected(optimizer(bernoulli_model, BERNOULLI_DATA, algorithm = algorithm, - tol_rel_grad = 0), "positive") - expected(optimizer(bernoulli_model, BERNOULLI_DATA, algorithm = algorithm, - tol_param = 0), "positive") - - if (algorithm == tinystan::OptimizationAlgorithm$LBFGS) { - expected <- expect_error - } else { - expected <- function(e, ...) { - expect_no_error(e) - } - } - - expected(optimizer(bernoulli_model, BERNOULLI_DATA, algorithm = algorithm, - max_history_size = 0), "positive") - - + for (algorithm in ALGORITHMS) { + expect_error( + optimizer(bernoulli_model, BERNOULLI_DATA, algorithm = algorithm, id = 0), + "positive" + ) + expect_error( + optimizer( + bernoulli_model, + BERNOULLI_DATA, + algorithm = algorithm, + num_iterations = 0 + ), + "positive" + ) + expect_error( + optimizer( + bernoulli_model, + BERNOULLI_DATA, + algorithm = algorithm, + init_radius = -0.1 + ), + "non-negative" + ) + + if (algorithm != tinystan::OptimizationAlgorithm$NEWTON) { + expected <- expect_error + } else { + expected <- function(e, ...) { + expect_no_error(e) + } } + expected( + optimizer( + bernoulli_model, + BERNOULLI_DATA, + algorithm = algorithm, + init_alpha = 0 + ), + "positive" + ) + expected( + optimizer( + bernoulli_model, + BERNOULLI_DATA, + algorithm = algorithm, + tol_obj = 0 + ), + "positive" + ) + expected( + optimizer( + bernoulli_model, + BERNOULLI_DATA, + algorithm = algorithm, + tol_rel_obj = 0 + ), + "positive" + ) + expected( + optimizer( + bernoulli_model, + BERNOULLI_DATA, + algorithm = algorithm, + tol_grad = 0 + ), + "positive" + ) + expected( + optimizer( + bernoulli_model, + BERNOULLI_DATA, + algorithm = algorithm, + tol_rel_grad = 0 + ), + "positive" + ) + expected( + optimizer( + bernoulli_model, + BERNOULLI_DATA, + algorithm = algorithm, + tol_param = 0 + ), + "positive" + ) + + if (algorithm == tinystan::OptimizationAlgorithm$LBFGS) { + expected <- expect_error + } else { + expected <- function(e, ...) { + expect_no_error(e) + } + } + expected( + optimizer( + bernoulli_model, + BERNOULLI_DATA, + algorithm = algorithm, + max_history_size = 0 + ), + "positive" + ) + } }) diff --git a/clients/R/tests/testthat/test_pathfinder.R b/clients/R/tests/testthat/test_pathfinder.R index e3aec83..0c0b491 100644 --- a/clients/R/tests/testthat/test_pathfinder.R +++ b/clients/R/tests/testthat/test_pathfinder.R @@ -1,144 +1,242 @@ test_that("data arguments work", { + out1 <- pathfinder(bernoulli_model, BERNOULLI_DATA) + expect_true(mean(out1$draws$theta) > 0.2 && mean(out1$draws$theta) < 0.3) - out1 <- pathfinder(bernoulli_model, BERNOULLI_DATA) - expect_true(mean(out1$draws$theta) > 0.2 && mean(out1$draws$theta) < 0.3) - - data_file <- file.path(stan_folder, "bernoulli", "bernoulli.data.json") - out2 <- pathfinder(bernoulli_model, data = data_file) - expect_true(mean(out2$draws$theta) > 0.2 && mean(out2$draws$theta) < 0.3) - + data_file <- file.path(stan_folder, "bernoulli", "bernoulli.data.json") + out2 <- pathfinder(bernoulli_model, data = data_file) + expect_true(mean(out2$draws$theta) > 0.2 && mean(out2$draws$theta) < 0.3) }) test_that("output sizes are correct", { - - out1 <- pathfinder(bernoulli_model, BERNOULLI_DATA, num_paths = 4, num_draws = 101, - num_multi_draws = 99) - expect_equal(posterior::ndraws(out1$draws), 99) - - out2 <- pathfinder(bernoulli_model, BERNOULLI_DATA, num_paths = 1, num_draws = 101, - num_multi_draws = 103) - expect_equal(posterior::ndraws(out2$draws), 103) - - out3 <- pathfinder(bernoulli_model, BERNOULLI_DATA, num_paths = 2, num_draws = 105, - num_multi_draws = 1, calculate_lp = FALSE) - expect_equal(posterior::ndraws(out3$draws), 2 * 105) - - out4 <- pathfinder(bernoulli_model, BERNOULLI_DATA, num_paths = 3, num_draws = 107, - num_multi_draws = 1, psis_resample = FALSE) - expect_equal(posterior::ndraws(out4$draws), 3 * 107) - - out5 <- pathfinder(bernoulli_model, BERNOULLI_DATA, num_paths = 1, num_draws = 109, - num_multi_draws = 1, psis_resample = FALSE) - expect_equal(posterior::ndraws(out5$draws), 109) + out1 <- pathfinder( + bernoulli_model, + BERNOULLI_DATA, + num_paths = 4, + num_draws = 101, + num_multi_draws = 99 + ) + expect_equal(posterior::ndraws(out1$draws), 99) + + out2 <- pathfinder( + bernoulli_model, + BERNOULLI_DATA, + num_paths = 1, + num_draws = 101, + num_multi_draws = 103 + ) + expect_equal(posterior::ndraws(out2$draws), 103) + + out3 <- pathfinder( + bernoulli_model, + BERNOULLI_DATA, + num_paths = 2, + num_draws = 105, + num_multi_draws = 1, + calculate_lp = FALSE + ) + expect_equal(posterior::ndraws(out3$draws), 2 * 105) + + out4 <- pathfinder( + bernoulli_model, + BERNOULLI_DATA, + num_paths = 3, + num_draws = 107, + num_multi_draws = 1, + psis_resample = FALSE + ) + expect_equal(posterior::ndraws(out4$draws), 3 * 107) + + out5 <- pathfinder( + bernoulli_model, + BERNOULLI_DATA, + num_paths = 1, + num_draws = 109, + num_multi_draws = 1, + psis_resample = FALSE + ) + expect_equal(posterior::ndraws(out5$draws), 109) }) test_that("calculate_lp works", { - out <- pathfinder(bernoulli_model, BERNOULLI_DATA, num_paths = 2, calculate_lp = FALSE) - - expect_gt(sum(is.nan(out$draws$lp__)), 0) - expect_lt(sum(is.nan(out$draws$lp__)), 2000) - - out_single <- pathfinder(bernoulli_model, BERNOULLI_DATA, num_paths = 1, calculate_lp = FALSE) - - expect_gt(sum(is.nan(out_single$draws$lp__)), 0) - expect_lt(sum(is.nan(out_single$draws$lp__)), 1000) + out <- pathfinder( + bernoulli_model, + BERNOULLI_DATA, + num_paths = 2, + calculate_lp = FALSE + ) + + expect_gt(sum(is.nan(out$draws$lp__)), 0) + expect_lt(sum(is.nan(out$draws$lp__)), 2000) + + out_single <- pathfinder( + bernoulli_model, + BERNOULLI_DATA, + num_paths = 1, + calculate_lp = FALSE + ) + + expect_gt(sum(is.nan(out_single$draws$lp__)), 0) + expect_lt(sum(is.nan(out_single$draws$lp__)), 1000) }) test_that("seed works", { - - out1 <- pathfinder(bernoulli_model, BERNOULLI_DATA, seed = 123) - out2 <- pathfinder(bernoulli_model, BERNOULLI_DATA, seed = 123) - - expect_equal(sort(posterior::draws_of(out1$draws$theta)), sort(posterior::draws_of(out2$draws$theta))) - - out3 <- pathfinder(bernoulli_model, BERNOULLI_DATA, seed = 456) - expect_error(expect_equal(sort(posterior::draws_of(out1$draws$theta)), sort(posterior::draws_of(out3$draws$theta)))) - + out1 <- pathfinder(bernoulli_model, BERNOULLI_DATA, seed = 123) + out2 <- pathfinder(bernoulli_model, BERNOULLI_DATA, seed = 123) + + expect_equal( + sort(posterior::draws_of(out1$draws$theta)), + sort(posterior::draws_of(out2$draws$theta)) + ) + + out3 <- pathfinder(bernoulli_model, BERNOULLI_DATA, seed = 456) + expect_error(expect_equal( + sort(posterior::draws_of(out1$draws$theta)), + sort(posterior::draws_of(out3$draws$theta)) + )) }) test_that("inits work", { - - init1 <- "{\"mu\": -1000}" - out1 <- pathfinder(multimodal_model, inits = init1) - expect_true(all(out1$draws$mu < 0)) - - init2 <- "{\"mu\": 1000}" - out2 <- pathfinder(multimodal_model, inits = init2) - expect_true(all(out2$draws$mu > 0)) - - temp_file <- tempfile(fileext = ".json") - write(init1, temp_file) - out3 <- pathfinder(multimodal_model, num_paths = 2, inits = c(temp_file, init1)) - - expect_true(all(out3$draws$mu < 0)) + init1 <- "{\"mu\": -1000}" + out1 <- pathfinder(multimodal_model, inits = init1) + expect_true(all(out1$draws$mu < 0)) + + init2 <- "{\"mu\": 1000}" + out2 <- pathfinder(multimodal_model, inits = init2) + expect_true(all(out2$draws$mu > 0)) + + temp_file <- tempfile(fileext = ".json") + write(init1, temp_file) + out3 <- pathfinder( + multimodal_model, + num_paths = 2, + inits = c(temp_file, init1) + ) + + expect_true(all(out3$draws$mu < 0)) }) test_that("bad data handled properly", { + data <- "{\"N\": -1}" + expect_error(pathfinder(bernoulli_model, data), "greater than or equal to 0") - data <- "{\"N\": -1}" - expect_error(pathfinder(bernoulli_model, data), "greater than or equal to 0") + data <- "{\"N\": 1, \"y\": [1,2]}" + expect_error(pathfinder(bernoulli_model, data), "mismatch in dimension") - data <- "{\"N\": 1, \"y\": [1,2]}" - expect_error(pathfinder(bernoulli_model, data), "mismatch in dimension") - - expect_error(pathfinder(bernoulli_model, "{\"bad\"}"), "Error in JSON parsing") - - expect_error(pathfinder(bernoulli_model, "not/real/path.json"), "Could not open data file") + expect_error( + pathfinder(bernoulli_model, "{\"bad\"}"), + "Error in JSON parsing" + ) + expect_error( + pathfinder(bernoulli_model, "not/real/path.json"), + "Could not open data file" + ) }) test_that("bad inits handled properly", { - init <- "{\"theta\": 2}" - expect_error(pathfinder(bernoulli_model, BERNOULLI_DATA, inits = init), "Initialization failed") - expect_error(pathfinder(bernoulli_model, BERNOULLI_DATA, num_paths = 1, inits = init), - "Initialization failed") - expect_error(pathfinder(bernoulli_model, BERNOULLI_DATA, inits = "bad/path.json"), - "Could not open data file") - inits <- c(init, init) - expect_error(pathfinder(bernoulli_model, BERNOULLI_DATA, num_paths = 2, inits = inits), - "Initialization failed") - - init2 <- "{\"theta\": 0.2}" - - # unlike sample, a failure of subset of inits is not fatal - inits <- rep(list(init), 10) - inits[[11]] <- init2 - pathfinder(bernoulli_model, BERNOULLI_DATA, num_paths = 11, num_multi_draws = 10, - inits = inits) - inits <- list(init2, init2) - expect_error(pathfinder(bernoulli_model, BERNOULLI_DATA, num_paths = 1, inits = inits), - "match the number of chains") - expect_error(pathfinder(bernoulli_model, BERNOULLI_DATA, num_paths = 3, inits = inits), - "match the number of chains") - + init <- "{\"theta\": 2}" + expect_error( + pathfinder(bernoulli_model, BERNOULLI_DATA, inits = init), + "Initialization failed" + ) + expect_error( + pathfinder(bernoulli_model, BERNOULLI_DATA, num_paths = 1, inits = init), + "Initialization failed" + ) + expect_error( + pathfinder(bernoulli_model, BERNOULLI_DATA, inits = "bad/path.json"), + "Could not open data file" + ) + inits <- c(init, init) + expect_error( + pathfinder(bernoulli_model, BERNOULLI_DATA, num_paths = 2, inits = inits), + "Initialization failed" + ) + + init2 <- "{\"theta\": 0.2}" + + # unlike sample, a failure of subset of inits is not fatal + inits <- rep(list(init), 10) + inits[[11]] <- init2 + pathfinder( + bernoulli_model, + BERNOULLI_DATA, + num_paths = 11, + num_multi_draws = 10, + inits = inits + ) + inits <- list(init2, init2) + expect_error( + pathfinder(bernoulli_model, BERNOULLI_DATA, num_paths = 1, inits = inits), + "match the number of chains" + ) + expect_error( + pathfinder(bernoulli_model, BERNOULLI_DATA, num_paths = 3, inits = inits), + "match the number of chains" + ) }) test_that("empty model ok", { - expect_error(pathfinder(empty_model), "no parameters") + expect_error(pathfinder(empty_model), "no parameters") }) test_that("bad args raise errors", { - expect_error(pathfinder(bernoulli_model, BERNOULLI_DATA, num_paths = 0), "at least 1") - expect_error(pathfinder(bernoulli_model, BERNOULLI_DATA, num_draws = 0), "at least 1") - expect_error(pathfinder(bernoulli_model, BERNOULLI_DATA, id = 0), "positive") - expect_error(pathfinder(bernoulli_model, BERNOULLI_DATA, init_radius = -0.1), - "non-negative") - expect_error(pathfinder(bernoulli_model, BERNOULLI_DATA, num_threads = 0), "positive") - expect_error(pathfinder(bernoulli_model, BERNOULLI_DATA, num_iterations = 0), - "positive") - expect_error(pathfinder(bernoulli_model, BERNOULLI_DATA, num_elbo_draws = 0), - "positive") - expect_error(pathfinder(bernoulli_model, BERNOULLI_DATA, num_multi_draws = 0), - "at least 1") - expect_error(pathfinder(bernoulli_model, BERNOULLI_DATA, max_history_size = 0), - "positive") - expect_error(pathfinder(bernoulli_model, BERNOULLI_DATA, init_alpha = 0), "positive") - expect_error(pathfinder(bernoulli_model, BERNOULLI_DATA, tol_obj = 0), "positive") - expect_error(pathfinder(bernoulli_model, BERNOULLI_DATA, tol_rel_obj = 0), "positive") - expect_error(pathfinder(bernoulli_model, BERNOULLI_DATA, tol_grad = 0), "positive") - expect_error(pathfinder(bernoulli_model, BERNOULLI_DATA, tol_rel_grad = 0), "positive") - expect_error(pathfinder(bernoulli_model, BERNOULLI_DATA, tol_param = 0), "positive") - + expect_error( + pathfinder(bernoulli_model, BERNOULLI_DATA, num_paths = 0), + "at least 1" + ) + expect_error( + pathfinder(bernoulli_model, BERNOULLI_DATA, num_draws = 0), + "at least 1" + ) + expect_error(pathfinder(bernoulli_model, BERNOULLI_DATA, id = 0), "positive") + expect_error( + pathfinder(bernoulli_model, BERNOULLI_DATA, init_radius = -0.1), + "non-negative" + ) + expect_error( + pathfinder(bernoulli_model, BERNOULLI_DATA, num_threads = 0), + "positive" + ) + expect_error( + pathfinder(bernoulli_model, BERNOULLI_DATA, num_iterations = 0), + "positive" + ) + expect_error( + pathfinder(bernoulli_model, BERNOULLI_DATA, num_elbo_draws = 0), + "positive" + ) + expect_error( + pathfinder(bernoulli_model, BERNOULLI_DATA, num_multi_draws = 0), + "at least 1" + ) + expect_error( + pathfinder(bernoulli_model, BERNOULLI_DATA, max_history_size = 0), + "positive" + ) + expect_error( + pathfinder(bernoulli_model, BERNOULLI_DATA, init_alpha = 0), + "positive" + ) + expect_error( + pathfinder(bernoulli_model, BERNOULLI_DATA, tol_obj = 0), + "positive" + ) + expect_error( + pathfinder(bernoulli_model, BERNOULLI_DATA, tol_rel_obj = 0), + "positive" + ) + expect_error( + pathfinder(bernoulli_model, BERNOULLI_DATA, tol_grad = 0), + "positive" + ) + expect_error( + pathfinder(bernoulli_model, BERNOULLI_DATA, tol_rel_grad = 0), + "positive" + ) + expect_error( + pathfinder(bernoulli_model, BERNOULLI_DATA, tol_param = 0), + "positive" + ) }) - diff --git a/clients/R/tests/testthat/test_sample.R b/clients/R/tests/testthat/test_sample.R index 323cf65..ec985c4 100644 --- a/clients/R/tests/testthat/test_sample.R +++ b/clients/R/tests/testthat/test_sample.R @@ -1,250 +1,445 @@ test_that("data arguments work", { - - out1 <- sampler(bernoulli_model, BERNOULLI_DATA, num_warmup = 100, num_samples = 100) - expect_true(mean(out1$draws$theta) > 0.2 && mean(out1$draws$theta) < 0.3) - data_file <- file.path(stan_folder, "bernoulli", "bernoulli.data.json") - out2 <- sampler(bernoulli_model, data = data_file, num_warmup = 100, num_samples = 100) - expect_true(mean(out2$draws$theta) > 0.2 && mean(out2$draws$theta) < 0.3) - + out1 <- sampler( + bernoulli_model, + BERNOULLI_DATA, + num_warmup = 100, + num_samples = 100 + ) + expect_true(mean(out1$draws$theta) > 0.2 && mean(out1$draws$theta) < 0.3) + data_file <- file.path(stan_folder, "bernoulli", "bernoulli.data.json") + out2 <- sampler( + bernoulli_model, + data = data_file, + num_warmup = 100, + num_samples = 100 + ) + expect_true(mean(out2$draws$theta) > 0.2 && mean(out2$draws$theta) < 0.3) }) test_that("save_warmup works", { - - out <- sampler(bernoulli_model, BERNOULLI_DATA, num_warmup = 12, num_samples = 34, - save_warmup = FALSE) - expect_equal(posterior::niterations(out$draws), 34) - - out <- sampler(bernoulli_model, BERNOULLI_DATA, num_warmup = 12, num_samples = 34, - save_warmup = TRUE) - expect_equal(posterior::niterations(out$draws), 12 + 34) - + out <- sampler( + bernoulli_model, + BERNOULLI_DATA, + num_warmup = 12, + num_samples = 34, + save_warmup = FALSE + ) + expect_equal(posterior::niterations(out$draws), 34) + + out <- sampler( + bernoulli_model, + BERNOULLI_DATA, + num_warmup = 12, + num_samples = 34, + save_warmup = TRUE + ) + expect_equal(posterior::niterations(out$draws), 12 + 34) }) test_that("seed works", { - - out1 <- sampler(bernoulli_model, BERNOULLI_DATA, seed = 123, num_warmup = 100, - num_samples = 100) - out2 <- sampler(bernoulli_model, BERNOULLI_DATA, seed = 123, num_warmup = 100, - num_samples = 100) - - expect_equal(out1$draws, out2$draws) - - out3 <- sampler(bernoulli_model, BERNOULLI_DATA, seed = 456, num_warmup = 100, - num_samples = 100) - - expect_error(expect_equal(out1$draws, out3$draws)) - + out1 <- sampler( + bernoulli_model, + BERNOULLI_DATA, + seed = 123, + num_warmup = 100, + num_samples = 100 + ) + out2 <- sampler( + bernoulli_model, + BERNOULLI_DATA, + seed = 123, + num_warmup = 100, + num_samples = 100 + ) + + expect_equal(out1$draws, out2$draws) + + out3 <- sampler( + bernoulli_model, + BERNOULLI_DATA, + seed = 456, + num_warmup = 100, + num_samples = 100 + ) + + expect_error(expect_equal(out1$draws, out3$draws)) }) test_that("stepsize is saved", { - out <- sampler(bernoulli_model, BERNOULLI_DATA, num_warmup = 100, num_samples = 100, - num_chains = 3, save_warmup = TRUE) - expect_true(exists("stepsize", out)) - expect_equal(length(out$stepsize), 3) - - out <- sampler(bernoulli_model, BERNOULLI_DATA, num_warmup = 100, num_samples = 100, - save_warmup = TRUE, adapt = FALSE) - expect_false(exists("stepsize", out)) + out <- sampler( + bernoulli_model, + BERNOULLI_DATA, + num_warmup = 100, + num_samples = 100, + num_chains = 3, + save_warmup = TRUE + ) + expect_true(exists("stepsize", out)) + expect_equal(length(out$stepsize), 3) + + out <- sampler( + bernoulli_model, + BERNOULLI_DATA, + num_warmup = 100, + num_samples = 100, + save_warmup = TRUE, + adapt = FALSE + ) + expect_false(exists("stepsize", out)) }) test_that("save_inv_metric works", { - - data <- "{\"N\": 5}" - - out_unit <- sampler(gaussian_model, data, num_warmup = 100, num_samples = 10, - save_inv_metric = TRUE, metric = tinystan::HMCMetric$UNIT) - expect_equal(dim(out_unit$inv_metric), c(4, 5)) - expect_equal(out_unit$inv_metric, matrix(1, ncol = 5, nrow = 4)) - - out_diag <- sampler(gaussian_model, data, num_warmup = 100, num_samples = 10, - save_inv_metric = TRUE, metric = tinystan::HMCMetric$DIAGONAL) - expect_equal(dim(out_diag$inv_metric), c(4, 5)) - expect_equal(out_diag$inv_metric, matrix(1, ncol = 5, nrow = 4)) - - out_dense <- sampler(gaussian_model, data, num_warmup = 100, num_samples = 10, - save_inv_metric = TRUE, metric = tinystan::HMCMetric$DENSE) - expect_equal(dim(out_dense$inv_metric), c(4, 5, 5)) - four_identities <- aperm(array(rep(diag(5), 4), c(5, 5, 4)), c(3, 2, 1)) - expect_equal(out_dense$inv_metric, four_identities) - - out_nometric <- sampler(gaussian_model, data, num_warmup = 10, num_samples = 10, - save_inv_metric = FALSE) - expect_false(exists("metric", out_nometric)) - + data <- "{\"N\": 5}" + + out_unit <- sampler( + gaussian_model, + data, + num_warmup = 100, + num_samples = 10, + save_inv_metric = TRUE, + metric = tinystan::HMCMetric$UNIT + ) + expect_equal(dim(out_unit$inv_metric), c(4, 5)) + expect_equal(out_unit$inv_metric, matrix(1, ncol = 5, nrow = 4)) + + out_diag <- sampler( + gaussian_model, + data, + num_warmup = 100, + num_samples = 10, + save_inv_metric = TRUE, + metric = tinystan::HMCMetric$DIAGONAL + ) + expect_equal(dim(out_diag$inv_metric), c(4, 5)) + expect_equal(out_diag$inv_metric, matrix(1, ncol = 5, nrow = 4)) + + out_dense <- sampler( + gaussian_model, + data, + num_warmup = 100, + num_samples = 10, + save_inv_metric = TRUE, + metric = tinystan::HMCMetric$DENSE + ) + expect_equal(dim(out_dense$inv_metric), c(4, 5, 5)) + four_identities <- aperm(array(rep(diag(5), 4), c(5, 5, 4)), c(3, 2, 1)) + expect_equal(out_dense$inv_metric, four_identities) + + out_nometric <- sampler( + gaussian_model, + data, + num_warmup = 10, + num_samples = 10, + save_inv_metric = FALSE + ) + expect_false(exists("metric", out_nometric)) }) test_that("init_inv_metric is used", { - data <- "{\"N\": 3}" - - for (adapt in c(TRUE, FALSE)) { - diag_metric <- matrix(1, 3, 2) - # multiply only first chain - diag_metric[, 1] <- diag_metric[, 1] * 1e+20 - out_diag <- sampler(gaussian_model, data, num_chains = 2, save_warmup = TRUE, - adapt = adapt, metric = HMCMetric$DIAGONAL, init_inv_metric = diag_metric, - save_inv_metric = TRUE, seed = 1234) - - divergent <- out_diag$draws$divergent__ - chain_one_divergences <- sum(posterior::subset_draws(divergent, chain = 1)) - expect_true(chain_one_divergences > ifelse(adapt, 12, 500)) - chain_two_divergences <- sum(posterior::subset_draws(divergent, chain = 2)) - expect_true(chain_two_divergences < 12) - expect_true(chain_two_divergences < chain_one_divergences) - if (adapt) { - expect_false(all(diag_metric == t(out_diag$inv_metric))) - } else { - expect_false(exists("inv_metric", out_diag)) - } - dense_metric <- array(0, c(3, 3, 2)) - dense_metric[, , 1] <- diag(3) * 1e+20 - dense_metric[, , 2] <- diag(3) - out_dense <- sampler(gaussian_model, data, num_chains = 2, save_warmup = TRUE, - adapt = adapt, metric = HMCMetric$DENSE, init_inv_metric = dense_metric, - save_inv_metric = TRUE, seed = 1234) - - divergent <- out_dense$draws$divergent__ - chain_one_divergences <- sum(posterior::subset_draws(divergent, chain = 1)) - expect_true(chain_one_divergences > ifelse(adapt, 12, 500)) - chain_two_divergences <- sum(posterior::subset_draws(divergent, chain = 2)) - expect_true(chain_two_divergences < 12) - expect_true(chain_two_divergences < chain_one_divergences) - if (adapt) { - expect_false(all(dense_metric == aperm(out_dense$inv_metric, c(3, 2, - 1)))) - } else { - expect_false(exists("inv_metric", out_dense)) - } + data <- "{\"N\": 3}" + + for (adapt in c(TRUE, FALSE)) { + diag_metric <- matrix(1, 3, 2) + # multiply only first chain + diag_metric[, 1] <- diag_metric[, 1] * 1e+20 + out_diag <- sampler( + gaussian_model, + data, + num_chains = 2, + save_warmup = TRUE, + adapt = adapt, + metric = HMCMetric$DIAGONAL, + init_inv_metric = diag_metric, + save_inv_metric = TRUE, + seed = 1234 + ) + + divergent <- out_diag$draws$divergent__ + chain_one_divergences <- sum(posterior::subset_draws(divergent, chain = 1)) + expect_true(chain_one_divergences > ifelse(adapt, 12, 500)) + chain_two_divergences <- sum(posterior::subset_draws(divergent, chain = 2)) + expect_true(chain_two_divergences < 12) + expect_true(chain_two_divergences < chain_one_divergences) + if (adapt) { + expect_false(all(diag_metric == t(out_diag$inv_metric))) + } else { + expect_false(exists("inv_metric", out_diag)) } + dense_metric <- array(0, c(3, 3, 2)) + dense_metric[,, 1] <- diag(3) * 1e+20 + dense_metric[,, 2] <- diag(3) + out_dense <- sampler( + gaussian_model, + data, + num_chains = 2, + save_warmup = TRUE, + adapt = adapt, + metric = HMCMetric$DENSE, + init_inv_metric = dense_metric, + save_inv_metric = TRUE, + seed = 1234 + ) + + divergent <- out_dense$draws$divergent__ + chain_one_divergences <- sum(posterior::subset_draws(divergent, chain = 1)) + expect_true(chain_one_divergences > ifelse(adapt, 12, 500)) + chain_two_divergences <- sum(posterior::subset_draws(divergent, chain = 2)) + expect_true(chain_two_divergences < 12) + expect_true(chain_two_divergences < chain_one_divergences) + if (adapt) { + expect_false(all(dense_metric == aperm(out_dense$inv_metric, c(3, 2, 1)))) + } else { + expect_false(exists("inv_metric", out_dense)) + } + } }) test_that("multiple inits work", { - - init1 <- "{\"mu\": -100}" - out1 <- sampler(multimodal_model, num_chains = 2, num_warmup = 100, num_samples = 100, - inits = init1) - expect_true(all(out1$draws$mu < 0)) - - init2 <- "{\"mu\": 100}" - out2 <- sampler(multimodal_model, num_chains = 2, num_warmup = 100, num_samples = 100, - inits = list(init1, init2)) - - expect_true(all(posterior::subset_draws(out2$draws$mu, chain = 1) < 0)) - expect_true(all(posterior::subset_draws(out2$draws$mu, chain = 2) > 0)) - - temp_file <- tempfile(fileext = ".json") - write(init1, temp_file) - out3 <- sampler(multimodal_model, num_chains = 2, num_warmup = 100, num_samples = 100, - inits = c(temp_file, init2)) - expect_true(all(posterior::subset_draws(out3$draws$mu, chain = 1) < 0)) - expect_true(all(posterior::subset_draws(out3$draws$mu, chain = 2) > 0)) - + init1 <- "{\"mu\": -100}" + out1 <- sampler( + multimodal_model, + num_chains = 2, + num_warmup = 100, + num_samples = 100, + inits = init1 + ) + expect_true(all(out1$draws$mu < 0)) + + init2 <- "{\"mu\": 100}" + out2 <- sampler( + multimodal_model, + num_chains = 2, + num_warmup = 100, + num_samples = 100, + inits = list(init1, init2) + ) + + expect_true(all(posterior::subset_draws(out2$draws$mu, chain = 1) < 0)) + expect_true(all(posterior::subset_draws(out2$draws$mu, chain = 2) > 0)) + + temp_file <- tempfile(fileext = ".json") + write(init1, temp_file) + out3 <- sampler( + multimodal_model, + num_chains = 2, + num_warmup = 100, + num_samples = 100, + inits = c(temp_file, init2) + ) + expect_true(all(posterior::subset_draws(out3$draws$mu, chain = 1) < 0)) + expect_true(all(posterior::subset_draws(out3$draws$mu, chain = 2) > 0)) }) test_that("bad data handled properly", { + data <- "{\"N\": -1}" + expect_error(sampler(bernoulli_model, data), "greater than or equal to 0") - data <- "{\"N\": -1}" - expect_error(sampler(bernoulli_model, data), "greater than or equal to 0") + data <- "{\"N\": 1, \"y\": [0, 1]}" + expect_error(sampler(bernoulli_model, data), "mismatch in dimension") - data <- "{\"N\": 1, \"y\": [0, 1]}" - expect_error(sampler(bernoulli_model, data), "mismatch in dimension") - - expect_error(sampler(bernoulli_model, "{\"bad\"}"), "Error in JSON parsing") - - expect_error(sampler(bernoulli_model, "not/real/path.json"), "Could not open data file") + expect_error(sampler(bernoulli_model, "{\"bad\"}"), "Error in JSON parsing") + expect_error( + sampler(bernoulli_model, "not/real/path.json"), + "Could not open data file" + ) }) test_that("bad inits handled properly", { - - init1 <- "{\"theta\": 2}" - expect_error(sampler(bernoulli_model, BERNOULLI_DATA, inits = init1), "Initialization failed") - - expect_error(sampler(bernoulli_model, BERNOULLI_DATA, inits = "bad/path.json"), - "Could not open data file") - - init2 <- "{\"theta\": 0.5}" - expect_error(sampler(bernoulli_model, BERNOULLI_DATA, num_chains = 2, inits = c(init2, - init1)), "Initialization failed") - - inits <- list(init2, init2) - expect_error(sampler(bernoulli_model, BERNOULLI_DATA, num_chains = 1, inits = inits), - "match the number of chains") - expect_error(sampler(bernoulli_model, BERNOULLI_DATA, num_chains = 3, inits = inits), - "match the number of chains") - + init1 <- "{\"theta\": 2}" + expect_error( + sampler(bernoulli_model, BERNOULLI_DATA, inits = init1), + "Initialization failed" + ) + + expect_error( + sampler(bernoulli_model, BERNOULLI_DATA, inits = "bad/path.json"), + "Could not open data file" + ) + + init2 <- "{\"theta\": 0.5}" + expect_error( + sampler( + bernoulli_model, + BERNOULLI_DATA, + num_chains = 2, + inits = c(init2, init1) + ), + "Initialization failed" + ) + + inits <- list(init2, init2) + expect_error( + sampler(bernoulli_model, BERNOULLI_DATA, num_chains = 1, inits = inits), + "match the number of chains" + ) + expect_error( + sampler(bernoulli_model, BERNOULLI_DATA, num_chains = 3, inits = inits), + "match the number of chains" + ) }) test_that("bad initial metric shape handled properly", { - data <- "{\"N\": 5}" - - expect_error(sampler(gaussian_model, data, metric = tinystan::HMCMetric$DENSE, - init_inv_metric = rep(1, 5)), "Invalid initial metric size") - - expect_error(sampler(gaussian_model, data, metric = tinystan::HMCMetric$DENSE, - init_inv_metric = matrix(1, 5, 4)), "Invalid initial metric size") - - expect_error(sampler(gaussian_model, data, num_chains = 4, metric = tinystan::HMCMetric$DENSE, - init_inv_metric = array(1, c(5, 5, 3))), "Invalid initial metric size") - - expect_error(sampler(gaussian_model, data, metric = tinystan::HMCMetric$DIAGONAL, - init_inv_metric = rep(1, 4)), "Invalid initial metric size") - - expect_error(sampler(gaussian_model, data, num_chains = 4, metric = tinystan::HMCMetric$DIAGONAL, - init_inv_metric = matrix(1, 5, 3)), "Invalid initial metric size") - - expect_error(sampler(gaussian_model, data, num_chains = 4, metric = tinystan::HMCMetric$DIAGONAL, - init_inv_metric = array(1, c(5, 5, 3))), "Invalid initial metric size") + data <- "{\"N\": 5}" + + expect_error( + sampler( + gaussian_model, + data, + metric = tinystan::HMCMetric$DENSE, + init_inv_metric = rep(1, 5) + ), + "Invalid initial metric size" + ) + + expect_error( + sampler( + gaussian_model, + data, + metric = tinystan::HMCMetric$DENSE, + init_inv_metric = matrix(1, 5, 4) + ), + "Invalid initial metric size" + ) + + expect_error( + sampler( + gaussian_model, + data, + num_chains = 4, + metric = tinystan::HMCMetric$DENSE, + init_inv_metric = array(1, c(5, 5, 3)) + ), + "Invalid initial metric size" + ) + + expect_error( + sampler( + gaussian_model, + data, + metric = tinystan::HMCMetric$DIAGONAL, + init_inv_metric = rep(1, 4) + ), + "Invalid initial metric size" + ) + + expect_error( + sampler( + gaussian_model, + data, + num_chains = 4, + metric = tinystan::HMCMetric$DIAGONAL, + init_inv_metric = matrix(1, 5, 3) + ), + "Invalid initial metric size" + ) + + expect_error( + sampler( + gaussian_model, + data, + num_chains = 4, + metric = tinystan::HMCMetric$DIAGONAL, + init_inv_metric = array(1, c(5, 5, 3)) + ), + "Invalid initial metric size" + ) }) test_that("bad initial metric handled properly", { - data <- "{\"N\": 3}" - - expect_error(sampler(gaussian_model, data, metric = tinystan::HMCMetric$DENSE, - init_inv_metric = matrix(1e+20, 3, 3)), "not positive definite") - - metric = array(0, c(3, 3, 2)) - metric[, , 1] <- 1e+20 - metric[, , 2] <- diag(3) - expect_error(sampler(gaussian_model, data, num_chains = 2, metric = tinystan::HMCMetric$DENSE, - init_inv_metric = metric), "not positive definite") - + data <- "{\"N\": 3}" + + expect_error( + sampler( + gaussian_model, + data, + metric = tinystan::HMCMetric$DENSE, + init_inv_metric = matrix(1e+20, 3, 3) + ), + "not positive definite" + ) + + metric = array(0, c(3, 3, 2)) + metric[,, 1] <- 1e+20 + metric[,, 2] <- diag(3) + expect_error( + sampler( + gaussian_model, + data, + num_chains = 2, + metric = tinystan::HMCMetric$DENSE, + init_inv_metric = metric + ), + "not positive definite" + ) }) test_that("bad num_warmup handled properly", { - - expect_error(sampler(bernoulli_model, BERNOULLI_DATA, num_warmup = -1), "non-negative") - expect_error(sampler(bernoulli_model, BERNOULLI_DATA, save_warmup = TRUE, num_warmup = -1), - "non-negative") - + expect_error( + sampler(bernoulli_model, BERNOULLI_DATA, num_warmup = -1), + "non-negative" + ) + expect_error( + sampler( + bernoulli_model, + BERNOULLI_DATA, + save_warmup = TRUE, + num_warmup = -1 + ), + "non-negative" + ) }) test_that("model with no params fails", { - - expect_no_error(sampler(empty_model)) - + expect_no_error(sampler(empty_model)) }) test_that("bad args raise errors", { - - expect_error(sampler(bernoulli_model, BERNOULLI_DATA, num_chains = 0), "at least 1") - expect_error(sampler(bernoulli_model, BERNOULLI_DATA, num_samples = 0), "at least 1") - expect_error(sampler(bernoulli_model, BERNOULLI_DATA, id = 0), "positive") - expect_error(sampler(bernoulli_model, BERNOULLI_DATA, init_radius = -0.1), "non-negative") - - expect_error(sampler(bernoulli_model, BERNOULLI_DATA, delta = -0.1), "between 0 and 1") - expect_error(sampler(bernoulli_model, BERNOULLI_DATA, delta = 1.1), "between 0 and 1") - - expect_error(sampler(bernoulli_model, BERNOULLI_DATA, gamma = 0), "positive") - expect_error(sampler(bernoulli_model, BERNOULLI_DATA, kappa = 0), "positive") - expect_error(sampler(bernoulli_model, BERNOULLI_DATA, t0 = 0), "positive") - - expect_error(sampler(bernoulli_model, BERNOULLI_DATA, stepsize = 0), "positive") - expect_error(sampler(bernoulli_model, BERNOULLI_DATA, stepsize_jitter = -0.1), - "between 0 and 1") - expect_error(sampler(bernoulli_model, BERNOULLI_DATA, stepsize_jitter = 1.1), - "between 0 and 1") - - expect_error(sampler(bernoulli_model, BERNOULLI_DATA, max_depth = 0), "positive") - expect_error(sampler(bernoulli_model, BERNOULLI_DATA, num_threads = 0), "positive") - + expect_error( + sampler(bernoulli_model, BERNOULLI_DATA, num_chains = 0), + "at least 1" + ) + expect_error( + sampler(bernoulli_model, BERNOULLI_DATA, num_samples = 0), + "at least 1" + ) + expect_error(sampler(bernoulli_model, BERNOULLI_DATA, id = 0), "positive") + expect_error( + sampler(bernoulli_model, BERNOULLI_DATA, init_radius = -0.1), + "non-negative" + ) + + expect_error( + sampler(bernoulli_model, BERNOULLI_DATA, delta = -0.1), + "between 0 and 1" + ) + expect_error( + sampler(bernoulli_model, BERNOULLI_DATA, delta = 1.1), + "between 0 and 1" + ) + + expect_error(sampler(bernoulli_model, BERNOULLI_DATA, gamma = 0), "positive") + expect_error(sampler(bernoulli_model, BERNOULLI_DATA, kappa = 0), "positive") + expect_error(sampler(bernoulli_model, BERNOULLI_DATA, t0 = 0), "positive") + + expect_error( + sampler(bernoulli_model, BERNOULLI_DATA, stepsize = 0), + "positive" + ) + expect_error( + sampler(bernoulli_model, BERNOULLI_DATA, stepsize_jitter = -0.1), + "between 0 and 1" + ) + expect_error( + sampler(bernoulli_model, BERNOULLI_DATA, stepsize_jitter = 1.1), + "between 0 and 1" + ) + + expect_error( + sampler(bernoulli_model, BERNOULLI_DATA, max_depth = 0), + "positive" + ) + expect_error( + sampler(bernoulli_model, BERNOULLI_DATA, num_threads = 0), + "positive" + ) }) - diff --git a/clients/julia/src/model.jl b/clients/julia/src/model.jl index abbb159..a9a8259 100644 --- a/clients/julia/src/model.jl +++ b/clients/julia/src/model.jl @@ -545,7 +545,10 @@ function laplace_sample( end with_model(model, data, seed) do model_ptr - required_params = @ccall $(dlsym(model.lib, :tinystan_model_num_constrained_params_for_unconstraining))( + required_params = @ccall $(dlsym( + model.lib, + :tinystan_model_num_constrained_params_for_unconstraining, + ))( model_ptr::Ptr{Cvoid}, )::Cint diff --git a/clients/julia/test/test_laplace.jl b/clients/julia/test/test_laplace.jl index e4c7238..88726eb 100644 --- a/clients/julia/test/test_laplace.jl +++ b/clients/julia/test/test_laplace.jl @@ -131,11 +131,7 @@ ) mode2 = [1.0] - @test_throws "incorrect length" laplace_sample( - gaussian_model, - mode2, - "{\"N\": 4 }", - ) + @test_throws "incorrect length" laplace_sample(gaussian_model, mode2, "{\"N\": 4 }") end @testset "Bad mode json" begin diff --git a/clients/julia/test/test_pathfinder.jl b/clients/julia/test/test_pathfinder.jl index 6d8e631..7ad5b88 100644 --- a/clients/julia/test/test_pathfinder.jl +++ b/clients/julia/test/test_pathfinder.jl @@ -79,10 +79,10 @@ @testset "Seed" begin out1 = pathfinder(bernoulli_model, BERNOULLI_DATA; seed = UInt32(123)) out2 = pathfinder(bernoulli_model, BERNOULLI_DATA; seed = UInt32(123)) - @test sortslices(out1.draws, dims=1) == sortslices(out2.draws, dims=1) + @test sortslices(out1.draws, dims = 1) == sortslices(out2.draws, dims = 1) out3 = pathfinder(bernoulli_model, BERNOULLI_DATA; seed = UInt32(456)) - @test sortslices(out1.draws, dims=1) != sortslices(out3.draws, dims=1) + @test sortslices(out1.draws, dims = 1) != sortslices(out3.draws, dims = 1) end @testset "Inits" begin diff --git a/clients/julia/test/test_sample.jl b/clients/julia/test/test_sample.jl index 46b7248..0f7b42b 100644 --- a/clients/julia/test/test_sample.jl +++ b/clients/julia/test/test_sample.jl @@ -220,8 +220,8 @@ num_chains = 2, inits = [init1, init2], ) - @test all(output2.draws[1, :, output2.names.=="mu", 1] .< 0) - @test all(output2.draws[2, :, output2.names.=="mu", 1] .> 0) + @test all(output2.draws[1, :, output2.names .== "mu", 1] .< 0) + @test all(output2.draws[2, :, output2.names .== "mu", 1] .> 0) init3 = tempname() * ".json" open(init3, "w") do io @@ -234,8 +234,8 @@ num_chains = 2, inits = [init3, init2], ) - @test all(output3.draws[1, :, output3.names.=="mu", 1] .< 0) - @test all(output3.draws[2, :, output3.names.=="mu", 1] .> 0) + @test all(output3.draws[1, :, output3.names .== "mu", 1] .< 0) + @test all(output3.draws[2, :, output3.names .== "mu", 1] .> 0) end diff --git a/clients/python/tinystan/model.py b/clients/python/tinystan/model.py index 594e7af..a197c88 100644 --- a/clients/python/tinystan/model.py +++ b/clients/python/tinystan/model.py @@ -1011,10 +1011,10 @@ def laplace_sample( with self._get_model(data, seed) as model: req_params = self._num_req_constrained_params(model) if mode_array is not None and len(mode_array) < req_params: - raise ValueError( + raise ValueError( "Mode array has incorrect length. " f"Expected at least {req_params} but got {len(mode_array)}" - ) + ) param_names = LAPLACE_VARIABLES + self._get_parameter_names(model) num_params = len(param_names) diff --git a/src/buffer.hpp b/src/buffer.hpp index 78b0e35..3558d90 100644 --- a/src/buffer.hpp +++ b/src/buffer.hpp @@ -41,8 +41,8 @@ namespace io { */ class buffer_writer : public stan::callbacks::writer { public: - buffer_writer(double *buf, size_t max) : buf(buf), pos(0), size(max) {}; - virtual ~buffer_writer() {}; + buffer_writer(double *buf, size_t max) : buf(buf), pos(0), size(max){}; + virtual ~buffer_writer(){}; /** * Primary method used by the Stan algorithms @@ -117,7 +117,7 @@ class buffer_writer : public stan::callbacks::writer { class filtered_writer : public stan::callbacks::structured_writer { public: filtered_writer() : filters{} {}; - virtual ~filtered_writer() {}; + virtual ~filtered_writer(){}; void add_key(const std::string &key_in, double *buf) { if (buf != nullptr) { @@ -172,8 +172,8 @@ class inv_metric_buffer_reader : public stan::io::empty_var_context { public: inv_metric_buffer_reader(const double *buf, size_t size, TinyStanMetric metric_choice) - : buf(buf), size(size), dense(metric_choice == TinyStanMetric::dense) {}; - virtual ~inv_metric_buffer_reader() {}; + : buf(buf), size(size), dense(metric_choice == TinyStanMetric::dense){}; + virtual ~inv_metric_buffer_reader(){}; bool contains_r(const std::string &name) const override { return name == "inv_metric"; diff --git a/src/errors.hpp b/src/errors.hpp index d857de6..ab53025 100644 --- a/src/errors.hpp +++ b/src/errors.hpp @@ -78,8 +78,8 @@ inline auto catch_exceptions(TinyStanError **err, F f) { class error_logger : public stan::callbacks::logger { public: error_logger(const TinyStanModel &model, bool print_non_errors) - : model(model), print(print_non_errors) {}; - virtual ~error_logger() {}; + : model(model), print(print_non_errors){}; + virtual ~error_logger(){}; void info(const std::string &s) override { if (print && !s.empty()) {