From 2a332c724523ef524522c6840c801a5634c96c79 Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Mon, 27 Mar 2023 15:16:03 -0400 Subject: [PATCH 01/20] Start C++ changes for moving RNGs outside model object --- src/bridgestan.cpp | 65 ++++++++++++++++++++++++++----------------- src/bridgestan.h | 68 ++++++++++++++++++++++++++------------------- src/bridgestanR.cpp | 33 +++++++++++----------- src/bridgestanR.h | 36 ++++++++++++------------ src/model_rng.cpp | 50 ++++++++++++++++----------------- src/model_rng.hpp | 20 ++++++++----- 6 files changed, 149 insertions(+), 123 deletions(-) diff --git a/src/bridgestan.cpp b/src/bridgestan.cpp index 81d80068..862d810b 100644 --- a/src/bridgestan.cpp +++ b/src/bridgestan.cpp @@ -9,10 +9,10 @@ int bs_patch_version = BRIDGESTAN_PATCH; #include -bs_model_rng* bs_construct(const char* data_file, unsigned int seed, - unsigned int chain_id, char** error_msg) { +bs_model* bs_construct(const char* data_file, unsigned int seed, + unsigned int chain_id, char** error_msg) { try { - return new bs_model_rng(data_file, seed, chain_id); + return new bs_model(data_file, seed); } catch (const std::exception& e) { if (error_msg) { std::stringstream error; @@ -34,34 +34,34 @@ bs_model_rng* bs_construct(const char* data_file, unsigned int seed, return nullptr; } -void bs_destruct(bs_model_rng* mr) { delete (mr); } +void bs_destruct(bs_model* m) { delete (m); } void bs_free_error_msg(char* error_msg) { free(error_msg); } -const char* bs_name(const bs_model_rng* mr) { return mr->name(); } +const char* bs_name(const bs_model* m) { return m->name(); } -const char* bs_model_info(const bs_model_rng* mr) { return mr->model_info(); } +const char* bs_model_info(const bs_model* m) { return m->model_info(); } -const char* bs_param_names(const bs_model_rng* mr, bool include_tp, +const char* bs_param_names(const bs_model* m, bool include_tp, bool include_gq) { - return mr->param_names(include_tp, include_gq); + return m->param_names(include_tp, include_gq); } -const char* bs_param_unc_names(const bs_model_rng* mr) { - return mr->param_unc_names(); +const char* bs_param_unc_names(const bs_model* m) { + return m->param_unc_names(); } -int bs_param_num(const bs_model_rng* mr, bool include_tp, bool include_gq) { - return mr->param_num(include_tp, include_gq); +int bs_param_num(const bs_model* m, bool include_tp, bool include_gq) { + return m->param_num(include_tp, include_gq); } -int bs_param_unc_num(const bs_model_rng* mr) { return mr->param_unc_num(); } +int bs_param_unc_num(const bs_model* m) { return m->param_unc_num(); } -int bs_param_constrain(bs_model_rng* mr, bool include_tp, bool include_gq, - const double* theta_unc, double* theta, +int bs_param_constrain(bs_model* m, bool include_tp, bool include_gq, + const double* theta_unc, double* theta, boost::ecuyer1988& rng, char** error_msg) { try { - mr->param_constrain(include_tp, include_gq, theta_unc, theta); + m->param_constrain(include_tp, include_gq, theta_unc, theta, rng); return 0; } catch (const std::exception& e) { if (error_msg) { @@ -80,10 +80,10 @@ int bs_param_constrain(bs_model_rng* mr, bool include_tp, bool include_gq, return 1; } -int bs_param_unconstrain(const bs_model_rng* mr, const double* theta, +int bs_param_unconstrain(const bs_model* m, const double* theta, double* theta_unc, char** error_msg) { try { - mr->param_unconstrain(theta, theta_unc); + m->param_unconstrain(theta, theta_unc); return 0; } catch (const std::exception& e) { if (error_msg) { @@ -102,10 +102,10 @@ int bs_param_unconstrain(const bs_model_rng* mr, const double* theta, return -1; } -int bs_param_unconstrain_json(const bs_model_rng* mr, const char* json, +int bs_param_unconstrain_json(const bs_model* m, const char* json, double* theta_unc, char** error_msg) { try { - mr->param_unconstrain_json(json, theta_unc); + m->param_unconstrain_json(json, theta_unc); return 0; } catch (const std::exception& e) { if (error_msg) { @@ -125,10 +125,10 @@ int bs_param_unconstrain_json(const bs_model_rng* mr, const char* json, return -1; } -int bs_log_density(const bs_model_rng* mr, bool propto, bool jacobian, +int bs_log_density(const bs_model* m, bool propto, bool jacobian, const double* theta_unc, double* val, char** error_msg) { try { - mr->log_density(propto, jacobian, theta_unc, val); + m->log_density(propto, jacobian, theta_unc, val); return 0; } catch (const std::exception& e) { if (error_msg) { @@ -146,11 +146,11 @@ int bs_log_density(const bs_model_rng* mr, bool propto, bool jacobian, return -1; } -int bs_log_density_gradient(const bs_model_rng* mr, bool propto, bool jacobian, +int bs_log_density_gradient(const bs_model* m, bool propto, bool jacobian, const double* theta_unc, double* val, double* grad, char** error_msg) { try { - mr->log_density_gradient(propto, jacobian, theta_unc, val, grad); + m->log_density_gradient(propto, jacobian, theta_unc, val, grad); return 0; } catch (const std::exception& e) { if (error_msg) { @@ -170,11 +170,11 @@ int bs_log_density_gradient(const bs_model_rng* mr, bool propto, bool jacobian, return -1; } -int bs_log_density_hessian(const bs_model_rng* mr, bool propto, bool jacobian, +int bs_log_density_hessian(const bs_model* m, bool propto, bool jacobian, const double* theta_unc, double* val, double* grad, double* hessian, char** error_msg) { try { - mr->log_density_hessian(propto, jacobian, theta_unc, val, grad, hessian); + m->log_density_hessian(propto, jacobian, theta_unc, val, grad, hessian); return 0; } catch (const std::exception& e) { if (error_msg) { @@ -193,3 +193,16 @@ int bs_log_density_hessian(const bs_model_rng* mr, bool propto, bool jacobian, } return -1; } + +bs_rng* bs_rng_construct(unsigned int seed) { + try { + return new bs_rng(seed); + } catch (const std::exception& e) { + std::cerr << "bs_rng_construct(" << seed + << ") failed with exception: " << e.what() << std::endl; + } catch (...) { + std::cerr << "bs_rng_construct(" << seed + << ") failed with unknown exception" << std::endl; + } + return nullptr; +} diff --git a/src/bridgestan.h b/src/bridgestan.h index d2d05d9d..ac5faebb 100644 --- a/src/bridgestan.h +++ b/src/bridgestan.h @@ -5,7 +5,8 @@ #include "model_rng.hpp" extern "C" { #else -typedef struct bs_model_rng bs_model_rng; +typedef struct bs_model bs_model; +typedef struct bs_rng bs_rng; typedef int bool; #endif @@ -29,15 +30,15 @@ extern int bs_patch_version; * @return pointer to constructed model or `nullptr` if construction * fails */ -bs_model_rng* bs_construct(const char* data_file, unsigned int seed, - unsigned int chain_id, char** error_msg); +bs_model* bs_construct(const char* data_file, unsigned int seed, + unsigned int chain_id, char** error_msg); /** * Destroy the model. * - * @param[in] mr pointer to model and RNG structure + * @param[in] m pointer to model and RNG structure */ -void bs_destruct(bs_model_rng* mr); +void bs_destruct(bs_model* m); /** * Free the error messages created by other methods. @@ -52,10 +53,10 @@ void bs_free_error_msg(char* error_msg); * The returned string should not be modified; it is freed when the * model and RNG wrapper is destroyed. * - * @param[in] mr pointer to model and RNG structure + * @param[in] m pointer to model and RNG structure * @return name of model */ -const char* bs_name(const bs_model_rng* mr); +const char* bs_name(const bs_model* m); /** * Return information about the compiled model as a C-style string. @@ -63,11 +64,11 @@ const char* bs_name(const bs_model_rng* mr); * The returned string should not be modified; it is freed when the * model and RNG wrapper is destroyed. * - * @param[in] mr pointer to model and RNG structure + * @param[in] m pointer to model and RNG structure * @return Information about the model including Stan version, Stan defines, and * compiler flags. */ -const char* bs_model_info(const bs_model_rng* mr); +const char* bs_model_info(const bs_model* m); /** * Return a comma-separated sequence of indexed parameter names, @@ -83,12 +84,12 @@ const char* bs_model_info(const bs_model_rng* mr); * The returned string should not be modified; it is freed when the * model and RNG wrapper is destroyed. * - * @param[in] mr pointer to model and RNG structure + * @param[in] m pointer to model and RNG structure * @param[in] include_tp `true` to include transformed parameters * @param[in] include_gq `true` to include generated quantities * @return CSV-separated, indexed, parameter names */ -const char* bs_param_names(const bs_model_rng* mr, bool include_tp, +const char* bs_param_names(const bs_model* m, bool include_tp, bool include_gq); /** @@ -105,22 +106,22 @@ const char* bs_param_names(const bs_model_rng* mr, bool include_tp, * The returned string should not be modified; it is freed when the * model and RNG wrapper is destroyed. * - * @param[in] mr pointer to model and RNG structure + * @param[in] m pointer to model and RNG structure * @return CSV-separated, indexed, unconstrained parameter names */ -const char* bs_param_unc_names(const bs_model_rng* mr); +const char* bs_param_unc_names(const bs_model* m); /** * Return the number of scalar parameters, optionally including the * number of transformed parameters and/or generated quantities. * For example, a 2 x 3 matrix counts as 6 scalar parameters. * - * @param[in] mr pointer to model and RNG structure + * @param[in] m pointer to model and RNG structure * @param[in] include_tp `true` to include transformed parameters * @param[in] include_gq `true` to include generated quantities * @return number of parameters */ -int bs_param_num(const bs_model_rng* mr, bool include_tp, bool include_gq); +int bs_param_num(const bs_model* m, bool include_tp, bool include_gq); /** * Return the number of unconstrained parameters. The number of @@ -128,10 +129,10 @@ int bs_param_num(const bs_model_rng* mr, bool include_tp, bool include_gq); * parameters if the unconstrained space has fewer dimensions than * the constrained (e.g., for simplexes or correlation matrices). * - * @param[in] mr pointer to model and RNG structure + * @param[in] m pointer to model and RNG structure * @return number of unconstrained parameters */ -int bs_param_unc_num(const bs_model_rng* mr); +int bs_param_unc_num(const bs_model* m); /** * Set the sequence of constrained parameters based on the specified @@ -141,7 +142,7 @@ int bs_param_unc_num(const bs_model_rng* mr); * in the Stan program, with multivariate parameters given in * last-index-major order. * - * @param[in] mr pointer to model and RNG structure + * @param[in] m pointer to model and RNG structure * @param[in] include_tp `true` to include transformed parameters * @param[in] include_gq `true` to include generated quantities * @param[in] theta_unc sequence of unconstrained parameters @@ -151,7 +152,7 @@ int bs_param_unc_num(const bs_model_rng* mr); * @return code 0 if successful and code -1 if there is an exception * in the underlying Stan code */ -int bs_param_constrain(bs_model_rng* mr, bool include_tp, bool include_gq, +int bs_param_constrain(bs_model* m, bool include_tp, bool include_gq, const double* theta_unc, double* theta, char** error_msg); @@ -162,7 +163,7 @@ int bs_param_constrain(bs_model_rng* mr, bool include_tp, bool include_gq, * in the Stan program, with multivariate parameters given in * last-index-major order. * - * @param[in] mr pointer to model and RNG structure + * @param[in] m pointer to model and RNG structure * @param[in] theta sequence of constrained parameters * @param[out] theta_unc sequence of unconstrained parameters * @param[out] error_msg a pointer to a string that will be allocated if there @@ -170,7 +171,7 @@ int bs_param_constrain(bs_model_rng* mr, bool include_tp, bool include_gq, * @return code 0 if successful and code -1 if there is an exception * in the underlying Stan code */ -int bs_param_unconstrain(const bs_model_rng* mr, const double* theta, +int bs_param_unconstrain(const bs_model* m, const double* theta, double* theta_unc, char** error_msg); /** @@ -181,7 +182,7 @@ int bs_param_unconstrain(const bs_model_rng* mr, const double* theta, * in last-index-major order. The JSON schema assumed is fully * defined in the *CmdStan Reference Manual*. * - * @param[in] mr pointer to model and RNG structure + * @param[in] m pointer to model and RNG structure * @param[in] json JSON-encoded constrained parameters * @param[out] theta_unc sequence of unconstrained parameters * @param[out] error_msg a pointer to a string that will be allocated if there @@ -189,7 +190,7 @@ int bs_param_unconstrain(const bs_model_rng* mr, const double* theta, * @return code 0 if successful and code -1 if there is an exception * in the underlying Stan code */ -int bs_param_unconstrain_json(const bs_model_rng* mr, const char* json, +int bs_param_unconstrain_json(const bs_model* m, const char* json, double* theta_unc, char** error_msg); /** @@ -199,7 +200,7 @@ int bs_param_unconstrain_json(const bs_model_rng* mr, const char* json, * and return a return code of 0 for success and -1 if there is an * exception executing the Stan program. * - * @param[in] mr pointer to model and RNG structure + * @param[in] m pointer to model and RNG structure * @param[in] propto `true` to discard constant terms * @param[in] jacobian `true` to include change-of-variables terms * @param[in] theta_unc unconstrained parameters @@ -209,7 +210,7 @@ int bs_param_unconstrain_json(const bs_model_rng* mr, const char* json, * @return code 0 if successful and code -1 if there is an exception * in the underlying Stan code */ -int bs_log_density(const bs_model_rng* mr, bool propto, bool jacobian, +int bs_log_density(const bs_model* m, bool propto, bool jacobian, const double* theta_unc, double* lp, char** error_msg); /** @@ -222,7 +223,7 @@ int bs_log_density(const bs_model_rng* mr, bool propto, bool jacobian, * * The gradients are computed using automatic differentiation. * - * @param[in] mr pointer to model and RNG structure + * @param[in] m pointer to model and RNG structure * @param[in] propto `true` to discard constant terms * @param[in] jacobian `true` to include change-of-variables terms * @param[in] theta_unc unconstrained parameters @@ -233,7 +234,7 @@ int bs_log_density(const bs_model_rng* mr, bool propto, bool jacobian, * @return code 0 if successful and code -1 if there is an exception * in the underlying Stan code */ -int bs_log_density_gradient(const bs_model_rng* mr, bool propto, bool jacobian, +int bs_log_density_gradient(const bs_model* m, bool propto, bool jacobian, const double* theta_unc, double* val, double* grad, char** error_msg); @@ -249,7 +250,7 @@ int bs_log_density_gradient(const bs_model_rng* mr, bool propto, bool jacobian, * The gradients are computed using automatic differentiation. the * Hessians are * - * @param[in] mr pointer to model and RNG structure + * @param[in] m pointer to model and RNG structure * @param[in] propto `true` to discard constant terms * @param[in] jacobian `true` to include change-of-variables terms * @param[in] theta_unc unconstrained parameters @@ -261,10 +262,19 @@ int bs_log_density_gradient(const bs_model_rng* mr, bool propto, bool jacobian, * @return code 0 if successful and code -1 if there is an exception * in the underlying Stan code */ -int bs_log_density_hessian(const bs_model_rng* mr, bool propto, bool jacobian, +int bs_log_density_hessian(const bs_model* m, bool propto, bool jacobian, const double* theta_unc, double* val, double* grad, double* hessian, char** error_msg); +/** + * Construct an RNG object to be used for param_constrain + * + * @param[in] seed seed for the RNG + */ +bs_rng* bs_construct_rng(unsigned int seed); + +// TODO destruct as well + #ifdef __cplusplus } #endif diff --git a/src/bridgestanR.cpp b/src/bridgestanR.cpp index b702c7df..94ecd833 100644 --- a/src/bridgestanR.cpp +++ b/src/bridgestanR.cpp @@ -1,7 +1,7 @@ #include "bridgestanR.h" #include "bridgestan.h" -void bs_construct_R(char** data, int* rng, int* chain, bs_model_rng** ptr_out, +void bs_construct_R(char** data, int* rng, int* chain, bs_model** ptr_out, char** err_msg, void** err_ptr) { *ptr_out = bs_construct(*data, *rng, *chain, err_msg); *err_ptr = static_cast(*err_msg); @@ -14,55 +14,54 @@ void bs_version_R(int* major, int* minor, int* patch) { *minor = bs_minor_version; *patch = bs_patch_version; } -void bs_destruct_R(bs_model_rng** model) { bs_destruct(*model); } -void bs_name_R(bs_model_rng** model, char const** name_out) { +void bs_destruct_R(bs_model** model) { bs_destruct(*model); } +void bs_name_R(bs_model** model, char const** name_out) { *name_out = bs_name(*model); } -void bs_model_info_R(bs_model_rng** model, char const** info_out) { +void bs_model_info_R(bs_model** model, char const** info_out) { *info_out = bs_model_info(*model); } -void bs_param_names_R(bs_model_rng** model, int* include_tp, int* include_gq, +void bs_param_names_R(bs_model** model, int* include_tp, int* include_gq, char const** names_out) { *names_out = bs_param_names(*model, *include_tp, *include_gq); } -void bs_param_unc_names_R(bs_model_rng** model, char const** names_out) { +void bs_param_unc_names_R(bs_model** model, char const** names_out) { *names_out = bs_param_unc_names(*model); } -void bs_param_num_R(bs_model_rng** model, int* include_tp, int* include_gq, +void bs_param_num_R(bs_model** model, int* include_tp, int* include_gq, int* num_out) { *num_out = bs_param_num(*model, *include_tp, *include_gq); } -void bs_param_unc_num_R(bs_model_rng** model, int* num_out) { +void bs_param_unc_num_R(bs_model** model, int* num_out) { *num_out = bs_param_unc_num(*model); } -void bs_param_constrain_R(bs_model_rng** model, int* include_tp, - int* include_gq, const double* theta_unc, - double* theta, int* return_code, char** err_msg, - void** err_ptr) { +void bs_param_constrain_R(bs_model** model, int* include_tp, int* include_gq, + const double* theta_unc, double* theta, + int* return_code, char** err_msg, void** err_ptr) { *return_code = bs_param_constrain(*model, *include_tp, *include_gq, theta_unc, theta, err_msg); *err_ptr = static_cast(*err_msg); } -void bs_param_unconstrain_R(bs_model_rng** model, const double* theta, +void bs_param_unconstrain_R(bs_model** model, const double* theta, double* theta_unc, int* return_code, char** err_msg, void** err_ptr) { *return_code = bs_param_unconstrain(*model, theta, theta_unc, err_msg); *err_ptr = static_cast(*err_msg); } -void bs_param_unconstrain_json_R(bs_model_rng** model, char const** json, +void bs_param_unconstrain_json_R(bs_model** model, char const** json, double* theta_unc, int* return_code, char** err_msg, void** err_ptr) { *return_code = bs_param_unconstrain_json(*model, *json, theta_unc, err_msg); *err_ptr = static_cast(*err_msg); } -void bs_log_density_R(bs_model_rng** model, int* propto, int* jacobian, +void bs_log_density_R(bs_model** model, int* propto, int* jacobian, const double* theta_unc, double* val, int* return_code, char** err_msg, void** err_ptr) { *return_code = bs_log_density(*model, *propto, *jacobian, theta_unc, val, err_msg); *err_ptr = static_cast(*err_msg); } -void bs_log_density_gradient_R(bs_model_rng** model, int* propto, int* jacobian, +void bs_log_density_gradient_R(bs_model** model, int* propto, int* jacobian, const double* theta_unc, double* val, double* grad, int* return_code, char** err_msg, void** err_ptr) { @@ -70,7 +69,7 @@ void bs_log_density_gradient_R(bs_model_rng** model, int* propto, int* jacobian, val, grad, err_msg); *err_ptr = static_cast(*err_msg); } -void bs_log_density_hessian_R(bs_model_rng** model, int* propto, int* jacobian, +void bs_log_density_hessian_R(bs_model** model, int* propto, int* jacobian, const double* theta_unc, double* val, double* grad, double* hess, int* return_code, char** err_msg, void** err_ptr) { diff --git a/src/bridgestanR.h b/src/bridgestanR.h index d3a54836..9eb6c16f 100644 --- a/src/bridgestanR.h +++ b/src/bridgestanR.h @@ -5,18 +5,19 @@ #include "model_rng.hpp" extern "C" { #else -typedef struct bs_model_rng bs_model_rng; +typedef struct bs_model bs_model; +typedef struct bs_rng bs_rng; typedef int bool; #endif // Shim to convert to R interface requirement of void with pointer args // All calls directly delegated to versions without _R suffix -void bs_construct_R(char** data, int* rng, int* chain, bs_model_rng** ptr_out, +void bs_construct_R(char** data, int* rng, int* chain, bs_model** ptr_out, char** err_msg, void** err_ptr); void bs_version_R(int* major, int* minor, int* patch); -void bs_destruct_R(bs_model_rng** model); +void bs_destruct_R(bs_model** model); /** * Free error message allocated in C++. Because R performs copies @@ -24,43 +25,42 @@ void bs_destruct_R(bs_model_rng** model); */ void bs_free_error_msg_R(void** err_msg); -void bs_name_R(bs_model_rng** model, char const** name_out); +void bs_name_R(bs_model** model, char const** name_out); -void bs_model_info_R(bs_model_rng** model, char const** info_out); +void bs_model_info_R(bs_model** model, char const** info_out); -void bs_param_names_R(bs_model_rng** model, int* include_tp, int* include_gq, +void bs_param_names_R(bs_model** model, int* include_tp, int* include_gq, char const** name_out); -void bs_param_unc_names_R(bs_model_rng** model, char const** name_out); +void bs_param_unc_names_R(bs_model** model, char const** name_out); -void bs_param_num_R(bs_model_rng** model, int* include_tp, int* include_gq, +void bs_param_num_R(bs_model** model, int* include_tp, int* include_gq, int* num_out); -void bs_param_unc_num_R(bs_model_rng** model, int* num_out); +void bs_param_unc_num_R(bs_model** model, int* num_out); -void bs_param_constrain_R(bs_model_rng** model, int* include_tp, - int* include_gq, const double* theta_unc, - double* theta, int* return_code, char** err_msg, - void** err_ptr); +void bs_param_constrain_R(bs_model** model, int* include_tp, int* include_gq, + const double* theta_unc, double* theta, + int* return_code, char** err_msg, void** err_ptr); -void bs_param_unconstrain_R(bs_model_rng** model, const double* theta, +void bs_param_unconstrain_R(bs_model** model, const double* theta, double* theta_unc, int* return_code, char** err_msg, void** err_ptr); -void bs_param_unconstrain_json_R(bs_model_rng** model, char const** json, +void bs_param_unconstrain_json_R(bs_model** model, char const** json, double* theta_unc, int* return_code, char** err_msg, void** err_ptr); -void bs_log_density_R(bs_model_rng** model, int* propto, int* jacobian, +void bs_log_density_R(bs_model** model, int* propto, int* jacobian, const double* theta_unc, double* val, int* return_code, char** err_msg, void** err_ptr); -void bs_log_density_gradient_R(bs_model_rng** model, int* propto, int* jacobian, +void bs_log_density_gradient_R(bs_model** model, int* propto, int* jacobian, const double* theta_unc, double* val, double* grad, int* return_code, char** err_msg, void** err_ptr); -void bs_log_density_hessian_R(bs_model_rng** model, int* propto, int* jacobian, +void bs_log_density_hessian_R(bs_model** model, int* propto, int* jacobian, const double* theta_unc, double* val, double* grad, double* hess, int* return_code, char** err_msg, void** err_ptr); diff --git a/src/model_rng.cpp b/src/model_rng.cpp index baf3794d..ce0428da 100644 --- a/src/model_rng.cpp +++ b/src/model_rng.cpp @@ -58,8 +58,7 @@ char* to_csv(const std::vector& names) { return strdup(s_c); } -bs_model_rng::bs_model_rng(const char* data_file, unsigned int seed, - unsigned int chain_id) { +bs_model::bs_model(const char* data_file, unsigned int seed) { std::string data(data_file); if (data.empty()) { auto data_context = stan::io::empty_var_context(); @@ -78,7 +77,6 @@ bs_model_rng::bs_model_rng(const char* data_file, unsigned int seed, model_ = &new_model(data_context, seed, &std::cout); } } - rng_ = stan::services::util::create_rng(seed, chain_id); std::string model_name = model_->model_name(); const char* model_name_c = model_name.c_str(); @@ -156,7 +154,7 @@ bs_model_rng::bs_model_rng(const char* data_file, unsigned int seed, param_tp_gq_num_ = names.size(); } -bs_model_rng::~bs_model_rng() noexcept { +bs_model::~bs_model() noexcept { delete (model_); free(name_); free(model_info_); @@ -167,11 +165,11 @@ bs_model_rng::~bs_model_rng() noexcept { free(param_tp_gq_names_); } -const char* bs_model_rng::name() const { return name_; } +const char* bs_model::name() const { return name_; } -const char* bs_model_rng::model_info() const { return model_info_; } +const char* bs_model::model_info() const { return model_info_; } -const char* bs_model_rng::param_names(bool include_tp, bool include_gq) const { +const char* bs_model::param_names(bool include_tp, bool include_gq) const { if (include_tp && include_gq) return param_tp_gq_names_; if (include_tp) @@ -181,11 +179,11 @@ const char* bs_model_rng::param_names(bool include_tp, bool include_gq) const { return param_names_; } -const char* bs_model_rng::param_unc_names() const { return param_unc_names_; } +const char* bs_model::param_unc_names() const { return param_unc_names_; } -int bs_model_rng::param_unc_num() const { return param_unc_num_; } +int bs_model::param_unc_num() const { return param_unc_num_; } -int bs_model_rng::param_num(bool include_tp, bool include_gq) const { +int bs_model::param_num(bool include_tp, bool include_gq) const { if (include_tp && include_gq) return param_tp_gq_num_; if (include_tp) @@ -195,8 +193,7 @@ int bs_model_rng::param_num(bool include_tp, bool include_gq) const { return param_num_; } -void bs_model_rng::param_unconstrain(const double* theta, - double* theta_unc) const { +void bs_model::param_unconstrain(const double* theta, double* theta_unc) const { using std::set; using std::string; using std::vector; @@ -229,8 +226,8 @@ void bs_model_rng::param_unconstrain(const double* theta, Eigen::VectorXd::Map(theta_unc, unc_params.size()) = unc_params; } -void bs_model_rng::param_unconstrain_json(const char* json, - double* theta_unc) const { +void bs_model::param_unconstrain_json(const char* json, + double* theta_unc) const { std::stringstream in(json); stan::json::json_data inits_context(in); Eigen::VectorXd params_unc; @@ -238,17 +235,18 @@ void bs_model_rng::param_unconstrain_json(const char* json, Eigen::VectorXd::Map(theta_unc, params_unc.size()) = params_unc; } -void bs_model_rng::param_constrain(bool include_tp, bool include_gq, - const double* theta_unc, double* theta) { +void bs_model::param_constrain(bool include_tp, bool include_gq, + const double* theta_unc, double* theta, + boost::ecuyer1988& rng) { using Eigen::VectorXd; VectorXd params_unc = VectorXd::Map(theta_unc, param_unc_num_); Eigen::VectorXd params; - model_->write_array(rng_, params_unc, params, include_tp, include_gq, + model_->write_array(rng, params_unc, params, include_tp, include_gq, &std::cout); Eigen::VectorXd::Map(theta, params.size()) = params; } -auto bs_model_rng::make_model_lambda(bool propto, bool jacobian) const { +auto bs_model::make_model_lambda(bool propto, bool jacobian) const { return [model = this->model_, propto, jacobian](auto& x) { // log_prob() requires non-const but doesn't modify its argument auto& params @@ -270,8 +268,8 @@ auto bs_model_rng::make_model_lambda(bool propto, bool jacobian) const { }; } -void bs_model_rng::log_density(bool propto, bool jacobian, - const double* theta_unc, double* val) const { +void bs_model::log_density(bool propto, bool jacobian, const double* theta_unc, + double* val) const { int N = param_unc_num_; if (propto) { Eigen::Map params_unc(theta_unc, N); @@ -291,9 +289,9 @@ void bs_model_rng::log_density(bool propto, bool jacobian, } } -void bs_model_rng::log_density_gradient(bool propto, bool jacobian, - const double* theta_unc, double* val, - double* grad) const { +void bs_model::log_density_gradient(bool propto, bool jacobian, + const double* theta_unc, double* val, + double* grad) const { #ifdef STAN_THREADS static thread_local stan::math::ChainableStack thread_instance; #endif @@ -303,9 +301,9 @@ void bs_model_rng::log_density_gradient(bool propto, bool jacobian, stan::math::gradient(logp, params_unc, *val, grad, grad + N); } -void bs_model_rng::log_density_hessian(bool propto, bool jacobian, - const double* theta_unc, double* val, - double* grad, double* hessian) const { +void bs_model::log_density_hessian(bool propto, bool jacobian, + const double* theta_unc, double* val, + double* grad, double* hessian) const { #ifdef STAN_THREADS static thread_local stan::math::ChainableStack thread_instance; #endif diff --git a/src/model_rng.hpp b/src/model_rng.hpp index 1ff7db78..59131937 100644 --- a/src/model_rng.hpp +++ b/src/model_rng.hpp @@ -5,6 +5,7 @@ #include #include #include +#include /** * This structure holds a pointer to a model, holds a pseudorandom @@ -12,7 +13,7 @@ * CSV format. Instances can be constructed with the C function * `construct()` and destroyed with the C function `destruct()`. */ -class bs_model_rng { +class bs_model { public: /** * Construct a model and random number generator with cached @@ -23,12 +24,12 @@ class bs_model_rng { * @param[in] chain_id number of gaps to skip in the pseudorandom * number generator for concurrent computations */ - bs_model_rng(const char* data_file, unsigned int seed, unsigned int chain_id); + bs_model(const char* data_file, unsigned int seed); /** * Destroy this object and free all of the memory allocated for it. */ - ~bs_model_rng() noexcept; + ~bs_model() noexcept; /** * Return the name of the model. This class manages the memory, @@ -114,7 +115,8 @@ class bs_model_rng { * @param[in,out] theta constrained parameters generated */ void param_constrain(bool include_tp, bool include_gq, - const double* theta_unc, double* theta); + const double* theta_unc, double* theta, + boost::ecuyer1988& rng); /** * Calculate the log density for the specified unconstrain @@ -181,9 +183,6 @@ class bs_model_rng { /** Stan model */ stan::model::model_base* model_; - /** pseudorandom number generator */ - boost::ecuyer1988 rng_; - /** name of the Stan model */ char* name_ = nullptr; @@ -224,4 +223,11 @@ class bs_model_rng { int param_unc_num_ = -1; }; +class bs_rng { + public: + bs_rng(unsigned int seed) : rng_(seed) { rng_.discard(1); } + + boost::ecuyer1988 rng_; +}; + #endif From a5d64db8978daae6b16327b7ca73a8091f33d6d2 Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Tue, 28 Mar 2023 14:56:35 -0400 Subject: [PATCH 02/20] Checkpoint --- python/bridgestan/model.py | 93 ++++++++++++++++++++++++++++++++------ src/bridgestan.cpp | 28 ++++++++++-- src/bridgestan.h | 36 ++++++++++++++- src/bridgestanR.cpp | 2 +- src/model_rng.hpp | 7 ++- 5 files changed, 141 insertions(+), 25 deletions(-) diff --git a/python/bridgestan/model.py b/python/bridgestan/model.py index 3e791d08..ecdb4ec1 100644 --- a/python/bridgestan/model.py +++ b/python/bridgestan/model.py @@ -41,7 +41,6 @@ def __init__( model_data: Optional[str] = None, *, seed: int = 1234, - chain_id: int = 0, ) -> None: """ Construct a StanModel object for a Stan model and data given @@ -50,9 +49,8 @@ def __init__( :param model_lib: A system path to compiled shared object. :param model_data: Either a JSON string literal or a system path to a data file in JSON format ending in ``.json``. - :param seed: A pseudo random number generator seed. - :param chain_id: A unique identifier for concurrent chains of - pseudorandom numbers. + :param seed: A pseudo random number generator seed. This is only used + if the model has RNG usage in the ``transformed data`` block. :raises FileNotFoundError or PermissionError: If ``model_lib`` is not readable or ``model_data`` is specified and not a path to a readable file. :raises RuntimeError: If there is an error instantiating the @@ -65,14 +63,12 @@ def __init__( self.stanlib = ctypes.CDLL(self.lib_path) self.data_path = model_data or "" self.seed = seed - self.chain_id = chain_id self._construct = self.stanlib.bs_construct self._construct.restype = ctypes.c_void_p self._construct.argtypes = [ ctypes.c_char_p, ctypes.c_uint, - ctypes.c_uint, star_star_char, ] @@ -82,7 +78,7 @@ def __init__( err = ctypes.pointer(ctypes.c_char_p()) self.model_rng = self._construct( - str.encode(self.data_path), self.seed, self.chain_id, err + str.encode(self.data_path), self.seed, err ) if not self.model_rng: @@ -132,6 +128,19 @@ def __init__( double_array, double_array, star_star_char, + ctypes.c_void_p, + ] + + self._param_constrain_seed = self.stanlib.bs_param_constrain_seed + self._param_constrain_seed.restype = ctypes.c_int + self._param_constrain_seed.argtypes = [ + ctypes.c_void_p, + ctypes.c_int, + ctypes.c_int, + double_array, + double_array, + star_star_char, + ctypes.c_uint, ] self._param_unconstrain = self.stanlib.bs_param_unconstrain @@ -201,7 +210,6 @@ def from_stan_file( stanc_args: List[str] = [], make_args: List[str] = [], seed: int = 1234, - chain_id: int = 0, ): """ Construct a StanModel instance from a ``.stan`` file, compiling if necessary. @@ -226,7 +234,7 @@ def from_stan_file( :raises RuntimeError: If compilation fails. """ result = compile_model(stan_file, stanc_args=stanc_args, make_args=make_args) - return cls(str(result), model_data, seed=seed, chain_id=chain_id) + return cls(str(result), model_data, seed=seed) def __del__(self) -> None: """ @@ -237,7 +245,7 @@ def __del__(self) -> None: def __repr__(self) -> str: data = f"{self.data_path!r}, " if self.data_path else "" - return f"StanModel({self.lib_path!r}, {data}seed={self.seed}, chain_id={self.chain_id})" + return f"StanModel({self.lib_path!r}, {data}seed={self.seed})" def name(self) -> str: """ @@ -326,6 +334,8 @@ def param_constrain( include_tp: bool = False, include_gq: bool = False, out: Optional[FloatArray] = None, + rng: Optional["StanRNG"] = None, + seed: Optional[int] = None, ) -> FloatArray: """ Return the constrained parameters derived from the specified @@ -354,13 +364,49 @@ def param_constrain( "Error: out must be same size as number of constrained parameters" ) err = ctypes.pointer(ctypes.c_char_p()) - rc = self._param_constrain( - self.model_rng, int(include_tp), int(include_gq), theta_unc, out, err - ) + + if seed is None: + if rng is None: + if include_gq: + raise ValueError( + "Error: must specify rng or seed when including generated quantities" + ) + rc = self._param_constrain( + self.model_rng, + int(include_tp), + int(include_gq), + theta_unc, + out, + None, + err, + ) + else: + rc = self._param_constrain( + self.model_rng, + int(include_tp), + int(include_gq), + theta_unc, + out, + rng.ptr, + err, + ) + else: + rc = self._param_constrain_seed( + self.model_rng, int(include_tp), int(include_gq), theta_unc, out, seed + ) if rc: raise self._handle_error(err.contents, "param_constrain") return out + def new_rng(self, seed: int) -> "StanRNG": + """ + Return a new PRNG for the model. + + :param seed: The seed for the PRNG. + :return: A new PRNG for the model. + """ + return StanRNG(self.stanlib, seed) + def param_unconstrain( self, theta: FloatArray, *, out: Optional[FloatArray] = None ) -> FloatArray: @@ -572,3 +618,24 @@ def _handle_error(self, err: ctypes.c_char_p, method: str) -> Exception: return RuntimeError(string) else: return RuntimeError(f"Unknown error in {method}. ") + + +class StanRNG: + def __init__(self, lib: ctypes.CDLL, seed: int) -> None: + self.stanlib = lib + + construct = self.stanlib.bs_construct_rng + construct.restype = ctypes.c_void_p + construct.argtypes = [ctypes.c_uint] + self.ptr = construct(seed) + + self._destruct = self.stanlib.bs_destruct_rng + self._destruct.restype = ctypes.c_int + self._destruct.argtypes = [ctypes.c_void_p] + + def __del__(self) -> None: + """ + Destroy the Stan model and free memory. + """ + if hasattr(self, "ptr") and hasattr(self, "_destruct"): + self._destruct(self.ptr) diff --git a/src/bridgestan.cpp b/src/bridgestan.cpp index 862d810b..8f5274dc 100644 --- a/src/bridgestan.cpp +++ b/src/bridgestan.cpp @@ -58,10 +58,10 @@ int bs_param_num(const bs_model* m, bool include_tp, bool include_gq) { int bs_param_unc_num(const bs_model* m) { return m->param_unc_num(); } int bs_param_constrain(bs_model* m, bool include_tp, bool include_gq, - const double* theta_unc, double* theta, boost::ecuyer1988& rng, + const double* theta_unc, double* theta, bs_rng* rng, char** error_msg) { try { - m->param_constrain(include_tp, include_gq, theta_unc, theta, rng); + m->param_constrain(include_tp, include_gq, theta_unc, theta, rng->rng_); return 0; } catch (const std::exception& e) { if (error_msg) { @@ -80,6 +80,14 @@ int bs_param_constrain(bs_model* m, bool include_tp, bool include_gq, return 1; } +int bs_param_constrain_seed(bs_model* m, bool include_tp, bool include_gq, + const double* theta_unc, double* theta, + unsigned int seed, char** error_msg) { + bs_rng rng(seed); + return bs_param_constrain(m, include_tp, include_gq, theta_unc, theta, &rng, + error_msg); +} + int bs_param_unconstrain(const bs_model* m, const double* theta, double* theta_unc, char** error_msg) { try { @@ -194,15 +202,25 @@ int bs_log_density_hessian(const bs_model* m, bool propto, bool jacobian, return -1; } -bs_rng* bs_rng_construct(unsigned int seed) { +bs_rng* bs_construct_rng(unsigned int seed) { try { return new bs_rng(seed); } catch (const std::exception& e) { - std::cerr << "bs_rng_construct(" << seed + std::cerr << "bs_construct_rng(" << seed << ") failed with exception: " << e.what() << std::endl; } catch (...) { - std::cerr << "bs_rng_construct(" << seed + std::cerr << "bs_construct_rng(" << seed << ") failed with unknown exception" << std::endl; } return nullptr; } + +int bs_destruct_rng(bs_rng* mr) { + try { + delete (mr); + return 0; + } catch (...) { + std::cerr << "destruct_rng() failed." << std::endl; + } + return -1; +} diff --git a/src/bridgestan.h b/src/bridgestan.h index ac5faebb..2c6f8eee 100644 --- a/src/bridgestan.h +++ b/src/bridgestan.h @@ -147,15 +147,42 @@ int bs_param_unc_num(const bs_model* m); * @param[in] include_gq `true` to include generated quantities * @param[in] theta_unc sequence of unconstrained parameters * @param[out] theta sequence of constrained parameters + * @param[in] rng pointer to pseudorandom number generator, should be created + * by `bs_construct_rng` * @param[out] error_msg a pointer to a string that will be allocated if there * is an error. This must later be freed by calling `bs_free_error_msg`. * @return code 0 if successful and code -1 if there is an exception * in the underlying Stan code */ int bs_param_constrain(bs_model* m, bool include_tp, bool include_gq, - const double* theta_unc, double* theta, + const double* theta_unc, double* theta, bs_rng* rng, char** error_msg); +/** + * Set the sequence of constrained parameters based on the specified + * unconstrained parameters, including transformed parameters and/or + * generated quantities as specified, and return a return code of 0 + * for success and -1 for failure. Parameter order is as declared + * in the Stan program, with multivariate parameters given in + * last-index-major order. + * + * @param[in] mr pointer to model and RNG structure + * @param[in] include_tp `true` to include transformed parameters + * @param[in] include_gq `true` to include generated quantities + * @param[in] theta_unc sequence of unconstrained parameters + * @param[out] theta sequence of constrained parameters + * @param[in] seed seed for pseudorandom number generator which will be created + * and destroyed during this call. See `bs_param_constrain` for an option with a + * persistent RNG. + * @param[out] error_msg a pointer to a string that will be allocated if there + * is an error. This must later be freed by calling `bs_free_error_msg`. + * @return code 0 if successful and code -1 if there is an exception + * in the underlying Stan code + */ +int bs_param_constrain_seed(bs_model* mr, bool include_tp, bool include_gq, + const double* theta_unc, double* theta, + unsigned int seed, char** error_msg); + /** * Set the sequence of unconstrained parameters based on the * specified constrained parameters, and return a return code of 0 @@ -273,7 +300,12 @@ int bs_log_density_hessian(const bs_model* m, bool propto, bool jacobian, */ bs_rng* bs_construct_rng(unsigned int seed); -// TODO destruct as well +/** + * Destruct an RNG object + * + * @param[in] rng pointer to RNG object + */ +int bs_destruct_rng(bs_rng* rng); #ifdef __cplusplus } diff --git a/src/bridgestanR.cpp b/src/bridgestanR.cpp index 94ecd833..fb1c01c7 100644 --- a/src/bridgestanR.cpp +++ b/src/bridgestanR.cpp @@ -39,7 +39,7 @@ void bs_param_constrain_R(bs_model** model, int* include_tp, int* include_gq, const double* theta_unc, double* theta, int* return_code, char** err_msg, void** err_ptr) { *return_code = bs_param_constrain(*model, *include_tp, *include_gq, theta_unc, - theta, err_msg); + theta, /* TODO(bmw) */ nullptr, err_msg); *err_ptr = static_cast(*err_msg); } void bs_param_unconstrain_R(bs_model** model, const double* theta, diff --git a/src/model_rng.hpp b/src/model_rng.hpp index 59131937..37034e1b 100644 --- a/src/model_rng.hpp +++ b/src/model_rng.hpp @@ -8,10 +8,9 @@ #include /** - * This structure holds a pointer to a model, holds a pseudorandom - * number generator, and holds pointers to the parameter names in - * CSV format. Instances can be constructed with the C function - * `construct()` and destroyed with the C function `destruct()`. + * This structure holds a pointer to a model and holds pointers to the parameter + * names in CSV format. Instances can be constructed with the C function + * `bs_construct()` and destroyed with the C function `destruct()`. */ class bs_model { public: From decdc1ebe2e6215ed256f7a23eb1be61d401701b Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Tue, 28 Mar 2023 16:45:04 -0400 Subject: [PATCH 03/20] Checkpointing; rebased on error msg PR to avoid stepping on myself --- python/bridgestan/model.py | 72 ++++++++++++++++++++------------------ src/bridgestan.cpp | 30 +++++++++++----- src/bridgestan.h | 12 ++++--- src/model_rng.cpp | 2 +- src/model_rng.hpp | 2 +- 5 files changed, 68 insertions(+), 50 deletions(-) diff --git a/python/bridgestan/model.py b/python/bridgestan/model.py index ecdb4ec1..f31643d2 100644 --- a/python/bridgestan/model.py +++ b/python/bridgestan/model.py @@ -77,9 +77,7 @@ def __init__( self._free_error.argtypes = [ctypes.c_char_p] err = ctypes.pointer(ctypes.c_char_p()) - self.model_rng = self._construct( - str.encode(self.data_path), self.seed, err - ) + self.model_rng = self._construct(str.encode(self.data_path), self.seed, err) if not self.model_rng: raise self._handle_error(err.contents, "bs_construct") @@ -127,8 +125,8 @@ def __init__( ctypes.c_int, double_array, double_array, - star_star_char, ctypes.c_void_p, + star_star_char, ] self._param_constrain_seed = self.stanlib.bs_param_constrain_seed @@ -139,8 +137,8 @@ def __init__( ctypes.c_int, double_array, double_array, - star_star_char, ctypes.c_uint, + star_star_char, ] self._param_unconstrain = self.stanlib.bs_param_unconstrain @@ -356,6 +354,15 @@ def param_constrain( shape as the return. :raises RuntimeError: If the C++ Stan model throws an exception. """ + if seed is None and rng is None: + if include_gq: + raise ValueError( + "Error: must specify rng or seed when including generated quantities" + ) + else: + # neither specified, but not doing gq, so use a fixed seed + seed = 0 + dims = self.param_num(include_tp=include_tp, include_gq=include_gq) if out is None: out = np.zeros(dims) @@ -365,35 +372,27 @@ def param_constrain( ) err = ctypes.pointer(ctypes.c_char_p()) - if seed is None: - if rng is None: - if include_gq: - raise ValueError( - "Error: must specify rng or seed when including generated quantities" - ) - rc = self._param_constrain( - self.model_rng, - int(include_tp), - int(include_gq), - theta_unc, - out, - None, - err, - ) - else: - rc = self._param_constrain( - self.model_rng, - int(include_tp), - int(include_gq), - theta_unc, - out, - rng.ptr, - err, - ) + if rng is not None: + rc = self._param_constrain( + self.model_rng, + int(include_tp), + int(include_gq), + theta_unc, + out, + rng.ptr, + err, + ) else: rc = self._param_constrain_seed( - self.model_rng, int(include_tp), int(include_gq), theta_unc, out, seed + self.model_rng, + int(include_tp), + int(include_gq), + theta_unc, + out, + seed, + err, ) + if rc: raise self._handle_error(err.contents, "param_constrain") return out @@ -626,16 +625,19 @@ def __init__(self, lib: ctypes.CDLL, seed: int) -> None: construct = self.stanlib.bs_construct_rng construct.restype = ctypes.c_void_p - construct.argtypes = [ctypes.c_uint] - self.ptr = construct(seed) + construct.argtypes = [ctypes.c_uint, star_star_char] + self.ptr = construct(seed, None) + + if not self.ptr: + raise RuntimeError("Failed to construct RNG.") self._destruct = self.stanlib.bs_destruct_rng self._destruct.restype = ctypes.c_int - self._destruct.argtypes = [ctypes.c_void_p] + self._destruct.argtypes = [ctypes.c_void_p, star_star_char] def __del__(self) -> None: """ Destroy the Stan model and free memory. """ if hasattr(self, "ptr") and hasattr(self, "_destruct"): - self._destruct(self.ptr) + self._destruct(self.ptr, None) diff --git a/src/bridgestan.cpp b/src/bridgestan.cpp index 8f5274dc..cf68d0ca 100644 --- a/src/bridgestan.cpp +++ b/src/bridgestan.cpp @@ -57,7 +57,7 @@ int bs_param_num(const bs_model* m, bool include_tp, bool include_gq) { int bs_param_unc_num(const bs_model* m) { return m->param_unc_num(); } -int bs_param_constrain(bs_model* m, bool include_tp, bool include_gq, +int bs_param_constrain(const bs_model* m, bool include_tp, bool include_gq, const double* theta_unc, double* theta, bs_rng* rng, char** error_msg) { try { @@ -80,7 +80,7 @@ int bs_param_constrain(bs_model* m, bool include_tp, bool include_gq, return 1; } -int bs_param_constrain_seed(bs_model* m, bool include_tp, bool include_gq, +int bs_param_constrain_seed(const bs_model* m, bool include_tp, bool include_gq, const double* theta_unc, double* theta, unsigned int seed, char** error_msg) { bs_rng rng(seed); @@ -202,25 +202,37 @@ int bs_log_density_hessian(const bs_model* m, bool propto, bool jacobian, return -1; } -bs_rng* bs_construct_rng(unsigned int seed) { +bs_rng* bs_construct_rng(unsigned int seed, char** error_msg) { try { return new bs_rng(seed); } catch (const std::exception& e) { - std::cerr << "bs_construct_rng(" << seed - << ") failed with exception: " << e.what() << std::endl; + if (error_msg) { + std::stringstream error; + error << "construct_rng() failed with exception: " << e.what() + << std::endl; + *error_msg = strdup(error.str().c_str()); + } } catch (...) { - std::cerr << "bs_construct_rng(" << seed - << ") failed with unknown exception" << std::endl; + if (error_msg) { + std::stringstream error; + error << "construct_rng() failed with unknown exception" << std::endl; + *error_msg = strdup(error.str().c_str()); + } } + return nullptr; } -int bs_destruct_rng(bs_rng* mr) { +int bs_destruct_rng(bs_rng* mr, char** error_msg) { try { delete (mr); return 0; } catch (...) { - std::cerr << "destruct_rng() failed." << std::endl; + if (error_msg) { + std::stringstream error; + error << "destruct_rng() failed." << std::endl; + *error_msg = strdup(error.str().c_str()); + } } return -1; } diff --git a/src/bridgestan.h b/src/bridgestan.h index 2c6f8eee..e330e78a 100644 --- a/src/bridgestan.h +++ b/src/bridgestan.h @@ -154,7 +154,7 @@ int bs_param_unc_num(const bs_model* m); * @return code 0 if successful and code -1 if there is an exception * in the underlying Stan code */ -int bs_param_constrain(bs_model* m, bool include_tp, bool include_gq, +int bs_param_constrain(const bs_model* m, bool include_tp, bool include_gq, const double* theta_unc, double* theta, bs_rng* rng, char** error_msg); @@ -179,7 +179,7 @@ int bs_param_constrain(bs_model* m, bool include_tp, bool include_gq, * @return code 0 if successful and code -1 if there is an exception * in the underlying Stan code */ -int bs_param_constrain_seed(bs_model* mr, bool include_tp, bool include_gq, +int bs_param_constrain_seed(const bs_model* mr, bool include_tp, bool include_gq, const double* theta_unc, double* theta, unsigned int seed, char** error_msg); @@ -297,15 +297,19 @@ int bs_log_density_hessian(const bs_model* m, bool propto, bool jacobian, * Construct an RNG object to be used for param_constrain * * @param[in] seed seed for the RNG + * @param[out] error_msg a pointer to a string that will be allocated if there + * is an error. This must later be freed by calling `bs_free_error_msg`. */ -bs_rng* bs_construct_rng(unsigned int seed); +bs_rng* bs_construct_rng(unsigned int seed, char** error_msg); /** * Destruct an RNG object * * @param[in] rng pointer to RNG object + * @param[out] error_msg a pointer to a string that will be allocated if there + * is an error. This must later be freed by calling `bs_free_error_msg`. */ -int bs_destruct_rng(bs_rng* rng); +int bs_destruct_rng(bs_rng* rng, char** error_msg); #ifdef __cplusplus } diff --git a/src/model_rng.cpp b/src/model_rng.cpp index ce0428da..48f20ec1 100644 --- a/src/model_rng.cpp +++ b/src/model_rng.cpp @@ -237,7 +237,7 @@ void bs_model::param_unconstrain_json(const char* json, void bs_model::param_constrain(bool include_tp, bool include_gq, const double* theta_unc, double* theta, - boost::ecuyer1988& rng) { + boost::ecuyer1988& rng) const { using Eigen::VectorXd; VectorXd params_unc = VectorXd::Map(theta_unc, param_unc_num_); Eigen::VectorXd params; diff --git a/src/model_rng.hpp b/src/model_rng.hpp index 37034e1b..2c360e00 100644 --- a/src/model_rng.hpp +++ b/src/model_rng.hpp @@ -115,7 +115,7 @@ class bs_model { */ void param_constrain(bool include_tp, bool include_gq, const double* theta_unc, double* theta, - boost::ecuyer1988& rng); + boost::ecuyer1988& rng) const; /** * Calculate the log density for the specified unconstrain From 5135cedc49c62233a895fa1c451bb02681a55122 Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Tue, 28 Mar 2023 16:51:37 -0400 Subject: [PATCH 04/20] Update doc comments --- python/bridgestan/model.py | 10 +++++++++- src/bridgestan.h | 12 +++++++----- src/bridgestanR.cpp | 4 +++- src/model_rng.hpp | 8 +++++++- 4 files changed, 26 insertions(+), 8 deletions(-) diff --git a/python/bridgestan/model.py b/python/bridgestan/model.py index f31643d2..eed9dc3c 100644 --- a/python/bridgestan/model.py +++ b/python/bridgestan/model.py @@ -349,9 +349,17 @@ def param_constrain( provided, it must have shape `(D, )`, where `D` is the number of constrained parameters. If not provided or `None`, a freshly allocated array is returned. + :param rng: A ``StanRNG`` object to use for generating random + numbers. One of ``rng`` or ``seed`` must be specified if + ``include_gq`` is ``True``. + :param seed: A pseudo random number generator seed. One of + ``rng`` or ``seed`` must be specified if ``include_gq`` + is ``True``. :return: The constrained parameter array. - :raises ValueError: If `out` is specified and is not the same + :raises ValueError: If ``out`` is specified and is not the same shape as the return. + :raises ValueError: If neither ``rng`` nor ``seed`` is specified + and ``include_gq`` is ``True``. :raises RuntimeError: If the C++ Stan model throws an exception. """ if seed is None and rng is None: diff --git a/src/bridgestan.h b/src/bridgestan.h index e330e78a..c7f1366c 100644 --- a/src/bridgestan.h +++ b/src/bridgestan.h @@ -179,9 +179,9 @@ int bs_param_constrain(const bs_model* m, bool include_tp, bool include_gq, * @return code 0 if successful and code -1 if there is an exception * in the underlying Stan code */ -int bs_param_constrain_seed(const bs_model* mr, bool include_tp, bool include_gq, - const double* theta_unc, double* theta, - unsigned int seed, char** error_msg); +int bs_param_constrain_seed(const bs_model* mr, bool include_tp, + bool include_gq, const double* theta_unc, + double* theta, unsigned int seed, char** error_msg); /** * Set the sequence of unconstrained parameters based on the @@ -294,7 +294,9 @@ int bs_log_density_hessian(const bs_model* m, bool propto, bool jacobian, double* hessian, char** error_msg); /** - * Construct an RNG object to be used for param_constrain + * Construct an RNG object to be used in `bs_param_constrain`. + * This object is not thread safe and should be constructed and + * destructed for each thread. * * @param[in] seed seed for the RNG * @param[out] error_msg a pointer to a string that will be allocated if there @@ -303,7 +305,7 @@ int bs_log_density_hessian(const bs_model* m, bool propto, bool jacobian, bs_rng* bs_construct_rng(unsigned int seed, char** error_msg); /** - * Destruct an RNG object + * Destruct an RNG object. * * @param[in] rng pointer to RNG object * @param[out] error_msg a pointer to a string that will be allocated if there diff --git a/src/bridgestanR.cpp b/src/bridgestanR.cpp index fb1c01c7..f78c84b6 100644 --- a/src/bridgestanR.cpp +++ b/src/bridgestanR.cpp @@ -38,8 +38,10 @@ void bs_param_unc_num_R(bs_model** model, int* num_out) { void bs_param_constrain_R(bs_model** model, int* include_tp, int* include_gq, const double* theta_unc, double* theta, int* return_code, char** err_msg, void** err_ptr) { + /* TODO(bmw): accept seed as argument, consider persistent version */ + bs_rng rng(0); *return_code = bs_param_constrain(*model, *include_tp, *include_gq, theta_unc, - theta, /* TODO(bmw) */ nullptr, err_msg); + theta, &rng, err_msg); *err_ptr = static_cast(*err_msg); } void bs_param_unconstrain_R(bs_model** model, const double* theta, diff --git a/src/model_rng.hpp b/src/model_rng.hpp index 2c360e00..6d1feb82 100644 --- a/src/model_rng.hpp +++ b/src/model_rng.hpp @@ -10,7 +10,7 @@ /** * This structure holds a pointer to a model and holds pointers to the parameter * names in CSV format. Instances can be constructed with the C function - * `bs_construct()` and destroyed with the C function `destruct()`. + * `bs_construct()` and destroyed with the C function `bs_destruct()`. */ class bs_model { public: @@ -222,6 +222,12 @@ class bs_model { int param_unc_num_ = -1; }; +/** + * A wrapper around the Boost random number generator required + * by the Stan model's write_array methods. Instances can be + * constructed with the C function `bs_construct_rng()` and destroyed + * with the C function `bs_destruct_rng()`. + */ class bs_rng { public: bs_rng(unsigned int seed) : rng_(seed) { rng_.discard(1); } From 4efe5d698f4abb346398d6b6778367d143de737e Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Wed, 29 Mar 2023 11:13:04 -0400 Subject: [PATCH 05/20] Update Python tests --- python/test/test_stanmodel.py | 21 ++++++++++++++++----- src/bridgestan.cpp | 15 +++++++-------- src/bridgestan.h | 5 +---- src/bridgestanR.cpp | 4 ++-- src/bridgestanR.h | 2 +- test_models/full/full.stan | 5 ++++- 6 files changed, 31 insertions(+), 21 deletions(-) diff --git a/python/test/test_stanmodel.py b/python/test/test_stanmodel.py index a641f63a..32914d77 100644 --- a/python/test/test_stanmodel.py +++ b/python/test/test_stanmodel.py @@ -43,7 +43,6 @@ def test_constructor(): b4 = bs.StanModel(throw_data_so) - def test_name(): std_so = str(STAN_FOLDER / "stdnormal" / "stdnormal_model.so") b = bs.StanModel(std_so) @@ -57,7 +56,6 @@ def test_model_info(): assert "BridgeStan version: 1." in b.model_info() - def test_param_num(): full_so = str(STAN_FOLDER / "full" / "full_model.so") b = bs.StanModel(full_so) @@ -210,13 +208,26 @@ def test_param_constrain(): full_so = str(STAN_FOLDER / "full" / "full_model.so") bridge2 = bs.StanModel(full_so) + rng = bridge2.new_rng(seed=1234) np.testing.assert_equal(1, bridge2.param_constrain(a).size) np.testing.assert_equal(2, bridge2.param_constrain(a, include_tp=True).size) - np.testing.assert_equal(3, bridge2.param_constrain(a, include_gq=True).size) np.testing.assert_equal( - 4, bridge2.param_constrain(a, include_tp=True, include_gq=True).size + 3, bridge2.param_constrain(a, include_gq=True, rng=rng).size ) + np.testing.assert_equal( + 4, bridge2.param_constrain(a, include_tp=True, include_gq=True, rng=rng).size + ) + + # reproducibility test + np.testing.assert_equal( + bridge2.param_constrain(a, include_gq=True, seed=1234), + bridge2.param_constrain(a, include_gq=True, seed=1234), + ) + + # test error if neither seed or rng is provided + with pytest.raises(ValueError): + bridge2.param_constrain(a, include_gq=True) # out tests, matched and mismatched scratch = np.zeros(16) @@ -240,7 +251,7 @@ def test_param_constrain(): bridge3 = bs.StanModel(throw_gq_so) bridge3.param_constrain(y, include_gq=False) with pytest.raises(RuntimeError, match="find this text: gqfails"): - bridge3.param_constrain(y, include_gq=True) + bridge3.param_constrain(y, include_gq=True, seed=123) def test_param_unconstrain(): diff --git a/src/bridgestan.cpp b/src/bridgestan.cpp index cf68d0ca..9f880262 100644 --- a/src/bridgestan.cpp +++ b/src/bridgestan.cpp @@ -10,14 +10,13 @@ int bs_patch_version = BRIDGESTAN_PATCH; #include bs_model* bs_construct(const char* data_file, unsigned int seed, - unsigned int chain_id, char** error_msg) { + char** error_msg) { try { return new bs_model(data_file, seed); } catch (const std::exception& e) { if (error_msg) { std::stringstream error; - error << "construct(" << data_file << ", " << seed << ", " << chain_id - << ")" + error << "construct(" << data_file << ", " << seed << ")" << " failed with exception: " << e.what() << std::endl; *error_msg = strdup(error.str().c_str()); } @@ -25,8 +24,7 @@ bs_model* bs_construct(const char* data_file, unsigned int seed, if (error_msg) { std::stringstream error; - error << "construct(" << data_file << ", " << seed << ", " << chain_id - << ")" + error << "construct(" << data_file << ", " << seed << ")" << " failed with unknown exception" << std::endl; *error_msg = strdup(error.str().c_str()); } @@ -80,9 +78,10 @@ int bs_param_constrain(const bs_model* m, bool include_tp, bool include_gq, return 1; } -int bs_param_constrain_seed(const bs_model* m, bool include_tp, bool include_gq, - const double* theta_unc, double* theta, - unsigned int seed, char** error_msg) { +int bs_param_constrain_seed(const bs_model* m, bool include_tp, + bool include_gq, const double* theta_unc, + double* theta, unsigned int seed, + char** error_msg) { bs_rng rng(seed); return bs_param_constrain(m, include_tp, include_gq, theta_unc, theta, &rng, error_msg); diff --git a/src/bridgestan.h b/src/bridgestan.h index c7f1366c..166f20c3 100644 --- a/src/bridgestan.h +++ b/src/bridgestan.h @@ -23,15 +23,12 @@ extern int bs_patch_version; * path to JSON-encoded data file (must end with ".json"), or * a JSON string literal. * @param[in] seed seed for PRNG - * @param[in] chain_id identifier for concurrent sequence of PRNG - * draws * @param[out] error_msg a pointer to a string that will be allocated if there * is an error. This must later be freed by calling `bs_free_error_msg`. * @return pointer to constructed model or `nullptr` if construction * fails */ -bs_model* bs_construct(const char* data_file, unsigned int seed, - unsigned int chain_id, char** error_msg); +bs_model* bs_construct(const char* data_file, unsigned int seed, char** error_msg); /** * Destroy the model. diff --git a/src/bridgestanR.cpp b/src/bridgestanR.cpp index f78c84b6..a0c43f2c 100644 --- a/src/bridgestanR.cpp +++ b/src/bridgestanR.cpp @@ -1,9 +1,9 @@ #include "bridgestanR.h" #include "bridgestan.h" -void bs_construct_R(char** data, int* rng, int* chain, bs_model** ptr_out, +void bs_construct_R(char** data, int* rng, bs_model** ptr_out, char** err_msg, void** err_ptr) { - *ptr_out = bs_construct(*data, *rng, *chain, err_msg); + *ptr_out = bs_construct(*data, *rng, err_msg); *err_ptr = static_cast(*err_msg); } void bs_free_error_msg_R(void** err_msg) { diff --git a/src/bridgestanR.h b/src/bridgestanR.h index 9eb6c16f..f3a04d53 100644 --- a/src/bridgestanR.h +++ b/src/bridgestanR.h @@ -12,7 +12,7 @@ typedef int bool; // Shim to convert to R interface requirement of void with pointer args // All calls directly delegated to versions without _R suffix -void bs_construct_R(char** data, int* rng, int* chain, bs_model** ptr_out, +void bs_construct_R(char** data, int* rng, bs_model** ptr_out, char** err_msg, void** err_ptr); void bs_version_R(int* major, int* minor, int* patch); diff --git a/test_models/full/full.stan b/test_models/full/full.stan index f5c64da9..6897448d 100644 --- a/test_models/full/full.stan +++ b/test_models/full/full.stan @@ -5,5 +5,8 @@ transformed parameters { real b = exp(a); } generated quantities { - vector[2] c = [1, 2]'; + array[2] int d; + for (i in 1:2) { + d[i] = categorical_rng([0.1, 0.2, 0.2, 0.3, 0.2]'); + } } From 690bf629b1da354c2f546b843c50e9ec7209e645 Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Wed, 29 Mar 2023 11:33:00 -0400 Subject: [PATCH 06/20] Start julia --- julia/src/BridgeStan.jl | 5 +- julia/src/model.jl | 101 ++++++++++++++++++++++++++++++++++------ 2 files changed, 89 insertions(+), 17 deletions(-) diff --git a/julia/src/BridgeStan.jl b/julia/src/BridgeStan.jl index 79ef77a7..091f9ee2 100644 --- a/julia/src/BridgeStan.jl +++ b/julia/src/BridgeStan.jl @@ -28,7 +28,7 @@ include("model.jl") include("compile.jl") """ - StanModel(;stan_file, stanc_args=[], make_args=[], data="", seed=204, chain_id=0) + StanModel(;stan_file, stanc_args=[], make_args=[], data="", seed=204) Construct a StanModel instance from a `.stan` file, compiling if necessary. @@ -40,7 +40,6 @@ StanModel(; make_args::AbstractVector{String} = String[], data::String = "", seed = 204, - chain_id = 0, -) = StanModel(compile_model(stan_file; stanc_args, make_args), data, seed, chain_id) +) = StanModel(compile_model(stan_file; stanc_args, make_args), data, seed) end diff --git a/julia/src/model.jl b/julia/src/model.jl index 284c41dc..96a6139c 100644 --- a/julia/src/model.jl +++ b/julia/src/model.jl @@ -1,5 +1,7 @@ -mutable struct StanModelStruct end +struct StanModelStruct end + +mutable struct StanRNGStruct end # utility macro to annotate a field as const only if supported @eval macro $(Symbol("const"))(x) @@ -11,19 +13,19 @@ mutable struct StanModelStruct end end """ - StanModel(lib, datafile="", seed=204, chain_id=0) + StanModel(lib, datafile="", seed=204) A StanModel instance encapsulates a Stan model instantiated with data. Construct a Stan model from the supplied library file path and data. Data should either be a string containing a JSON string literal or a path to a data file ending in `.json`. -If seed or chain_id are supplied, these are used to initialize the RNG used by the model. +If seed is supplied, it is used to initialize the RNG used by the model's constructor. - StanModel(;stan_file, data="", seed=204, chain_id=0) + StanModel(;stan_file, data="", seed=204) Construct a `StanModel` instance from a `.stan` file, compiling if necessary. - StanModel(;stan_file, stanc_args=[], make_args=[], data="", seed=204, chain_id=0) + StanModel(;stan_file, stanc_args=[], make_args=[], data="", seed=204) Construct a `StanModel` instance from a `.stan` file. Compilation occurs if no shared object file exists for the supplied Stan file or @@ -36,11 +38,9 @@ mutable struct StanModel stanmodel::Ptr{StanModelStruct} @const data::String @const seed::UInt32 - @const chain_id::UInt32 - function StanModel(lib::String, data::String = "", seed = 204, chain_id = 0) + function StanModel(lib::String, data::String = "", seed = 204) seed = convert(UInt32, seed) - chain_id = convert(UInt32, chain_id) if !isfile(lib) throw(SystemError("Dynamic library file not found")) @@ -64,17 +64,16 @@ mutable struct StanModel stanmodel = ccall( Libc.Libdl.dlsym(lib, "bs_construct"), Ptr{StanModelStruct}, - (Cstring, UInt32, UInt32, Ref{Cstring}), + (Cstring, UInt32, Ref{Cstring}), data, seed, - chain_id, err, ) if stanmodel == C_NULL error(handle_error(lib, err, "bs_construct")) end - sm = new(lib, stanmodel, data, seed, chain_id) + sm = new(lib, stanmodel, data, seed) function f(sm) ccall( @@ -89,6 +88,44 @@ mutable struct StanModel end end + +mutable struct StanRNG + lib::Ptr{Nothing} + rng::Ptr{StanRNGStruct} + seed::UInt32 + + function StanRNG(lib::Ptr{Nothing}, seed) + seed = convert(UInt32, seed) + + + err = Ref{Cstring}() + + rng = ccall( + Libc.Libdl.dlsym(lib, "bs_construct_rng"), + Ptr{StanModelStruct}, + (UInt32, Ref{Cstring}), + seed, + err, + ) + if rng == C_NULL + error(_handle_error(lib, err, "bs_construct_rng")) + end + + stanrng = new(lib, rng, data, seed) + + function f(stanrng) + ccall( + Libc.Libdl.dlsym(stanrng.lib, "bs_destruct_rng"), + UInt32, + (Ptr{StanModelStruct},Ref{Cstring}), + stanrng.rng, C_NULL + ) + end + + finalizer(f, stanrng) + end +end + """ name(sm) @@ -234,6 +271,8 @@ function param_constrain!( out::Vector{Float64}; include_tp = false, include_gq = false, + seed::Union{Int, Nothing} = nothing, + rng::Union{StanRNG, Nothing} = nothing, ) dims = param_num(sm; include_tp = include_tp, include_gq = include_gq) if length(out) != dims @@ -241,18 +280,50 @@ function param_constrain!( DimensionMismatch("out must be same size as number of constrained parameters"), ) end + + if seed === nothing && rng === nothing + if include_gq + throw( + ArgumentError( + "Must provide either a seed or an RNG when including generated quantities", + ), + ) + else + seed = 0 + end + end + + err = Ref{Cstring}() - rc = ccall( + if rng !== nothing + + rc = ccall( Libc.Libdl.dlsym(sm.lib, "bs_param_constrain"), Cint, - (Ptr{StanModelStruct}, Cint, Cint, Ref{Cdouble}, Ref{Cdouble}, Ref{Cstring}), + (Ptr{StanModelStruct}, Cint, Cint, Ref{Cdouble}, Ref{Cdouble}, Ptr{StanRNGStruct},Ref{Cstring}), sm.stanmodel, include_tp, include_gq, theta_unc, out, + rng, err, ) + else + seed = convert(UInt32, seed) + rc = ccall( + Libc.Libdl.dlsym(sm.lib, "bs_param_constrain"), + Cint, + (Ptr{StanModelStruct}, Cint, Cint, Ref{Cdouble}, Ref{Cdouble}, Cuint, Ref{Cstring}), + sm.stanmodel, + include_tp, + include_gq, + theta_unc, + out, + seed, + err, + ) + end if rc != 0 error(handle_error(sm.lib, err, "param_constrain")) end @@ -277,9 +348,11 @@ function param_constrain( theta_unc::Vector{Float64}; include_tp = false, include_gq = false, + seed::Union{Int, Nothing} = nothing, + rng::Union{StanRNG, Nothing} = nothing, ) out = zeros(param_num(sm, include_tp = include_tp, include_gq = include_gq)) - param_constrain!(sm, theta_unc, out; include_tp = include_tp, include_gq = include_gq) + param_constrain!(sm, theta_unc, out; include_tp = include_tp, include_gq = include_gq, seed=seed, rng=rng) end """ From a82cd50b1fb7cb47a9be0db7d5506cb07b2cbf8b Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Wed, 29 Mar 2023 11:59:06 -0400 Subject: [PATCH 07/20] Update Julia tests, doc --- docs/languages/julia.md | 68 +++++++++++++++++++++++++-------------- julia/docs/src/julia.md | 1 + julia/src/BridgeStan.jl | 3 +- julia/src/model.jl | 29 +++++++++++------ julia/test/model_tests.jl | 13 ++++++-- 5 files changed, 78 insertions(+), 36 deletions(-) diff --git a/docs/languages/julia.md b/docs/languages/julia.md index 85f685e9..96a85412 100644 --- a/docs/languages/julia.md +++ b/docs/languages/julia.md @@ -116,27 +116,27 @@ An example program is provided alongside the Julia interface code in `example.jl ```julia -StanModel(lib, datafile="", seed=204, chain_id=0) +StanModel(lib, datafile="", seed=204) ``` A StanModel instance encapsulates a Stan model instantiated with data. -Construct a Stan model from the supplied library file path and data. Data should either be a string containing a JSON string literal or a path to a data file ending in `.json`. If seed or chain_id are supplied, these are used to initialize the RNG used by the model. +Construct a Stan model from the supplied library file path and data. Data should either be a string containing a JSON string literal or a path to a data file ending in `.json`. If seed is supplied, it is used to initialize the RNG used by the model's constructor. ``` -StanModel(;stan_file, data="", seed=204, chain_id=0) +StanModel(;stan_file, data="", seed=204) ``` Construct a `StanModel` instance from a `.stan` file, compiling if necessary. ``` -StanModel(;stan_file, stanc_args=[], make_args=[], data="", seed=204, chain_id=0) +StanModel(;stan_file, stanc_args=[], make_args=[], data="", seed=204) ``` Construct a `StanModel` instance from a `.stan` file. Compilation occurs if no shared object file exists for the supplied Stan file or if a shared object file exists and the Stan file has changed since last compilation. This is equivalent to calling `compile_model` and then the original constructor of `StanModel`. -source
+source
# **`BridgeStan.log_density`** — *Function*. @@ -152,7 +152,7 @@ Return the log density of the specified unconstrained parameters. This calculation drops constant terms that do not depend on the parameters if `propto` is `true` and includes change of variables terms for constrained parameters if `jacobian` is `true`. -source
+source
# **`BridgeStan.log_density_gradient`** — *Function*. @@ -170,7 +170,7 @@ This calculation drops constant terms that do not depend on the parameters if `p This allocates new memory for the gradient output each call. See `log_density_gradient!` for a version which allows re-using existing memory. -source
+source
# **`BridgeStan.log_density_hessian`** — *Function*. @@ -188,7 +188,7 @@ This calculation drops constant terms that do not depend on the parameters if `p This allocates new memory for the gradient and Hessian output each call. See `log_density_gradient!` for a version which allows re-using existing memory. -source
+source
# **`BridgeStan.param_constrain`** — *Function*. @@ -196,17 +196,19 @@ This allocates new memory for the gradient and Hessian output each call. See `lo ```julia -param_constrain(sm, theta_unc, out; include_tp=false, include_gq=false) +param_constrain(sm, theta_unc, out; include_tp=false, include_gq=false, seed=nothing, rng=nothing) ``` Returns a vector constrained parameters given unconstrained parameters. Additionally (if `include_tp` and `include_gq` are set, respectively) returns transformed parameters and generated quantities. +If `include_gq` is set, then either `seed` or `rng` must be provided. + This allocates new memory for the output each call. See `param_constrain!` for a version which allows re-using existing memory. This is the inverse of `param_unconstrain`. -source
+source
# **`BridgeStan.param_unconstrain`** — *Function*. @@ -226,7 +228,7 @@ This allocates new memory for the output each call. See `param_unconstrain!` for This is the inverse of `param_constrain`. -source
+source
# **`BridgeStan.param_unconstrain_json`** — *Function*. @@ -244,7 +246,7 @@ The JSON is expected to be in the [JSON Format for CmdStan](https://mc-stan.org/ This allocates new memory for the output each call. See `param_unconstrain_json!` for a version which allows re-using existing memory. -source
+source
# **`BridgeStan.name`** — *Function*. @@ -258,7 +260,7 @@ name(sm) Return the name of the model `sm` -source
+source
# **`BridgeStan.model_info`** — *Function*. @@ -274,7 +276,7 @@ Return information about the model `sm`. This includes the Stan version and important compiler flags. -source
+source
# **`BridgeStan.param_num`** — *Function*. @@ -290,7 +292,7 @@ Return the number of (constrained) parameters in the model. This is the total of all the sizes of items declared in the `parameters` block of the model. If `include_tp` or `include_gq` are true, items declared in the `transformed parameters` and `generate quantities` blocks are included, respectively. -source
+source
# **`BridgeStan.param_unc_num`** — *Function*. @@ -306,7 +308,7 @@ Return the number of unconstrained parameters in the model. This function is mainly different from `param_num` when variables are declared with constraints. For example, `simplex[5]` has a constrained size of 5, but an unconstrained size of 4. -source
+source
# **`BridgeStan.param_names`** — *Function*. @@ -324,7 +326,7 @@ For containers, indexes are separated by periods (.). For example, the scalar `a` has indexed name `"a"`, the vector entry `a[1]` has indexed name `"a.1"` and the matrix entry `a[2, 3]` has indexed names `"a.2.3"`. Parameter order of the output is column major and more generally last-index major for containers. -source
+source
# **`BridgeStan.param_unc_names`** — *Function*. @@ -340,7 +342,7 @@ Return the indexed names of the unconstrained parameters. For example, a scalar unconstrained parameter `b` has indexed name `b` and a vector entry `b[3]` has indexed name `b.3`. -source
+source
# **`BridgeStan.log_density_gradient!`** — *Function*. @@ -358,7 +360,7 @@ This calculation drops constant terms that do not depend on the parameters if `p The gradient is stored in the vector `out`, and a reference is returned. See `log_density_gradient` for a version which allocates fresh memory. -source
+source
# **`BridgeStan.log_density_hessian!`** — *Function*. @@ -376,7 +378,7 @@ This calculation drops constant terms that do not depend on the parameters if `p The gradient is stored in the vector `out_grad` and the Hessian is stored in `out_hess` and references are returned. See `log_density_hessian` for a version which allocates fresh memory. -source
+source
# **`BridgeStan.param_constrain!`** — *Function*. @@ -384,17 +386,19 @@ The gradient is stored in the vector `out_grad` and the Hessian is stored in `ou ```julia -param_constrain!(sm, theta_unc, out; include_tp=false, include_gq=false) +param_constrain!(sm, theta_unc, out; include_tp=false, include_gq=false, seed=nothing, rng=nothing) ``` Returns a vector constrained parameters given unconstrained parameters. Additionally (if `include_tp` and `include_gq` are set, respectively) returns transformed parameters and generated quantities. +If `include_gq` is set, then either `seed` or `rng` must be provided. + The result is stored in the vector `out`, and a reference is returned. See `param_constrain` for a version which allocates fresh memory. This is the inverse of `param_unconstrain!`. -source
+source
# **`BridgeStan.param_unconstrain!`** — *Function*. @@ -414,7 +418,7 @@ The result is stored in the vector `out`, and a reference is returned. See `para This is the inverse of `param_constrain!`. -source
+source
# **`BridgeStan.param_unconstrain_json!`** — *Function*. @@ -432,7 +436,23 @@ The JSON is expected to be in the [JSON Format for CmdStan](https://mc-stan.org/ The result is stored in the vector `out`, and a reference is returned. See `param_unconstrain_json` for a version which allocates fresh memory. -source
+source
+ +# +**`BridgeStan.StanRNG`** — *Type*. + + + +```julia +StanRNG(sm::StanModel, seed) +``` + +Construct a StanRNG instance from a `StanModel` instance and a seed. This can be used in the `param_constrain` and `param_constrain!` methods when using the generated quantities block. + +This object is not thread-safe, one should be created per thread. + + +source
diff --git a/julia/docs/src/julia.md b/julia/docs/src/julia.md index 6a8f55f2..5f99b949 100644 --- a/julia/docs/src/julia.md +++ b/julia/docs/src/julia.md @@ -89,6 +89,7 @@ log_density_hessian! param_constrain! param_unconstrain! param_unconstrain_json! +StanRNG ``` ### Compilation utilities diff --git a/julia/src/BridgeStan.jl b/julia/src/BridgeStan.jl index 091f9ee2..5ad0322d 100644 --- a/julia/src/BridgeStan.jl +++ b/julia/src/BridgeStan.jl @@ -22,7 +22,8 @@ export StanModel, log_density_hessian, get_bridgestan_path, set_bridgestan_path!, - compile_model + compile_model, + StanRNG include("model.jl") include("compile.jl") diff --git a/julia/src/model.jl b/julia/src/model.jl index 96a6139c..fb6008dd 100644 --- a/julia/src/model.jl +++ b/julia/src/model.jl @@ -88,30 +88,37 @@ mutable struct StanModel end end +""" + StanRNG(sm::StanModel, seed) + +Construct a StanRNG instance from a `StanModel` instance and a seed. +This can be used in the `param_constrain` and `param_constrain!` methods +when using the generated quantities block. +This object is not thread-safe, one should be created per thread. +""" mutable struct StanRNG lib::Ptr{Nothing} rng::Ptr{StanRNGStruct} seed::UInt32 - function StanRNG(lib::Ptr{Nothing}, seed) + function StanRNG(sm::StanModel, seed) seed = convert(UInt32, seed) err = Ref{Cstring}() - rng = ccall( - Libc.Libdl.dlsym(lib, "bs_construct_rng"), + Libc.Libdl.dlsym(sm.lib, "bs_construct_rng"), Ptr{StanModelStruct}, (UInt32, Ref{Cstring}), seed, err, ) if rng == C_NULL - error(_handle_error(lib, err, "bs_construct_rng")) + error(_handle_error(sm.lib, err, "bs_construct_rng")) end - stanrng = new(lib, rng, data, seed) + stanrng = new(sm.lib, rng, seed) function f(stanrng) ccall( @@ -254,12 +261,14 @@ function param_unc_names(sm::StanModel) end """ - param_constrain!(sm, theta_unc, out; include_tp=false, include_gq=false) + param_constrain!(sm, theta_unc, out; include_tp=false, include_gq=false, seed=nothing, rng=nothing) Returns a vector constrained parameters given unconstrained parameters. Additionally (if `include_tp` and `include_gq` are set, respectively) returns transformed parameters and generated quantities. +If `include_gq` is set, then either `seed` or `rng` must be provided. + The result is stored in the vector `out`, and a reference is returned. See `param_constrain` for a version which allocates fresh memory. @@ -306,13 +315,13 @@ function param_constrain!( include_gq, theta_unc, out, - rng, + rng.rng, err, ) else seed = convert(UInt32, seed) rc = ccall( - Libc.Libdl.dlsym(sm.lib, "bs_param_constrain"), + Libc.Libdl.dlsym(sm.lib, "bs_param_constrain_seed"), Cint, (Ptr{StanModelStruct}, Cint, Cint, Ref{Cdouble}, Ref{Cdouble}, Cuint, Ref{Cstring}), sm.stanmodel, @@ -331,12 +340,14 @@ function param_constrain!( end """ - param_constrain(sm, theta_unc, out; include_tp=false, include_gq=false) + param_constrain(sm, theta_unc, out; include_tp=false, include_gq=false, seed=nothing, rng=nothing) Returns a vector constrained parameters given unconstrained parameters. Additionally (if `include_tp` and `include_gq` are set, respectively) returns transformed parameters and generated quantities. +If `include_gq` is set, then either `seed` or `rng` must be provided. + This allocates new memory for the output each call. See `param_constrain!` for a version which allows re-using existing memory. diff --git a/julia/test/model_tests.jl b/julia/test/model_tests.jl index ed76b871..feed1660 100644 --- a/julia/test/model_tests.jl +++ b/julia/test/model_tests.jl @@ -176,13 +176,21 @@ end model2 = load_test_model("full", false) + rng = StanRNG(model2, 1234) @test 1 == length(BridgeStan.param_constrain(model2, a)) @test 2 == length(BridgeStan.param_constrain(model2, a; include_tp = true)) - @test 3 == length(BridgeStan.param_constrain(model2, a; include_gq = true)) + @test 3 == length(BridgeStan.param_constrain(model2, a; include_gq = true, rng=rng)) @test 4 == length( - BridgeStan.param_constrain(model2, a; include_tp = true, include_gq = true), + BridgeStan.param_constrain(model2, a; include_tp = true, include_gq = true, rng=rng), ) + # reproducibility + @test isapprox(BridgeStan.param_constrain(model2, a; include_gq = true, seed=1234), + BridgeStan.param_constrain(model2, a; include_gq = true, seed=1234)) + + # no seed or rng provided + @test_throws ArgumentError BridgeStan.param_constrain(model2, a; include_gq = true) + # exception handling model3 = load_test_model("throw_tp", false) y = rand(1) @@ -199,6 +207,7 @@ end model4, y; include_gq = true, + seed=1234 ) end From 2e942982df1413cea088d96249effd4c6cec634b Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Wed, 29 Mar 2023 12:44:32 -0400 Subject: [PATCH 08/20] R interface work --- R/R/bridgestan.R | 99 ++++++++++++++++++++++++++---- R/man/StanModel.Rd | 39 ++++++++++-- R/tests/testthat/test-bridgestan.R | 23 ++++++- docs/languages/r.md | 34 ++++++++-- python/bridgestan/model.py | 4 -- src/bridgestanR.cpp | 26 ++++++-- src/bridgestanR.h | 16 ++++- 7 files changed, 204 insertions(+), 37 deletions(-) diff --git a/R/R/bridgestan.R b/R/R/bridgestan.R index fd9ca3eb..e3d0c792 100755 --- a/R/R/bridgestan.R +++ b/R/R/bridgestan.R @@ -12,9 +12,8 @@ StanModel <- R6::R6Class("StanModel", #' @param lib A path to a compiled BridgeStan Shared Object file. #' @param data Either a JSON string literal or a path to a data file in JSON format ending in ".json". #' @param rng_seed Seed for the RNG in the model object. - #' @param chain_id Used to offset the RNG by a fixed amount. #' @return A new StanModel. - initialize = function(lib, data, rng_seed, chain_id) { + initialize = function(lib, data, rng_seed) { if (.Platform$OS.type == "windows"){ lib_old <- lib lib <- paste0(tools::file_path_sans_ext(lib), ".dll") @@ -33,7 +32,7 @@ StanModel <- R6::R6Class("StanModel", dyn.load(private$lib, PACKAGE = private$lib_name) ret <- .C("bs_construct_R", - as.character(data), as.integer(rng_seed), as.integer(chain_id), + as.character(data), as.integer(rng_seed), ptr_out = raw(8), err_msg = as.character(""), err_ptr = raw(8), @@ -140,16 +139,41 @@ StanModel <- R6::R6Class("StanModel", #' @param theta_unc The vector of unconstrained parameters. #' @param include_tp Whether to also output the transformed parameters of the model. #' @param include_gq Whether to also output the generated quantities of the model. + #' @param seed The seed for the random number generator. One of + #' `rng` or `seed` must be specified if `include_gq` is `True`. + #' @param rng The random number generator to use. One of + #' `rng` or `seed` must be specified if `include_gq` is `True`. #' @return The constrained parameters of the model. - param_constrain = function(theta_unc, include_tp = FALSE, include_gq = FALSE) { - vars <- .C("bs_param_constrain_R", as.raw(private$model), - as.logical(include_tp), as.logical(include_gq), as.double(theta_unc), - theta = double(self$param_num(include_tp = include_tp, include_gq = include_gq)), - return_code = as.integer(0), - err_msg = as.character(""), - err_ptr = raw(8), - PACKAGE = private$lib_name - ) + param_constrain = function(theta_unc, include_tp = FALSE, include_gq = FALSE, seed, rng) { + if (missing(seed) && missing(rng)){ + if (include_gq) { + stop("Either seed or rng must be specified if include_gq is True.") + } else { + seed = 0 + } + } + if (!missing(rng)){ + vars <- .C("bs_param_constrain_R", as.raw(private$model), + as.logical(include_tp), as.logical(include_gq), as.double(theta_unc), + theta = double(self$param_num(include_tp = include_tp, include_gq = include_gq)), + rng = as.raw(rng$rng), + return_code = as.integer(0), + err_msg = as.character(""), + err_ptr = raw(8), + PACKAGE = private$lib_name + ) + } else { + vars <- .C("bs_param_constrain_seed_R", as.raw(private$model), + as.logical(include_tp), as.logical(include_gq), as.double(theta_unc), + theta = double(self$param_num(include_tp = include_tp, include_gq = include_gq)), + seed = as.integer(seed), + return_code = as.integer(0), + err_msg = as.character(""), + err_ptr = raw(8), + PACKAGE = private$lib_name + ) + } + if (vars$return_code) { stop(handle_error(private$lib_name, vars$err_msg, vars$err_ptr, "param_constrain")) } @@ -262,6 +286,13 @@ StanModel <- R6::R6Class("StanModel", stop(handle_error(private$lib_name, vars$err_msg, vars$err_ptr, "log_density_hessian")) } list(val = vars$val, gradient = vars$gradient, hessian = matrix(vars$hess, nrow = dims, byrow = TRUE)) + }, + #' @description + #' Create a new persistent RNG object for use in `param_constrain()`. + #' @param seed The seed for the RNG. + #' @return A `StanRNG` object. + new_rng = function(seed) { + StanRNG$new(private$lib_name, seed) } ), private = list( @@ -287,3 +318,47 @@ handle_error <- function(lib_name, err_msg, err_ptr, function_name) { return(err_msg) } } + +#' @description +#' RNG object for use with `StanModel$param_constrain()` +#' @field rng The pointer to the RNG object. +#' @keywords internal +StanRNG <- R6::R6Class("StanRNG", + public = list( + #' @description + #' Create a StanRng + #' @param lib_name The name of the Stan dynamic library. + #' @param rng_seed The seed for the RNG. + #' @return A new StanRNG. + initialize = function(lib_name, rng_seed) { + private$lib_name <- lib_name + + vars <- .C("bs_construct_rng_R", + as.integer(rng_seed), + ptr_out = raw(8), + err_msg = as.character(""), + err_ptr = raw(8), + PACKAGE = private$lib_name + ) + + if (all(vars$ptr_out == 0)) { + stop(handle_error("construct_rng", vars$err_msg, vars$err_ptr, private$lib_name)) + } else { + self$rng <- vars$ptr_out + } + }, + rng = NA + ), + private = list( + lib = NA, + lib_name = NA, + finalize = function() { + .C("bs_destruct_rng_R", + as.raw(self$rng), + return_code = as.integer(0), + PACKAGE = private$lib_name + ) + } + ), + cloneable=FALSE +) diff --git a/R/man/StanModel.Rd b/R/man/StanModel.Rd index 19c4f85a..37b21038 100644 --- a/R/man/StanModel.Rd +++ b/R/man/StanModel.Rd @@ -30,6 +30,7 @@ as well as constraining and unconstraining transforms. \item \href{#method-StanModel-log_density}{\code{StanModel$log_density()}} \item \href{#method-StanModel-log_density_gradient}{\code{StanModel$log_density_gradient()}} \item \href{#method-StanModel-log_density_hessian}{\code{StanModel$log_density_hessian()}} +\item \href{#method-StanModel-new_rng}{\code{StanModel$new_rng()}} } } \if{html}{\out{
}} @@ -38,7 +39,7 @@ as well as constraining and unconstraining transforms. \subsection{Method \code{new()}}{ Create a Stan Model instance. \subsection{Usage}{ -\if{html}{\out{
}}\preformatted{StanModel$new(lib, data, rng_seed, chain_id)}\if{html}{\out{
}} +\if{html}{\out{
}}\preformatted{StanModel$new(lib, data, rng_seed)}\if{html}{\out{
}} } \subsection{Arguments}{ @@ -49,8 +50,6 @@ Create a Stan Model instance. \item{\code{data}}{Either a JSON string literal or a path to a data file in JSON format ending in ".json".} \item{\code{rng_seed}}{Seed for the RNG in the model object.} - -\item{\code{chain_id}}{Used to offset the RNG by a fixed amount.} } \if{html}{\out{}} } @@ -174,7 +173,13 @@ The number of parameters in the model. Returns a vector of constrained parameters given the unconstrained parameters. See also \code{StanModel$param_unconstrain()}, the inverse of this function. \subsection{Usage}{ -\if{html}{\out{
}}\preformatted{StanModel$param_constrain(theta_unc, include_tp = FALSE, include_gq = FALSE)}\if{html}{\out{
}} +\if{html}{\out{
}}\preformatted{StanModel$param_constrain( + theta_unc, + include_tp = FALSE, + include_gq = FALSE, + seed, + rng +)}\if{html}{\out{
}} } \subsection{Arguments}{ @@ -185,6 +190,12 @@ See also \code{StanModel$param_unconstrain()}, the inverse of this function. \item{\code{include_tp}}{Whether to also output the transformed parameters of the model.} \item{\code{include_gq}}{Whether to also output the generated quantities of the model.} + +\item{\code{seed}}{The seed for the random number generator. One of +\code{rng} or \code{seed} must be specified if \code{include_gq} is \code{True}.} + +\item{\code{rng}}{The random number generator to use. One of +\code{rng} or \code{seed} must be specified if \code{include_gq} is \code{True}.} } \if{html}{\out{}} } @@ -314,4 +325,24 @@ See also \code{StanModel$param_unconstrain()}, the inverse of this function. List containing entries \code{val} (the log density), \code{gradient} (the gradient), and \code{hessian} (the Hessian). } } +\if{html}{\out{
}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-StanModel-new_rng}{}}} +\subsection{Method \code{new_rng()}}{ +Create a new persistent RNG object for use in \code{param_constrain()}. +\subsection{Usage}{ +\if{html}{\out{
}}\preformatted{StanModel$new_rng(seed)}\if{html}{\out{
}} +} + +\subsection{Arguments}{ +\if{html}{\out{
}} +\describe{ +\item{\code{seed}}{The seed for the RNG.} +} +\if{html}{\out{
}} +} +\subsection{Returns}{ +A \code{StanRNG} object. +} +} } diff --git a/R/tests/testthat/test-bridgestan.R b/R/tests/testthat/test-bridgestan.R index 2bf6215a..e1b9da08 100644 --- a/R/tests/testthat/test-bridgestan.R +++ b/R/tests/testthat/test-bridgestan.R @@ -6,12 +6,12 @@ load_model <- function(name, include_data=TRUE) { } else { data = "" } - model <- StanModel$new(file.path(base,paste0("/test_models/", name, "/",name,"_model.so")), data, 1234, 0) + model <- StanModel$new(file.path(base,paste0("/test_models/", name, "/",name,"_model.so")), data, 1234) return(model) } test_that("missing data throws error", { - expect_error(StanModel$new(file.path(base,paste0("/test_models/simple/simple_model.so")), "", 1234, 0)) + expect_error(load_model("simple",include_data=FALSE)) }) simple <- load_model("simple") @@ -106,6 +106,23 @@ test_that("param_unconstrain works for a nontrivial case", { expect_equal(c, a) }) +test_that("param_constrain handles rng arguments", { + full <- load_model("full", include_data=FALSE) + expect_equal(1, length(full$param_constrain(c(1.2)))) + expect_equal(2, length(full$param_constrain(c(1.2), include_tp=TRUE))) + rng <- full$new_rng(123) + expect_equal(3, length(full$param_constrain(c(1.2), include_gq=TRUE, rng=rng))) + expect_equal(4, length(full$param_constrain(c(1.2), include_tp=TRUE, include_gq=TRUE, rng=rng))) + + # check reproducibility + expect_equal(full$param_constrain(c(1.2), include_gq=TRUE, seed=1234), + full$param_constrain(c(1.2), include_gq=TRUE, seed=1234)) + + # require at least one present + expect_error(full$param_constrain(c(1.2), include_gq=TRUE), "seed or rng must be specified") +}) + + test_that("constructor propagates errors", { expect_error(load_model("throw_data",include_data=FALSE), "find this text: datafails") }) @@ -125,5 +142,5 @@ test_that("param_constrain propagates errors", { m2 <- load_model("throw_gq",include_data=FALSE) m2$param_constrain(c(1.2)) # no error - expect_error(m2$param_constrain(c(1.2), include_gq=TRUE), "find this text: gqfails") + expect_error(m2$param_constrain(c(1.2), include_gq=TRUE, seed=123), "find this text: gqfails") }) diff --git a/docs/languages/r.md b/docs/languages/r.md index b82dc911..df9dd167 100644 --- a/docs/languages/r.md +++ b/docs/languages/r.md @@ -62,7 +62,7 @@ Create a Stan Model instance. _Usage_ ```R -StanModel$new(lib, data, rng_seed, chain_id) +StanModel$new(lib, data, rng_seed) ``` @@ -74,9 +74,6 @@ _Arguments_ - `rng_seed` Seed for the RNG in the model object. - - `chain_id` Used to offset the RNG by a fixed amount. - - _Returns_ A new StanModel. @@ -220,7 +217,7 @@ of this function. _Usage_ ```R -StanModel$param_constrain(theta_unc, include_tp = FALSE, include_gq = FALSE) +StanModel$param_constrain(theta_unc, include_tp = FALSE, include_gq = FALSE, seed, rng) ``` @@ -234,6 +231,12 @@ _Arguments_ - `include_gq` Whether to also output the generated quantities of the model. + - `seed` Seed for the RNG in the model object. One of + `rng` or `seed` must be specified if `include_gq` is `True` + + - `rng` StanRNG to use in the model object. One of + `rng` or `seed` must be specified if `include_gq` is `True` + _Returns_ @@ -380,3 +383,24 @@ _Returns_ List containing entries `val` (the log density), `gradient` (the gradient), and `hessian` (the Hessian). + + +**Method** `new_rng()`: + +Create a new persistent RNG object for use in `param_constrain()`. + + +_Usage_ + +```R +StanModel$new_rng(seed) +``` + + +_Arguments_ + + - `seed` The seed for the RNG. + +_Returns_ + + A `StanRNG` object. diff --git a/python/bridgestan/model.py b/python/bridgestan/model.py index eed9dc3c..7691c97f 100644 --- a/python/bridgestan/model.py +++ b/python/bridgestan/model.py @@ -28,8 +28,6 @@ class StanModel: :param model_data: Either a JSON string literal or a path to a data file in JSON format ending in ``.json``. :param seed: A pseudo random number generator seed. - :param chain_id: A unique identifier for concurrent chains of - pseudorandom numbers. :raises FileNotFoundError or PermissionError: If ``model_lib`` is not readable or ``model_data`` is specified and not a path to a readable file. :raises RuntimeError: If there is an error instantiating the Stan model. @@ -224,8 +222,6 @@ def from_stan_file( threading for the compiled model. If the same flags are defined in ``make/local``, the versions passed here will take precedent. :param seed: A pseudo random number generator seed. - :param chain_id: A unique identifier for concurrent chains of - pseudorandom numbers. :raises FileNotFoundError or PermissionError: If `stan_file` does not exist or is not readable. :raises ValueError: If BridgeStan cannot be located. diff --git a/src/bridgestanR.cpp b/src/bridgestanR.cpp index a0c43f2c..6da40ef3 100644 --- a/src/bridgestanR.cpp +++ b/src/bridgestanR.cpp @@ -1,8 +1,8 @@ #include "bridgestanR.h" #include "bridgestan.h" -void bs_construct_R(char** data, int* rng, bs_model** ptr_out, - char** err_msg, void** err_ptr) { +void bs_construct_R(char** data, int* rng, bs_model** ptr_out, char** err_msg, + void** err_ptr) { *ptr_out = bs_construct(*data, *rng, err_msg); *err_ptr = static_cast(*err_msg); } @@ -36,12 +36,18 @@ void bs_param_unc_num_R(bs_model** model, int* num_out) { *num_out = bs_param_unc_num(*model); } void bs_param_constrain_R(bs_model** model, int* include_tp, int* include_gq, - const double* theta_unc, double* theta, + const double* theta_unc, double* theta, bs_rng** rng, int* return_code, char** err_msg, void** err_ptr) { - /* TODO(bmw): accept seed as argument, consider persistent version */ - bs_rng rng(0); *return_code = bs_param_constrain(*model, *include_tp, *include_gq, theta_unc, - theta, &rng, err_msg); + theta, *rng, err_msg); + *err_ptr = static_cast(*err_msg); +} +void bs_param_constrain_seed_R(bs_model** model, int* include_tp, + int* include_gq, const double* theta_unc, + double* theta, int* seed, int* return_code, + char** err_msg, void** err_ptr) { + *return_code = bs_param_constrain_seed(*model, *include_tp, *include_gq, + theta_unc, theta, *seed, err_msg); *err_ptr = static_cast(*err_msg); } void bs_param_unconstrain_R(bs_model** model, const double* theta, @@ -79,3 +85,11 @@ void bs_log_density_hessian_R(bs_model** model, int* propto, int* jacobian, val, grad, hess, err_msg); *err_ptr = static_cast(*err_msg); } +void bs_construct_rng_R(int* seed, bs_rng** ptr_out, char** err_msg, + void** err_ptr) { + *ptr_out = bs_construct_rng(*seed, err_msg); + *err_ptr = static_cast(*err_msg); +} +void bs_destruct_rng_R(bs_rng** rng, int* return_code) { + *return_code = bs_destruct_rng(*rng, nullptr); +} diff --git a/src/bridgestanR.h b/src/bridgestanR.h index f3a04d53..e0c2d4a2 100644 --- a/src/bridgestanR.h +++ b/src/bridgestanR.h @@ -12,8 +12,8 @@ typedef int bool; // Shim to convert to R interface requirement of void with pointer args // All calls directly delegated to versions without _R suffix -void bs_construct_R(char** data, int* rng, bs_model** ptr_out, - char** err_msg, void** err_ptr); +void bs_construct_R(char** data, int* rng, bs_model** ptr_out, char** err_msg, + void** err_ptr); void bs_version_R(int* major, int* minor, int* patch); @@ -40,9 +40,14 @@ void bs_param_num_R(bs_model** model, int* include_tp, int* include_gq, void bs_param_unc_num_R(bs_model** model, int* num_out); void bs_param_constrain_R(bs_model** model, int* include_tp, int* include_gq, - const double* theta_unc, double* theta, + const double* theta_unc, double* theta, bs_rng** rng, int* return_code, char** err_msg, void** err_ptr); +void bs_param_constrain_seed_R(bs_model** model, int* include_tp, + int* include_gq, const double* theta_unc, + double* theta, int* seed, int* return_code, + char** err_msg, void** err_ptr); + void bs_param_unconstrain_R(bs_model** model, const double* theta, double* theta_unc, int* return_code, char** err_msg, void** err_ptr); @@ -65,6 +70,11 @@ void bs_log_density_hessian_R(bs_model** model, int* propto, int* jacobian, double* grad, double* hess, int* return_code, char** err_msg, void** err_ptr); +void bs_construct_rng_R(int* seed, bs_rng** ptr_out, char** err_msg, + void** err_ptr); + +void bs_destruct_rng_R(bs_rng** rng, int* return_code); + #ifdef __cplusplus } #endif From 155671eb6b4559524d95d3adabf40b1ff769eec6 Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Wed, 29 Mar 2023 12:58:59 -0400 Subject: [PATCH 09/20] Documentation --- R/R/bridgestan.R | 21 ++++++++-------- R/man/StanModel.Rd | 46 +++++++++++++++++------------------ R/man/StanRNG.Rd | 50 ++++++++++++++++++++++++++++++++++++++ docs/languages/julia.md | 26 ++++++++++---------- docs/languages/r.md | 45 +++++++++++++++++----------------- julia/src/model.jl | 2 ++ python/bridgestan/model.py | 11 +++++---- src/bridgestan.h | 3 +++ 8 files changed, 131 insertions(+), 73 deletions(-) create mode 100644 R/man/StanRNG.Rd diff --git a/R/R/bridgestan.R b/R/R/bridgestan.R index e3d0c792..c13a1fd0 100755 --- a/R/R/bridgestan.R +++ b/R/R/bridgestan.R @@ -141,8 +141,8 @@ StanModel <- R6::R6Class("StanModel", #' @param include_gq Whether to also output the generated quantities of the model. #' @param seed The seed for the random number generator. One of #' `rng` or `seed` must be specified if `include_gq` is `True`. - #' @param rng The random number generator to use. One of - #' `rng` or `seed` must be specified if `include_gq` is `True`. + #' @param rng The random number generator to use. See `StanModel$new_rng()`. + #' One of `rng` or `seed` must be specified if `include_gq` is `True`. #' @return The constrained parameters of the model. param_constrain = function(theta_unc, include_tp = FALSE, include_gq = FALSE, seed, rng) { if (missing(seed) && missing(rng)){ @@ -180,6 +180,13 @@ StanModel <- R6::R6Class("StanModel", vars$theta }, #' @description + #' Create a new persistent PRNG object for use in `param_constrain()`. + #' @param seed The seed for the PRNG. + #' @return A `StanRNG` object. + new_rng = function(seed) { + StanRNG$new(private$lib_name, seed) + }, + #' @description #' Returns a vector of unconstrained parameters give the constrained parameters. #' #' It is assumed that these will be in the same order as internally represented by @@ -286,13 +293,6 @@ StanModel <- R6::R6Class("StanModel", stop(handle_error(private$lib_name, vars$err_msg, vars$err_ptr, "log_density_hessian")) } list(val = vars$val, gradient = vars$gradient, hessian = matrix(vars$hess, nrow = dims, byrow = TRUE)) - }, - #' @description - #' Create a new persistent RNG object for use in `param_constrain()`. - #' @param seed The seed for the RNG. - #' @return A `StanRNG` object. - new_rng = function(seed) { - StanRNG$new(private$lib_name, seed) } ), private = list( @@ -319,7 +319,8 @@ handle_error <- function(lib_name, err_msg, err_ptr, function_name) { } } -#' @description +#' StanRNG +#' #' RNG object for use with `StanModel$param_constrain()` #' @field rng The pointer to the RNG object. #' @keywords internal diff --git a/R/man/StanModel.Rd b/R/man/StanModel.Rd index 37b21038..466916cf 100644 --- a/R/man/StanModel.Rd +++ b/R/man/StanModel.Rd @@ -25,12 +25,12 @@ as well as constraining and unconstraining transforms. \item \href{#method-StanModel-param_num}{\code{StanModel$param_num()}} \item \href{#method-StanModel-param_unc_num}{\code{StanModel$param_unc_num()}} \item \href{#method-StanModel-param_constrain}{\code{StanModel$param_constrain()}} +\item \href{#method-StanModel-new_rng}{\code{StanModel$new_rng()}} \item \href{#method-StanModel-param_unconstrain}{\code{StanModel$param_unconstrain()}} \item \href{#method-StanModel-param_unconstrain_json}{\code{StanModel$param_unconstrain_json()}} \item \href{#method-StanModel-log_density}{\code{StanModel$log_density()}} \item \href{#method-StanModel-log_density_gradient}{\code{StanModel$log_density_gradient()}} \item \href{#method-StanModel-log_density_hessian}{\code{StanModel$log_density_hessian()}} -\item \href{#method-StanModel-new_rng}{\code{StanModel$new_rng()}} } } \if{html}{\out{
}} @@ -194,8 +194,8 @@ See also \code{StanModel$param_unconstrain()}, the inverse of this function. \item{\code{seed}}{The seed for the random number generator. One of \code{rng} or \code{seed} must be specified if \code{include_gq} is \code{True}.} -\item{\code{rng}}{The random number generator to use. One of -\code{rng} or \code{seed} must be specified if \code{include_gq} is \code{True}.} +\item{\code{rng}}{The random number generator to use. See \code{StanModel$new_rng()}. +One of \code{rng} or \code{seed} must be specified if \code{include_gq} is \code{True}.} } \if{html}{\out{}} } @@ -204,6 +204,26 @@ The constrained parameters of the model. } } \if{html}{\out{
}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-StanModel-new_rng}{}}} +\subsection{Method \code{new_rng()}}{ +Create a new persistent PRNG object for use in \code{param_constrain()}. +\subsection{Usage}{ +\if{html}{\out{
}}\preformatted{StanModel$new_rng(seed)}\if{html}{\out{
}} +} + +\subsection{Arguments}{ +\if{html}{\out{
}} +\describe{ +\item{\code{seed}}{The seed for the PRNG.} +} +\if{html}{\out{
}} +} +\subsection{Returns}{ +A \code{StanRNG} object. +} +} +\if{html}{\out{
}} \if{html}{\out{}} \if{latex}{\out{\hypertarget{method-StanModel-param_unconstrain}{}}} \subsection{Method \code{param_unconstrain()}}{ @@ -325,24 +345,4 @@ See also \code{StanModel$param_unconstrain()}, the inverse of this function. List containing entries \code{val} (the log density), \code{gradient} (the gradient), and \code{hessian} (the Hessian). } } -\if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-StanModel-new_rng}{}}} -\subsection{Method \code{new_rng()}}{ -Create a new persistent RNG object for use in \code{param_constrain()}. -\subsection{Usage}{ -\if{html}{\out{
}}\preformatted{StanModel$new_rng(seed)}\if{html}{\out{
}} -} - -\subsection{Arguments}{ -\if{html}{\out{
}} -\describe{ -\item{\code{seed}}{The seed for the RNG.} -} -\if{html}{\out{
}} -} -\subsection{Returns}{ -A \code{StanRNG} object. -} -} } diff --git a/R/man/StanRNG.Rd b/R/man/StanRNG.Rd new file mode 100644 index 00000000..d0219d38 --- /dev/null +++ b/R/man/StanRNG.Rd @@ -0,0 +1,50 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/bridgestan.R +\name{StanRNG} +\alias{StanRNG} +\title{StanRNG} +\description{ +StanRNG + +StanRNG +} +\details{ +RNG object for use with \code{StanModel$param_constrain()} +} +\keyword{internal} +\section{Public fields}{ +\if{html}{\out{
}} +\describe{ +\item{\code{rng}}{The pointer to the RNG object.} +} +\if{html}{\out{
}} +} +\section{Methods}{ +\subsection{Public methods}{ +\itemize{ +\item \href{#method-StanRNG-new}{\code{StanRNG$new()}} +} +} +\if{html}{\out{
}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-StanRNG-new}{}}} +\subsection{Method \code{new()}}{ +Create a StanRng +\subsection{Usage}{ +\if{html}{\out{
}}\preformatted{StanRNG$new(lib_name, rng_seed)}\if{html}{\out{
}} +} + +\subsection{Arguments}{ +\if{html}{\out{
}} +\describe{ +\item{\code{lib_name}}{The name of the Stan dynamic library.} + +\item{\code{rng_seed}}{The seed for the RNG.} +} +\if{html}{\out{
}} +} +\subsection{Returns}{ +A new StanRNG. +} +} +} diff --git a/docs/languages/julia.md b/docs/languages/julia.md index 96a85412..71af5d95 100644 --- a/docs/languages/julia.md +++ b/docs/languages/julia.md @@ -152,7 +152,7 @@ Return the log density of the specified unconstrained parameters. This calculation drops constant terms that do not depend on the parameters if `propto` is `true` and includes change of variables terms for constrained parameters if `jacobian` is `true`. -source
+source
# **`BridgeStan.log_density_gradient`** — *Function*. @@ -170,7 +170,7 @@ This calculation drops constant terms that do not depend on the parameters if `p This allocates new memory for the gradient output each call. See `log_density_gradient!` for a version which allows re-using existing memory. -source
+source
# **`BridgeStan.log_density_hessian`** — *Function*. @@ -188,7 +188,7 @@ This calculation drops constant terms that do not depend on the parameters if `p This allocates new memory for the gradient and Hessian output each call. See `log_density_gradient!` for a version which allows re-using existing memory. -source
+source
# **`BridgeStan.param_constrain`** — *Function*. @@ -201,14 +201,14 @@ param_constrain(sm, theta_unc, out; include_tp=false, include_gq=false, seed=not Returns a vector constrained parameters given unconstrained parameters. Additionally (if `include_tp` and `include_gq` are set, respectively) returns transformed parameters and generated quantities. -If `include_gq` is set, then either `seed` or `rng` must be provided. +If `include_gq` is set, then either `seed` or `rng` must be provided. See `StanRNG` for details on how to construct persistent RNGs. This allocates new memory for the output each call. See `param_constrain!` for a version which allows re-using existing memory. This is the inverse of `param_unconstrain`. -source
+source
# **`BridgeStan.param_unconstrain`** — *Function*. @@ -228,7 +228,7 @@ This allocates new memory for the output each call. See `param_unconstrain!` for This is the inverse of `param_constrain`. -source
+source
# **`BridgeStan.param_unconstrain_json`** — *Function*. @@ -246,7 +246,7 @@ The JSON is expected to be in the [JSON Format for CmdStan](https://mc-stan.org/ This allocates new memory for the output each call. See `param_unconstrain_json!` for a version which allows re-using existing memory. -source
+source
# **`BridgeStan.name`** — *Function*. @@ -360,7 +360,7 @@ This calculation drops constant terms that do not depend on the parameters if `p The gradient is stored in the vector `out`, and a reference is returned. See `log_density_gradient` for a version which allocates fresh memory. -source
+source
# **`BridgeStan.log_density_hessian!`** — *Function*. @@ -378,7 +378,7 @@ This calculation drops constant terms that do not depend on the parameters if `p The gradient is stored in the vector `out_grad` and the Hessian is stored in `out_hess` and references are returned. See `log_density_hessian` for a version which allocates fresh memory. -source
+source
# **`BridgeStan.param_constrain!`** — *Function*. @@ -391,14 +391,14 @@ param_constrain!(sm, theta_unc, out; include_tp=false, include_gq=false, seed=no Returns a vector constrained parameters given unconstrained parameters. Additionally (if `include_tp` and `include_gq` are set, respectively) returns transformed parameters and generated quantities. -If `include_gq` is set, then either `seed` or `rng` must be provided. +If `include_gq` is set, then either `seed` or `rng` must be provided. See `StanRNG` for details on how to construct persistent RNGs. The result is stored in the vector `out`, and a reference is returned. See `param_constrain` for a version which allocates fresh memory. This is the inverse of `param_unconstrain!`. -source
+source
# **`BridgeStan.param_unconstrain!`** — *Function*. @@ -418,7 +418,7 @@ The result is stored in the vector `out`, and a reference is returned. See `para This is the inverse of `param_constrain!`. -source
+source
# **`BridgeStan.param_unconstrain_json!`** — *Function*. @@ -436,7 +436,7 @@ The JSON is expected to be in the [JSON Format for CmdStan](https://mc-stan.org/ The result is stored in the vector `out`, and a reference is returned. See `param_unconstrain_json` for a version which allocates fresh memory. -source
+source
# **`BridgeStan.StanRNG`** — *Type*. diff --git a/docs/languages/r.md b/docs/languages/r.md index df9dd167..3d39b854 100644 --- a/docs/languages/r.md +++ b/docs/languages/r.md @@ -234,8 +234,8 @@ _Arguments_ - `seed` Seed for the RNG in the model object. One of `rng` or `seed` must be specified if `include_gq` is `True` - - `rng` StanRNG to use in the model object. One of - `rng` or `seed` must be specified if `include_gq` is `True` + - `rng` StanRNG to use in the model object. See `StanModel$new_rng()`. + One of `rng` or `seed` must be specified if `include_gq` is `True` _Returns_ @@ -243,6 +243,27 @@ _Returns_ The constrained parameters of the model. +**Method** `new_rng()`: + +Create a new persistent PRNG object for use in `param_constrain()`. + + +_Usage_ + +```R +StanModel$new_rng(seed) +``` + + +_Arguments_ + + - `seed` The seed for the PRNG. + +_Returns_ + + A `StanRNG` object. + + **Method** `param_unconstrain()`: Returns a vector of unconstrained parameters given the constrained parameters. @@ -384,23 +405,3 @@ _Returns_ List containing entries `val` (the log density), `gradient` (the gradient), and `hessian` (the Hessian). - -**Method** `new_rng()`: - -Create a new persistent RNG object for use in `param_constrain()`. - - -_Usage_ - -```R -StanModel$new_rng(seed) -``` - - -_Arguments_ - - - `seed` The seed for the RNG. - -_Returns_ - - A `StanRNG` object. diff --git a/julia/src/model.jl b/julia/src/model.jl index fb6008dd..47aa22e6 100644 --- a/julia/src/model.jl +++ b/julia/src/model.jl @@ -268,6 +268,7 @@ Additionally (if `include_tp` and `include_gq` are set, respectively) returns transformed parameters and generated quantities. If `include_gq` is set, then either `seed` or `rng` must be provided. +See `StanRNG` for details on how to construct persistent RNGs. The result is stored in the vector `out`, and a reference is returned. See `param_constrain` for a version which allocates fresh memory. @@ -347,6 +348,7 @@ Additionally (if `include_tp` and `include_gq` are set, respectively) returns transformed parameters and generated quantities. If `include_gq` is set, then either `seed` or `rng` must be provided. +See `StanRNG` for details on how to construct persistent RNGs. This allocates new memory for the output each call. See `param_constrain!` for a version which allows diff --git a/python/bridgestan/model.py b/python/bridgestan/model.py index 7691c97f..856bf629 100644 --- a/python/bridgestan/model.py +++ b/python/bridgestan/model.py @@ -328,8 +328,8 @@ def param_constrain( include_tp: bool = False, include_gq: bool = False, out: Optional[FloatArray] = None, - rng: Optional["StanRNG"] = None, seed: Optional[int] = None, + rng: Optional["StanRNG"] = None, ) -> FloatArray: """ Return the constrained parameters derived from the specified @@ -345,12 +345,13 @@ def param_constrain( provided, it must have shape `(D, )`, where `D` is the number of constrained parameters. If not provided or `None`, a freshly allocated array is returned. - :param rng: A ``StanRNG`` object to use for generating random - numbers. One of ``rng`` or ``seed`` must be specified if - ``include_gq`` is ``True``. :param seed: A pseudo random number generator seed. One of ``rng`` or ``seed`` must be specified if ``include_gq`` is ``True``. + :param rng: A ``StanRNG`` object to use for generating random + numbers, see :meth:`~StanModel.new_rng`. + One of ``rng`` or ``seed`` must be specified if ``include_gq`` + is ``True``. :return: The constrained parameter array. :raises ValueError: If ``out`` is specified and is not the same shape as the return. @@ -403,7 +404,7 @@ def param_constrain( def new_rng(self, seed: int) -> "StanRNG": """ - Return a new PRNG for the model. + Return a new PRNG for use in :meth:`~StanModel.param_constrain`. :param seed: The seed for the PRNG. :return: A new PRNG for the model. diff --git a/src/bridgestan.h b/src/bridgestan.h index 166f20c3..13c455c1 100644 --- a/src/bridgestan.h +++ b/src/bridgestan.h @@ -163,6 +163,9 @@ int bs_param_constrain(const bs_model* m, bool include_tp, bool include_gq, * in the Stan program, with multivariate parameters given in * last-index-major order. * + * This version accepts a seed which is used to create a fresh PRNG + * which lives only for the duration of this call. + * * @param[in] mr pointer to model and RNG structure * @param[in] include_tp `true` to include transformed parameters * @param[in] include_gq `true` to include generated quantities From 64c28ad6d94f8d26810efdb636b39dd74ad00523 Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Wed, 29 Mar 2023 14:25:10 -0400 Subject: [PATCH 10/20] Update c-example --- c-example/example.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/c-example/example.c b/c-example/example.c index b30c161f..a955bce5 100644 --- a/c-example/example.c +++ b/c-example/example.c @@ -14,7 +14,7 @@ int main(int argc, char** argv) { // this could potentially error, and we may get information back about why. char* err; - bs_model_rng* model = bs_construct(data, 123, 0, &err); + bs_model* model = bs_construct(data, 123, &err); if (!model) { if (err) { printf("Error: %s", err); From 7304b04282518b2bddbab7dfc02226221d16c67b Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Mon, 3 Apr 2023 11:56:39 -0400 Subject: [PATCH 11/20] Rework to use model seed (by default) and chain offset --- R/R/bridgestan.R | 45 ++++++++++++++++---------- R/man/StanModel.Rd | 19 +++++++---- R/man/StanRNG.Rd | 6 ++-- R/tests/testthat/test-bridgestan.R | 8 ++--- docs/languages/julia.md | 50 +++++++++++++++-------------- docs/languages/r.md | 15 +++++---- julia/src/model.jl | 50 +++++++++++++++++------------ julia/test/model_tests.jl | 8 ++--- python/bridgestan/model.py | 51 ++++++++++++++++-------------- python/test/test_stanmodel.py | 8 ++--- src/bridgestan.cpp | 10 +++--- src/bridgestan.h | 26 ++++++++------- src/bridgestanR.cpp | 12 +++---- src/bridgestanR.h | 6 ++-- src/model_rng.cpp | 8 +++++ src/model_rng.hpp | 13 +++++++- 16 files changed, 197 insertions(+), 138 deletions(-) diff --git a/R/R/bridgestan.R b/R/R/bridgestan.R index c13a1fd0..7fa64a8d 100755 --- a/R/R/bridgestan.R +++ b/R/R/bridgestan.R @@ -11,15 +11,16 @@ StanModel <- R6::R6Class("StanModel", #' Create a Stan Model instance. #' @param lib A path to a compiled BridgeStan Shared Object file. #' @param data Either a JSON string literal or a path to a data file in JSON format ending in ".json". - #' @param rng_seed Seed for the RNG in the model object. + #' @param seed Seed for the RNG in the model object. #' @return A new StanModel. - initialize = function(lib, data, rng_seed) { + initialize = function(lib, data, seed) { if (.Platform$OS.type == "windows"){ lib_old <- lib lib <- paste0(tools::file_path_sans_ext(lib), ".dll") file.copy(from=lib_old, to=lib) } + private$seed <- seed private$lib <- tools::file_path_as_absolute(lib) private$lib_name <- tools::file_path_sans_ext(basename(lib)) if (is.loaded("construct_R", PACKAGE = private$lib_name)) { @@ -32,7 +33,7 @@ StanModel <- R6::R6Class("StanModel", dyn.load(private$lib, PACKAGE = private$lib_name) ret <- .C("bs_construct_R", - as.character(data), as.integer(rng_seed), + as.character(data), as.integer(seed), ptr_out = raw(8), err_msg = as.character(""), err_ptr = raw(8), @@ -139,17 +140,18 @@ StanModel <- R6::R6Class("StanModel", #' @param theta_unc The vector of unconstrained parameters. #' @param include_tp Whether to also output the transformed parameters of the model. #' @param include_gq Whether to also output the generated quantities of the model. - #' @param seed The seed for the random number generator. One of - #' `rng` or `seed` must be specified if `include_gq` is `True`. + #' @param chain_id A chain ID used to offset a PRNG seeded with the model's seed + #' which should be unique between calls. One of `rng` or `seed` must be specified + #' if `include_gq` is `True`. #' @param rng The random number generator to use. See `StanModel$new_rng()`. #' One of `rng` or `seed` must be specified if `include_gq` is `True`. #' @return The constrained parameters of the model. - param_constrain = function(theta_unc, include_tp = FALSE, include_gq = FALSE, seed, rng) { - if (missing(seed) && missing(rng)){ + param_constrain = function(theta_unc, include_tp = FALSE, include_gq = FALSE, chain_id, rng) { + if (missing(chain_id) && missing(rng)){ if (include_gq) { - stop("Either seed or rng must be specified if include_gq is True.") + stop("Either chain_id or rng must be specified if include_gq is True.") } else { - seed = 0 + chain_id <- 0 } } if (!missing(rng)){ @@ -163,10 +165,10 @@ StanModel <- R6::R6Class("StanModel", PACKAGE = private$lib_name ) } else { - vars <- .C("bs_param_constrain_seed_R", as.raw(private$model), + vars <- .C("bs_param_constrain_id_R", as.raw(private$model), as.logical(include_tp), as.logical(include_gq), as.double(theta_unc), theta = double(self$param_num(include_tp = include_tp, include_gq = include_gq)), - seed = as.integer(seed), + chain_id = as.integer(chain_id), return_code = as.integer(0), err_msg = as.character(""), err_ptr = raw(8), @@ -181,10 +183,16 @@ StanModel <- R6::R6Class("StanModel", }, #' @description #' Create a new persistent PRNG object for use in `param_constrain()`. - #' @param seed The seed for the PRNG. + #' @param chain_id Identifier for a sequence in the RNG. This should be made + #' a distinct number for each PRNG created with the same seed + #' (for example, 1:N for N PRNGS). + #' @param seed The seed for the PRNG. If this is not specified, the model's seed is used. #' @return A `StanRNG` object. - new_rng = function(seed) { - StanRNG$new(private$lib_name, seed) + new_rng = function(chain_id, seed) { + if (missing(seed)){ + seed <- private$seed + } + StanRNG$new(private$lib_name, seed, chain_id) }, #' @description #' Returns a vector of unconstrained parameters give the constrained parameters. @@ -299,6 +307,7 @@ StanModel <- R6::R6Class("StanModel", lib = NA, lib_name = NA, model = NA, + seed = NA, finalize = function() { .C("bs_destruct_R", as.raw(private$model), @@ -329,13 +338,15 @@ StanRNG <- R6::R6Class("StanRNG", #' @description #' Create a StanRng #' @param lib_name The name of the Stan dynamic library. - #' @param rng_seed The seed for the RNG. + #' @param seed The seed for the RNG. + #' @param chain_id The chain ID for the RNG. #' @return A new StanRNG. - initialize = function(lib_name, rng_seed) { + initialize = function(lib_name, seed, chain_id) { private$lib_name <- lib_name vars <- .C("bs_construct_rng_R", - as.integer(rng_seed), + as.integer(seed), + as.integer(chain_id), ptr_out = raw(8), err_msg = as.character(""), err_ptr = raw(8), diff --git a/R/man/StanModel.Rd b/R/man/StanModel.Rd index 466916cf..89894be9 100644 --- a/R/man/StanModel.Rd +++ b/R/man/StanModel.Rd @@ -39,7 +39,7 @@ as well as constraining and unconstraining transforms. \subsection{Method \code{new()}}{ Create a Stan Model instance. \subsection{Usage}{ -\if{html}{\out{
}}\preformatted{StanModel$new(lib, data, rng_seed)}\if{html}{\out{
}} +\if{html}{\out{
}}\preformatted{StanModel$new(lib, data, seed)}\if{html}{\out{
}} } \subsection{Arguments}{ @@ -49,7 +49,7 @@ Create a Stan Model instance. \item{\code{data}}{Either a JSON string literal or a path to a data file in JSON format ending in ".json".} -\item{\code{rng_seed}}{Seed for the RNG in the model object.} +\item{\code{seed}}{Seed for the RNG in the model object.} } \if{html}{\out{}} } @@ -177,7 +177,7 @@ See also \code{StanModel$param_unconstrain()}, the inverse of this function. theta_unc, include_tp = FALSE, include_gq = FALSE, - seed, + chain_id, rng )}\if{html}{\out{}} } @@ -191,8 +191,9 @@ See also \code{StanModel$param_unconstrain()}, the inverse of this function. \item{\code{include_gq}}{Whether to also output the generated quantities of the model.} -\item{\code{seed}}{The seed for the random number generator. One of -\code{rng} or \code{seed} must be specified if \code{include_gq} is \code{True}.} +\item{\code{chain_id}}{A chain ID used to offset a PRNG seeded with the model's seed +which should be unique between calls. One of \code{rng} or \code{seed} must be specified +if \code{include_gq} is \code{True}.} \item{\code{rng}}{The random number generator to use. See \code{StanModel$new_rng()}. One of \code{rng} or \code{seed} must be specified if \code{include_gq} is \code{True}.} @@ -209,13 +210,17 @@ The constrained parameters of the model. \subsection{Method \code{new_rng()}}{ Create a new persistent PRNG object for use in \code{param_constrain()}. \subsection{Usage}{ -\if{html}{\out{
}}\preformatted{StanModel$new_rng(seed)}\if{html}{\out{
}} +\if{html}{\out{
}}\preformatted{StanModel$new_rng(chain_id, seed)}\if{html}{\out{
}} } \subsection{Arguments}{ \if{html}{\out{
}} \describe{ -\item{\code{seed}}{The seed for the PRNG.} +\item{\code{chain_id}}{Identifier for a sequence in the RNG. This should be made +a distinct number for each PRNG created with the same seed +(for example, 1:N for N PRNGS).} + +\item{\code{seed}}{The seed for the PRNG. If this is not specified, the model's seed is used.} } \if{html}{\out{
}} } diff --git a/R/man/StanRNG.Rd b/R/man/StanRNG.Rd index d0219d38..f0af6940 100644 --- a/R/man/StanRNG.Rd +++ b/R/man/StanRNG.Rd @@ -31,7 +31,7 @@ RNG object for use with \code{StanModel$param_constrain()} \subsection{Method \code{new()}}{ Create a StanRng \subsection{Usage}{ -\if{html}{\out{
}}\preformatted{StanRNG$new(lib_name, rng_seed)}\if{html}{\out{
}} +\if{html}{\out{
}}\preformatted{StanRNG$new(lib_name, seed, chain_id)}\if{html}{\out{
}} } \subsection{Arguments}{ @@ -39,7 +39,9 @@ Create a StanRng \describe{ \item{\code{lib_name}}{The name of the Stan dynamic library.} -\item{\code{rng_seed}}{The seed for the RNG.} +\item{\code{seed}}{The seed for the RNG.} + +\item{\code{chain_id}}{The chain ID for the RNG.} } \if{html}{\out{}} } diff --git a/R/tests/testthat/test-bridgestan.R b/R/tests/testthat/test-bridgestan.R index e1b9da08..97cfa683 100644 --- a/R/tests/testthat/test-bridgestan.R +++ b/R/tests/testthat/test-bridgestan.R @@ -115,11 +115,11 @@ test_that("param_constrain handles rng arguments", { expect_equal(4, length(full$param_constrain(c(1.2), include_tp=TRUE, include_gq=TRUE, rng=rng))) # check reproducibility - expect_equal(full$param_constrain(c(1.2), include_gq=TRUE, seed=1234), - full$param_constrain(c(1.2), include_gq=TRUE, seed=1234)) + expect_equal(full$param_constrain(c(1.2), include_gq=TRUE, chain_id=3), + full$param_constrain(c(1.2), include_gq=TRUE, chain_id=3)) # require at least one present - expect_error(full$param_constrain(c(1.2), include_gq=TRUE), "seed or rng must be specified") + expect_error(full$param_constrain(c(1.2), include_gq=TRUE), "chain_id or rng must be specified") }) @@ -142,5 +142,5 @@ test_that("param_constrain propagates errors", { m2 <- load_model("throw_gq",include_data=FALSE) m2$param_constrain(c(1.2)) # no error - expect_error(m2$param_constrain(c(1.2), include_gq=TRUE, seed=123), "find this text: gqfails") + expect_error(m2$param_constrain(c(1.2), include_gq=TRUE, chain_id=1), "find this text: gqfails") }) diff --git a/docs/languages/julia.md b/docs/languages/julia.md index 71af5d95..e31b2b29 100644 --- a/docs/languages/julia.md +++ b/docs/languages/julia.md @@ -152,7 +152,7 @@ Return the log density of the specified unconstrained parameters. This calculation drops constant terms that do not depend on the parameters if `propto` is `true` and includes change of variables terms for constrained parameters if `jacobian` is `true`. -source
+source
# **`BridgeStan.log_density_gradient`** — *Function*. @@ -170,7 +170,7 @@ This calculation drops constant terms that do not depend on the parameters if `p This allocates new memory for the gradient output each call. See `log_density_gradient!` for a version which allows re-using existing memory. -source
+source
# **`BridgeStan.log_density_hessian`** — *Function*. @@ -188,7 +188,7 @@ This calculation drops constant terms that do not depend on the parameters if `p This allocates new memory for the gradient and Hessian output each call. See `log_density_gradient!` for a version which allows re-using existing memory. -source
+source
# **`BridgeStan.param_constrain`** — *Function*. @@ -196,19 +196,19 @@ This allocates new memory for the gradient and Hessian output each call. See `lo ```julia -param_constrain(sm, theta_unc, out; include_tp=false, include_gq=false, seed=nothing, rng=nothing) +param_constrain(sm, theta_unc, out; include_tp=false, include_gq=false, chain_id=nothing, rng=nothing) ``` Returns a vector constrained parameters given unconstrained parameters. Additionally (if `include_tp` and `include_gq` are set, respectively) returns transformed parameters and generated quantities. -If `include_gq` is set, then either `seed` or `rng` must be provided. See `StanRNG` for details on how to construct persistent RNGs. +If `include_gq` is set, then either `chain_id` or `rng` must be provided. `chain_id` specifies an offset in a PRNG seeded with the model's seed which should be unique between calls. See `StanRNG` for details on how to construct persistent RNGs. This allocates new memory for the output each call. See `param_constrain!` for a version which allows re-using existing memory. This is the inverse of `param_unconstrain`. -source
+source
# **`BridgeStan.param_unconstrain`** — *Function*. @@ -228,7 +228,7 @@ This allocates new memory for the output each call. See `param_unconstrain!` for This is the inverse of `param_constrain`. -source
+source
# **`BridgeStan.param_unconstrain_json`** — *Function*. @@ -246,7 +246,7 @@ The JSON is expected to be in the [JSON Format for CmdStan](https://mc-stan.org/ This allocates new memory for the output each call. See `param_unconstrain_json!` for a version which allows re-using existing memory. -source
+source
# **`BridgeStan.name`** — *Function*. @@ -260,7 +260,7 @@ name(sm) Return the name of the model `sm` -source
+source
# **`BridgeStan.model_info`** — *Function*. @@ -276,7 +276,7 @@ Return information about the model `sm`. This includes the Stan version and important compiler flags. -source
+source
# **`BridgeStan.param_num`** — *Function*. @@ -292,7 +292,7 @@ Return the number of (constrained) parameters in the model. This is the total of all the sizes of items declared in the `parameters` block of the model. If `include_tp` or `include_gq` are true, items declared in the `transformed parameters` and `generate quantities` blocks are included, respectively. -source
+source
# **`BridgeStan.param_unc_num`** — *Function*. @@ -308,7 +308,7 @@ Return the number of unconstrained parameters in the model. This function is mainly different from `param_num` when variables are declared with constraints. For example, `simplex[5]` has a constrained size of 5, but an unconstrained size of 4. -source
+source
# **`BridgeStan.param_names`** — *Function*. @@ -326,7 +326,7 @@ For containers, indexes are separated by periods (.). For example, the scalar `a` has indexed name `"a"`, the vector entry `a[1]` has indexed name `"a.1"` and the matrix entry `a[2, 3]` has indexed names `"a.2.3"`. Parameter order of the output is column major and more generally last-index major for containers. -source
+source
# **`BridgeStan.param_unc_names`** — *Function*. @@ -342,7 +342,7 @@ Return the indexed names of the unconstrained parameters. For example, a scalar unconstrained parameter `b` has indexed name `b` and a vector entry `b[3]` has indexed name `b.3`. -source
+source
# **`BridgeStan.log_density_gradient!`** — *Function*. @@ -360,7 +360,7 @@ This calculation drops constant terms that do not depend on the parameters if `p The gradient is stored in the vector `out`, and a reference is returned. See `log_density_gradient` for a version which allocates fresh memory. -source
+source
# **`BridgeStan.log_density_hessian!`** — *Function*. @@ -378,7 +378,7 @@ This calculation drops constant terms that do not depend on the parameters if `p The gradient is stored in the vector `out_grad` and the Hessian is stored in `out_hess` and references are returned. See `log_density_hessian` for a version which allocates fresh memory. -source
+source
# **`BridgeStan.param_constrain!`** — *Function*. @@ -386,19 +386,19 @@ The gradient is stored in the vector `out_grad` and the Hessian is stored in `ou ```julia -param_constrain!(sm, theta_unc, out; include_tp=false, include_gq=false, seed=nothing, rng=nothing) +param_constrain!(sm, theta_unc, out; include_tp=false, include_gq=false, chain_id=nothing, rng=nothing) ``` Returns a vector constrained parameters given unconstrained parameters. Additionally (if `include_tp` and `include_gq` are set, respectively) returns transformed parameters and generated quantities. -If `include_gq` is set, then either `seed` or `rng` must be provided. See `StanRNG` for details on how to construct persistent RNGs. +If `include_gq` is set, then either `chain_id` or `rng` must be provided. `chain_id` specifies an offset in a PRNG seeded with the model's seed which should be unique between calls. See `StanRNG` for details on how to construct persistent RNGs. The result is stored in the vector `out`, and a reference is returned. See `param_constrain` for a version which allocates fresh memory. This is the inverse of `param_unconstrain!`. -source
+source
# **`BridgeStan.param_unconstrain!`** — *Function*. @@ -418,7 +418,7 @@ The result is stored in the vector `out`, and a reference is returned. See `para This is the inverse of `param_constrain!`. -source
+source
# **`BridgeStan.param_unconstrain_json!`** — *Function*. @@ -436,7 +436,7 @@ The JSON is expected to be in the [JSON Format for CmdStan](https://mc-stan.org/ The result is stored in the vector `out`, and a reference is returned. See `param_unconstrain_json` for a version which allocates fresh memory. -source
+source
# **`BridgeStan.StanRNG`** — *Type*. @@ -444,15 +444,17 @@ The result is stored in the vector `out`, and a reference is returned. See `para ```julia -StanRNG(sm::StanModel, seed) +StanRNG(sm::StanModel, chain_id; seed=nothing) ``` -Construct a StanRNG instance from a `StanModel` instance and a seed. This can be used in the `param_constrain` and `param_constrain!` methods when using the generated quantities block. +Construct a StanRNG instance from a `StanModel` instance and a chain ID. This ID serves as an offset in the PRNG stream and should be unique for each PRNG created with the same seed. A seed can be supplied to use in place of the model's seed, which is used by default. + +This can be used in the `param_constrain` and `param_constrain!` methods when using the generated quantities block. This object is not thread-safe, one should be created per thread. -source
+source
diff --git a/docs/languages/r.md b/docs/languages/r.md index 3d39b854..20c5199e 100644 --- a/docs/languages/r.md +++ b/docs/languages/r.md @@ -217,7 +217,7 @@ of this function. _Usage_ ```R -StanModel$param_constrain(theta_unc, include_tp = FALSE, include_gq = FALSE, seed, rng) +StanModel$param_constrain(theta_unc, include_tp = FALSE, include_gq = FALSE, chain_id, rng) ``` @@ -231,11 +231,11 @@ _Arguments_ - `include_gq` Whether to also output the generated quantities of the model. - - `seed` Seed for the RNG in the model object. One of - `rng` or `seed` must be specified if `include_gq` is `True` + - `chain_id` A chain ID used to offset a PRNG seeded with the model's seed + which should be unique between calls. One of `rng` or `chain_id` must be specified if `include_gq` is `True` - `rng` StanRNG to use in the model object. See `StanModel$new_rng()`. - One of `rng` or `seed` must be specified if `include_gq` is `True` + One of `rng` or `chain_id` must be specified if `include_gq` is `True` _Returns_ @@ -251,13 +251,16 @@ Create a new persistent PRNG object for use in `param_constrain()`. _Usage_ ```R -StanModel$new_rng(seed) +StanModel$new_rng(chain_id, seed) ``` _Arguments_ - - `seed` The seed for the PRNG. + - `chain_id` Identifier for a sequence in the RNG. This should be made a distinct number for each PRNG created with the sane seed (for example, 1:N for N PRNGS). + + - `seed` The seed for the PRNG. If this is not specified, the model's seed is used. + _Returns_ diff --git a/julia/src/model.jl b/julia/src/model.jl index 47aa22e6..ef9d716c 100644 --- a/julia/src/model.jl +++ b/julia/src/model.jl @@ -89,9 +89,13 @@ mutable struct StanModel end """ - StanRNG(sm::StanModel, seed) + StanRNG(sm::StanModel, chain_id; seed=nothing) + +Construct a StanRNG instance from a `StanModel` instance and a chain ID. +This ID serves as an offset in the PRNG stream and should be unique for each +PRNG created with the same seed. +A seed can be supplied to use in place of the model's seed, which is used by default. -Construct a StanRNG instance from a `StanModel` instance and a seed. This can be used in the `param_constrain` and `param_constrain!` methods when using the generated quantities block. @@ -101,24 +105,26 @@ mutable struct StanRNG lib::Ptr{Nothing} rng::Ptr{StanRNGStruct} seed::UInt32 + chain_id::UInt32 - function StanRNG(sm::StanModel, seed) - seed = convert(UInt32, seed) - + function StanRNG(sm::StanModel, chain_id; seed=nothing) + seed = convert(UInt32, if isnothing(seed) sm.seed else seed end) + chain_id = convert(UInt32, chain_id) err = Ref{Cstring}() rng = ccall( Libc.Libdl.dlsym(sm.lib, "bs_construct_rng"), Ptr{StanModelStruct}, - (UInt32, Ref{Cstring}), + (UInt32, UInt32, Ref{Cstring}), seed, + chain_id, err, ) if rng == C_NULL error(_handle_error(sm.lib, err, "bs_construct_rng")) end - stanrng = new(sm.lib, rng, seed) + stanrng = new(sm.lib, rng, seed, chain_id) function f(stanrng) ccall( @@ -261,13 +267,15 @@ function param_unc_names(sm::StanModel) end """ - param_constrain!(sm, theta_unc, out; include_tp=false, include_gq=false, seed=nothing, rng=nothing) + param_constrain!(sm, theta_unc, out; include_tp=false, include_gq=false, chain_id=nothing, rng=nothing) Returns a vector constrained parameters given unconstrained parameters. Additionally (if `include_tp` and `include_gq` are set, respectively) returns transformed parameters and generated quantities. -If `include_gq` is set, then either `seed` or `rng` must be provided. +If `include_gq` is set, then either `chain_id` or `rng` must be provided. +`chain_id` specifies an offset in a PRNG seeded with the model's +seed which should be unique between calls. See `StanRNG` for details on how to construct persistent RNGs. The result is stored in the vector `out`, and a reference is returned. See @@ -281,7 +289,7 @@ function param_constrain!( out::Vector{Float64}; include_tp = false, include_gq = false, - seed::Union{Int, Nothing} = nothing, + chain_id::Union{Int, Nothing} = nothing, rng::Union{StanRNG, Nothing} = nothing, ) dims = param_num(sm; include_tp = include_tp, include_gq = include_gq) @@ -291,15 +299,15 @@ function param_constrain!( ) end - if seed === nothing && rng === nothing + if chain_id === nothing && rng === nothing if include_gq throw( ArgumentError( - "Must provide either a seed or an RNG when including generated quantities", + "Must provide either a chain_id or an RNG when including generated quantities", ), ) else - seed = 0 + chain_id = 0 end end @@ -320,9 +328,9 @@ function param_constrain!( err, ) else - seed = convert(UInt32, seed) + chain_id = convert(UInt32, chain_id) rc = ccall( - Libc.Libdl.dlsym(sm.lib, "bs_param_constrain_seed"), + Libc.Libdl.dlsym(sm.lib, "bs_param_constrain_id"), Cint, (Ptr{StanModelStruct}, Cint, Cint, Ref{Cdouble}, Ref{Cdouble}, Cuint, Ref{Cstring}), sm.stanmodel, @@ -330,7 +338,7 @@ function param_constrain!( include_gq, theta_unc, out, - seed, + chain_id, err, ) end @@ -341,13 +349,15 @@ function param_constrain!( end """ - param_constrain(sm, theta_unc, out; include_tp=false, include_gq=false, seed=nothing, rng=nothing) + param_constrain(sm, theta_unc, out; include_tp=false, include_gq=false, chain_id=nothing, rng=nothing) Returns a vector constrained parameters given unconstrained parameters. Additionally (if `include_tp` and `include_gq` are set, respectively) returns transformed parameters and generated quantities. -If `include_gq` is set, then either `seed` or `rng` must be provided. +If `include_gq` is set, then either `chain_id` or `rng` must be provided. +`chain_id` specifies an offset in a PRNG seeded with the model's +seed which should be unique between calls. See `StanRNG` for details on how to construct persistent RNGs. This allocates new memory for the output each call. @@ -361,11 +371,11 @@ function param_constrain( theta_unc::Vector{Float64}; include_tp = false, include_gq = false, - seed::Union{Int, Nothing} = nothing, + chain_id::Union{Int, Nothing} = nothing, rng::Union{StanRNG, Nothing} = nothing, ) out = zeros(param_num(sm, include_tp = include_tp, include_gq = include_gq)) - param_constrain!(sm, theta_unc, out; include_tp = include_tp, include_gq = include_gq, seed=seed, rng=rng) + param_constrain!(sm, theta_unc, out; include_tp = include_tp, include_gq = include_gq, chain_id=chain_id, rng=rng) end """ diff --git a/julia/test/model_tests.jl b/julia/test/model_tests.jl index feed1660..48248ac0 100644 --- a/julia/test/model_tests.jl +++ b/julia/test/model_tests.jl @@ -176,7 +176,7 @@ end model2 = load_test_model("full", false) - rng = StanRNG(model2, 1234) + rng = StanRNG(model2, 2; seed=1234) @test 1 == length(BridgeStan.param_constrain(model2, a)) @test 2 == length(BridgeStan.param_constrain(model2, a; include_tp = true)) @test 3 == length(BridgeStan.param_constrain(model2, a; include_gq = true, rng=rng)) @@ -185,8 +185,8 @@ end ) # reproducibility - @test isapprox(BridgeStan.param_constrain(model2, a; include_gq = true, seed=1234), - BridgeStan.param_constrain(model2, a; include_gq = true, seed=1234)) + @test isapprox(BridgeStan.param_constrain(model2, a; include_gq = true, chain_id=3), + BridgeStan.param_constrain(model2, a; include_gq = true, chain_id=3)) # no seed or rng provided @test_throws ArgumentError BridgeStan.param_constrain(model2, a; include_gq = true) @@ -207,7 +207,7 @@ end model4, y; include_gq = true, - seed=1234 + chain_id=1, ) end diff --git a/python/bridgestan/model.py b/python/bridgestan/model.py index 856bf629..5a2d6034 100644 --- a/python/bridgestan/model.py +++ b/python/bridgestan/model.py @@ -47,8 +47,7 @@ def __init__( :param model_lib: A system path to compiled shared object. :param model_data: Either a JSON string literal or a system path to a data file in JSON format ending in ``.json``. - :param seed: A pseudo random number generator seed. This is only used - if the model has RNG usage in the ``transformed data`` block. + :param seed: A pseudo random number generator seed. :raises FileNotFoundError or PermissionError: If ``model_lib`` is not readable or ``model_data`` is specified and not a path to a readable file. :raises RuntimeError: If there is an error instantiating the @@ -127,9 +126,9 @@ def __init__( star_star_char, ] - self._param_constrain_seed = self.stanlib.bs_param_constrain_seed - self._param_constrain_seed.restype = ctypes.c_int - self._param_constrain_seed.argtypes = [ + self._param_constrain_id = self.stanlib.bs_param_constrain_id + self._param_constrain_id.restype = ctypes.c_int + self._param_constrain_id.argtypes = [ ctypes.c_void_p, ctypes.c_int, ctypes.c_int, @@ -239,7 +238,7 @@ def __del__(self) -> None: def __repr__(self) -> str: data = f"{self.data_path!r}, " if self.data_path else "" - return f"StanModel({self.lib_path!r}, {data}seed={self.seed})" + return f"StanModel({self.lib_path!r}, {data}, seed={self.seed})" def name(self) -> str: """ @@ -328,7 +327,7 @@ def param_constrain( include_tp: bool = False, include_gq: bool = False, out: Optional[FloatArray] = None, - seed: Optional[int] = None, + chain_id: Optional[int] = None, rng: Optional["StanRNG"] = None, ) -> FloatArray: """ @@ -345,28 +344,28 @@ def param_constrain( provided, it must have shape `(D, )`, where `D` is the number of constrained parameters. If not provided or `None`, a freshly allocated array is returned. - :param seed: A pseudo random number generator seed. One of - ``rng`` or ``seed`` must be specified if ``include_gq`` - is ``True``. + :param chain_id: A chain ID used to offset a PRNG seeded with the model's + seed which should be unique between calls. One of ``rng`` or + ``chain_id`` must be specified if ``include_gq`` is ``True``. :param rng: A ``StanRNG`` object to use for generating random numbers, see :meth:`~StanModel.new_rng`. - One of ``rng`` or ``seed`` must be specified if ``include_gq`` + One of ``rng`` or ``chain_id`` must be specified if ``include_gq`` is ``True``. :return: The constrained parameter array. :raises ValueError: If ``out`` is specified and is not the same shape as the return. - :raises ValueError: If neither ``rng`` nor ``seed`` is specified + :raises ValueError: If neither ``rng`` nor ``chain_id`` is specified and ``include_gq`` is ``True``. :raises RuntimeError: If the C++ Stan model throws an exception. """ - if seed is None and rng is None: + if chain_id is None and rng is None: if include_gq: raise ValueError( - "Error: must specify rng or seed when including generated quantities" + "Error: must specify rng or chain_id when including generated quantities" ) else: - # neither specified, but not doing gq, so use a fixed seed - seed = 0 + # neither specified, but not doing gq, so use a fixed chain id + chain_id = 0 dims = self.param_num(include_tp=include_tp, include_gq=include_gq) if out is None: @@ -388,13 +387,13 @@ def param_constrain( err, ) else: - rc = self._param_constrain_seed( + rc = self._param_constrain_id( self.model_rng, int(include_tp), int(include_gq), theta_unc, out, - seed, + chain_id, err, ) @@ -402,14 +401,18 @@ def param_constrain( raise self._handle_error(err.contents, "param_constrain") return out - def new_rng(self, seed: int) -> "StanRNG": + def new_rng(self, chain_id: int, *, seed=None) -> "StanRNG": """ Return a new PRNG for use in :meth:`~StanModel.param_constrain`. - :param seed: The seed for the PRNG. + :param chain_id: Identifier for a sequence in the RNG. This should be made + a distinct number for each PRNG created with the same seed + (for example, 1:N for N PRNGS). + :param seed: A seed for the PRNG. If not specified, the + model's seed is used. :return: A new PRNG for the model. """ - return StanRNG(self.stanlib, seed) + return StanRNG(self.stanlib, seed or self.seed, chain_id) def param_unconstrain( self, theta: FloatArray, *, out: Optional[FloatArray] = None @@ -625,13 +628,13 @@ def _handle_error(self, err: ctypes.c_char_p, method: str) -> Exception: class StanRNG: - def __init__(self, lib: ctypes.CDLL, seed: int) -> None: + def __init__(self, lib: ctypes.CDLL, seed: int, chain_id: int) -> None: self.stanlib = lib construct = self.stanlib.bs_construct_rng construct.restype = ctypes.c_void_p - construct.argtypes = [ctypes.c_uint, star_star_char] - self.ptr = construct(seed, None) + construct.argtypes = [ctypes.c_uint, ctypes.c_uint, star_star_char] + self.ptr = construct(seed, chain_id, None) if not self.ptr: raise RuntimeError("Failed to construct RNG.") diff --git a/python/test/test_stanmodel.py b/python/test/test_stanmodel.py index 32914d77..4d319ac1 100644 --- a/python/test/test_stanmodel.py +++ b/python/test/test_stanmodel.py @@ -208,7 +208,7 @@ def test_param_constrain(): full_so = str(STAN_FOLDER / "full" / "full_model.so") bridge2 = bs.StanModel(full_so) - rng = bridge2.new_rng(seed=1234) + rng = bridge2.new_rng(chain_id=1, seed=1234) np.testing.assert_equal(1, bridge2.param_constrain(a).size) np.testing.assert_equal(2, bridge2.param_constrain(a, include_tp=True).size) @@ -221,8 +221,8 @@ def test_param_constrain(): # reproducibility test np.testing.assert_equal( - bridge2.param_constrain(a, include_gq=True, seed=1234), - bridge2.param_constrain(a, include_gq=True, seed=1234), + bridge2.param_constrain(a, include_gq=True, chain_id=4), + bridge2.param_constrain(a, include_gq=True, chain_id=4), ) # test error if neither seed or rng is provided @@ -251,7 +251,7 @@ def test_param_constrain(): bridge3 = bs.StanModel(throw_gq_so) bridge3.param_constrain(y, include_gq=False) with pytest.raises(RuntimeError, match="find this text: gqfails"): - bridge3.param_constrain(y, include_gq=True, seed=123) + bridge3.param_constrain(y, include_gq=True, chain_id=1) def test_param_unconstrain(): diff --git a/src/bridgestan.cpp b/src/bridgestan.cpp index 9f880262..1bf200c7 100644 --- a/src/bridgestan.cpp +++ b/src/bridgestan.cpp @@ -78,11 +78,11 @@ int bs_param_constrain(const bs_model* m, bool include_tp, bool include_gq, return 1; } -int bs_param_constrain_seed(const bs_model* m, bool include_tp, +int bs_param_constrain_id(const bs_model* m, bool include_tp, bool include_gq, const double* theta_unc, - double* theta, unsigned int seed, + double* theta, unsigned int chain_id, char** error_msg) { - bs_rng rng(seed); + bs_rng rng(m->seed(), chain_id); return bs_param_constrain(m, include_tp, include_gq, theta_unc, theta, &rng, error_msg); } @@ -201,9 +201,9 @@ int bs_log_density_hessian(const bs_model* m, bool propto, bool jacobian, return -1; } -bs_rng* bs_construct_rng(unsigned int seed, char** error_msg) { +bs_rng* bs_construct_rng(unsigned int seed, unsigned int chain_id, char** error_msg) { try { - return new bs_rng(seed); + return new bs_rng(seed, chain_id); } catch (const std::exception& e) { if (error_msg) { std::stringstream error; diff --git a/src/bridgestan.h b/src/bridgestan.h index 13c455c1..be8a9582 100644 --- a/src/bridgestan.h +++ b/src/bridgestan.h @@ -28,7 +28,8 @@ extern int bs_patch_version; * @return pointer to constructed model or `nullptr` if construction * fails */ -bs_model* bs_construct(const char* data_file, unsigned int seed, char** error_msg); +bs_model* bs_construct(const char* data_file, unsigned int seed, + char** error_msg); /** * Destroy the model. @@ -163,25 +164,26 @@ int bs_param_constrain(const bs_model* m, bool include_tp, bool include_gq, * in the Stan program, with multivariate parameters given in * last-index-major order. * - * This version accepts a seed which is used to create a fresh PRNG - * which lives only for the duration of this call. + * This version accepts a chain_id which is used to create a PRNG + * offset from the model's seed which lives only for the duration + * of this call. * * @param[in] mr pointer to model and RNG structure * @param[in] include_tp `true` to include transformed parameters * @param[in] include_gq `true` to include generated quantities * @param[in] theta_unc sequence of unconstrained parameters * @param[out] theta sequence of constrained parameters - * @param[in] seed seed for pseudorandom number generator which will be created - * and destroyed during this call. See `bs_param_constrain` for an option with a - * persistent RNG. + * @param[in] chain_id offset for pseudorandom number generator which will be + * created and destroyed during this call. seeded with model seed. See + * `bs_param_constrain` for an option with a persistent RNG. * @param[out] error_msg a pointer to a string that will be allocated if there * is an error. This must later be freed by calling `bs_free_error_msg`. * @return code 0 if successful and code -1 if there is an exception * in the underlying Stan code */ -int bs_param_constrain_seed(const bs_model* mr, bool include_tp, - bool include_gq, const double* theta_unc, - double* theta, unsigned int seed, char** error_msg); +int bs_param_constrain_id(const bs_model* mr, bool include_tp, bool include_gq, + const double* theta_unc, double* theta, + unsigned int chain_id, char** error_msg); /** * Set the sequence of unconstrained parameters based on the @@ -294,15 +296,17 @@ int bs_log_density_hessian(const bs_model* m, bool propto, bool jacobian, double* hessian, char** error_msg); /** - * Construct an RNG object to be used in `bs_param_constrain`. + * Construct an PRNG object to be used in `bs_param_constrain`. * This object is not thread safe and should be constructed and * destructed for each thread. * * @param[in] seed seed for the RNG + * @param[in] chain_id identifier for the current sequence of PRNG draws * @param[out] error_msg a pointer to a string that will be allocated if there * is an error. This must later be freed by calling `bs_free_error_msg`. */ -bs_rng* bs_construct_rng(unsigned int seed, char** error_msg); +bs_rng* bs_construct_rng(unsigned int seed, unsigned int chain_id, + char** error_msg); /** * Destruct an RNG object. diff --git a/src/bridgestanR.cpp b/src/bridgestanR.cpp index 6da40ef3..b6e8ff42 100644 --- a/src/bridgestanR.cpp +++ b/src/bridgestanR.cpp @@ -42,12 +42,12 @@ void bs_param_constrain_R(bs_model** model, int* include_tp, int* include_gq, theta, *rng, err_msg); *err_ptr = static_cast(*err_msg); } -void bs_param_constrain_seed_R(bs_model** model, int* include_tp, +void bs_param_constrain_id_R(bs_model** model, int* include_tp, int* include_gq, const double* theta_unc, - double* theta, int* seed, int* return_code, + double* theta, int* chain_id, int* return_code, char** err_msg, void** err_ptr) { - *return_code = bs_param_constrain_seed(*model, *include_tp, *include_gq, - theta_unc, theta, *seed, err_msg); + *return_code = bs_param_constrain_id(*model, *include_tp, *include_gq, + theta_unc, theta, *chain_id, err_msg); *err_ptr = static_cast(*err_msg); } void bs_param_unconstrain_R(bs_model** model, const double* theta, @@ -85,9 +85,9 @@ void bs_log_density_hessian_R(bs_model** model, int* propto, int* jacobian, val, grad, hess, err_msg); *err_ptr = static_cast(*err_msg); } -void bs_construct_rng_R(int* seed, bs_rng** ptr_out, char** err_msg, +void bs_construct_rng_R(int* seed, int* chain_id, bs_rng** ptr_out, char** err_msg, void** err_ptr) { - *ptr_out = bs_construct_rng(*seed, err_msg); + *ptr_out = bs_construct_rng(*seed, *chain_id, err_msg); *err_ptr = static_cast(*err_msg); } void bs_destruct_rng_R(bs_rng** rng, int* return_code) { diff --git a/src/bridgestanR.h b/src/bridgestanR.h index e0c2d4a2..381df9fe 100644 --- a/src/bridgestanR.h +++ b/src/bridgestanR.h @@ -43,9 +43,9 @@ void bs_param_constrain_R(bs_model** model, int* include_tp, int* include_gq, const double* theta_unc, double* theta, bs_rng** rng, int* return_code, char** err_msg, void** err_ptr); -void bs_param_constrain_seed_R(bs_model** model, int* include_tp, +void bs_param_constrain_id_R(bs_model** model, int* include_tp, int* include_gq, const double* theta_unc, - double* theta, int* seed, int* return_code, + double* theta, int* chain_id, int* return_code, char** err_msg, void** err_ptr); void bs_param_unconstrain_R(bs_model** model, const double* theta, @@ -70,7 +70,7 @@ void bs_log_density_hessian_R(bs_model** model, int* propto, int* jacobian, double* grad, double* hess, int* return_code, char** err_msg, void** err_ptr); -void bs_construct_rng_R(int* seed, bs_rng** ptr_out, char** err_msg, +void bs_construct_rng_R(int* seed, int* chain_id, bs_rng** ptr_out, char** err_msg, void** err_ptr); void bs_destruct_rng_R(bs_rng** rng, int* return_code); diff --git a/src/model_rng.cpp b/src/model_rng.cpp index 48f20ec1..375f14ad 100644 --- a/src/model_rng.cpp +++ b/src/model_rng.cpp @@ -78,6 +78,8 @@ bs_model::bs_model(const char* data_file, unsigned int seed) { } } + seed_ = seed; + std::string model_name = model_->model_name(); const char* model_name_c = model_name.c_str(); name_ = strdup(model_name_c); @@ -167,6 +169,8 @@ bs_model::~bs_model() noexcept { const char* bs_model::name() const { return name_; } +unsigned int bs_model::seed() const { return seed_; } + const char* bs_model::model_info() const { return model_info_; } const char* bs_model::param_names(bool include_tp, bool include_gq) const { @@ -323,3 +327,7 @@ void bs_model::log_density_hessian(bool propto, bool jacobian, Eigen::VectorXd::Map(grad, N) = grad_vec; Eigen::MatrixXd::Map(hessian, N, N) = hess_mat; } + +bs_rng::bs_rng(unsigned int seed, unsigned int chain_id) { + rng_ = stan::services::util::create_rng(seed, chain_id + 1); +} diff --git a/src/model_rng.hpp b/src/model_rng.hpp index 6d1feb82..0509a98a 100644 --- a/src/model_rng.hpp +++ b/src/model_rng.hpp @@ -38,6 +38,14 @@ class bs_model { */ const char* name() const; + /** + * Return the pseudorandom number generator seed + * used during model construction. + * + * @return seed + */ + unsigned int seed() const; + /** * Return information about the compiled model. This class manages the * memory, so the returned string should not be freed. @@ -182,6 +190,9 @@ class bs_model { /** Stan model */ stan::model::model_base* model_; + /** RNG seed provided during model creation */ + unsigned int seed_; + /** name of the Stan model */ char* name_ = nullptr; @@ -230,7 +241,7 @@ class bs_model { */ class bs_rng { public: - bs_rng(unsigned int seed) : rng_(seed) { rng_.discard(1); } + bs_rng(unsigned int seed, unsigned int chain_id); boost::ecuyer1988 rng_; }; From 65bb3d7d6f2cfd6694770d786b0cb678bf15eca7 Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Tue, 4 Apr 2023 14:59:45 -0400 Subject: [PATCH 12/20] Expose seed AND chain in C API, keep higher-level APIs as prior commit --- R/R/bridgestan.R | 3 ++- julia/src/model.jl | 5 +++-- python/bridgestan/model.py | 10 ++++++---- src/bridgestan.cpp | 6 +++--- src/bridgestan.h | 4 ++-- src/bridgestanR.cpp | 8 ++++---- src/bridgestanR.h | 4 ++-- src/model_rng.cpp | 4 ---- src/model_rng.hpp | 11 ----------- 9 files changed, 22 insertions(+), 33 deletions(-) diff --git a/R/R/bridgestan.R b/R/R/bridgestan.R index 7fa64a8d..a8bd000a 100755 --- a/R/R/bridgestan.R +++ b/R/R/bridgestan.R @@ -165,9 +165,10 @@ StanModel <- R6::R6Class("StanModel", PACKAGE = private$lib_name ) } else { - vars <- .C("bs_param_constrain_id_R", as.raw(private$model), + vars <- .C("bs_param_constrain_seeded_R", as.raw(private$model), as.logical(include_tp), as.logical(include_gq), as.double(theta_unc), theta = double(self$param_num(include_tp = include_tp, include_gq = include_gq)), + seed = as.integer(private$seed), chain_id = as.integer(chain_id), return_code = as.integer(0), err_msg = as.character(""), diff --git a/julia/src/model.jl b/julia/src/model.jl index ef9d716c..eacb398c 100644 --- a/julia/src/model.jl +++ b/julia/src/model.jl @@ -330,14 +330,15 @@ function param_constrain!( else chain_id = convert(UInt32, chain_id) rc = ccall( - Libc.Libdl.dlsym(sm.lib, "bs_param_constrain_id"), + Libc.Libdl.dlsym(sm.lib, "bs_param_constrain_seeded"), Cint, - (Ptr{StanModelStruct}, Cint, Cint, Ref{Cdouble}, Ref{Cdouble}, Cuint, Ref{Cstring}), + (Ptr{StanModelStruct}, Cint, Cint, Ref{Cdouble}, Ref{Cdouble}, Cuint, Cuint, Ref{Cstring}), sm.stanmodel, include_tp, include_gq, theta_unc, out, + sm.seed, chain_id, err, ) diff --git a/python/bridgestan/model.py b/python/bridgestan/model.py index 5a2d6034..ab2659ed 100644 --- a/python/bridgestan/model.py +++ b/python/bridgestan/model.py @@ -126,15 +126,16 @@ def __init__( star_star_char, ] - self._param_constrain_id = self.stanlib.bs_param_constrain_id - self._param_constrain_id.restype = ctypes.c_int - self._param_constrain_id.argtypes = [ + self._param_constrain_seeded = self.stanlib.bs_param_constrain_seeded + self._param_constrain_seeded.restype = ctypes.c_int + self._param_constrain_seeded.argtypes = [ ctypes.c_void_p, ctypes.c_int, ctypes.c_int, double_array, double_array, ctypes.c_uint, + ctypes.c_uint, star_star_char, ] @@ -387,12 +388,13 @@ def param_constrain( err, ) else: - rc = self._param_constrain_id( + rc = self._param_constrain_seeded( self.model_rng, int(include_tp), int(include_gq), theta_unc, out, + self.seed, chain_id, err, ) diff --git a/src/bridgestan.cpp b/src/bridgestan.cpp index 1bf200c7..0d525626 100644 --- a/src/bridgestan.cpp +++ b/src/bridgestan.cpp @@ -78,11 +78,11 @@ int bs_param_constrain(const bs_model* m, bool include_tp, bool include_gq, return 1; } -int bs_param_constrain_id(const bs_model* m, bool include_tp, +int bs_param_constrain_seeded(const bs_model* m, bool include_tp, bool include_gq, const double* theta_unc, - double* theta, unsigned int chain_id, + double* theta, unsigned int seed, unsigned int chain_id, char** error_msg) { - bs_rng rng(m->seed(), chain_id); + bs_rng rng(seed, chain_id); return bs_param_constrain(m, include_tp, include_gq, theta_unc, theta, &rng, error_msg); } diff --git a/src/bridgestan.h b/src/bridgestan.h index be8a9582..08400b04 100644 --- a/src/bridgestan.h +++ b/src/bridgestan.h @@ -181,8 +181,8 @@ int bs_param_constrain(const bs_model* m, bool include_tp, bool include_gq, * @return code 0 if successful and code -1 if there is an exception * in the underlying Stan code */ -int bs_param_constrain_id(const bs_model* mr, bool include_tp, bool include_gq, - const double* theta_unc, double* theta, +int bs_param_constrain_seeded(const bs_model* mr, bool include_tp, bool include_gq, + const double* theta_unc, double* theta, unsigned int seed, unsigned int chain_id, char** error_msg); /** diff --git a/src/bridgestanR.cpp b/src/bridgestanR.cpp index b6e8ff42..69c2372b 100644 --- a/src/bridgestanR.cpp +++ b/src/bridgestanR.cpp @@ -42,12 +42,12 @@ void bs_param_constrain_R(bs_model** model, int* include_tp, int* include_gq, theta, *rng, err_msg); *err_ptr = static_cast(*err_msg); } -void bs_param_constrain_id_R(bs_model** model, int* include_tp, +void bs_param_constrain_seeded_R(bs_model** model, int* include_tp, int* include_gq, const double* theta_unc, - double* theta, int* chain_id, int* return_code, + double* theta, int* seed, int* chain_id, int* return_code, char** err_msg, void** err_ptr) { - *return_code = bs_param_constrain_id(*model, *include_tp, *include_gq, - theta_unc, theta, *chain_id, err_msg); + *return_code = bs_param_constrain_seeded(*model, *include_tp, *include_gq, + theta_unc, theta, *seed, *chain_id, err_msg); *err_ptr = static_cast(*err_msg); } void bs_param_unconstrain_R(bs_model** model, const double* theta, diff --git a/src/bridgestanR.h b/src/bridgestanR.h index 381df9fe..ec07826d 100644 --- a/src/bridgestanR.h +++ b/src/bridgestanR.h @@ -43,9 +43,9 @@ void bs_param_constrain_R(bs_model** model, int* include_tp, int* include_gq, const double* theta_unc, double* theta, bs_rng** rng, int* return_code, char** err_msg, void** err_ptr); -void bs_param_constrain_id_R(bs_model** model, int* include_tp, +void bs_param_constrain_seeded_R(bs_model** model, int* include_tp, int* include_gq, const double* theta_unc, - double* theta, int* chain_id, int* return_code, + double* theta, int* seed, int* chain_id, int* return_code, char** err_msg, void** err_ptr); void bs_param_unconstrain_R(bs_model** model, const double* theta, diff --git a/src/model_rng.cpp b/src/model_rng.cpp index 375f14ad..3ef28186 100644 --- a/src/model_rng.cpp +++ b/src/model_rng.cpp @@ -78,8 +78,6 @@ bs_model::bs_model(const char* data_file, unsigned int seed) { } } - seed_ = seed; - std::string model_name = model_->model_name(); const char* model_name_c = model_name.c_str(); name_ = strdup(model_name_c); @@ -169,8 +167,6 @@ bs_model::~bs_model() noexcept { const char* bs_model::name() const { return name_; } -unsigned int bs_model::seed() const { return seed_; } - const char* bs_model::model_info() const { return model_info_; } const char* bs_model::param_names(bool include_tp, bool include_gq) const { diff --git a/src/model_rng.hpp b/src/model_rng.hpp index 0509a98a..482cf2f0 100644 --- a/src/model_rng.hpp +++ b/src/model_rng.hpp @@ -38,14 +38,6 @@ class bs_model { */ const char* name() const; - /** - * Return the pseudorandom number generator seed - * used during model construction. - * - * @return seed - */ - unsigned int seed() const; - /** * Return information about the compiled model. This class manages the * memory, so the returned string should not be freed. @@ -190,9 +182,6 @@ class bs_model { /** Stan model */ stan::model::model_base* model_; - /** RNG seed provided during model creation */ - unsigned int seed_; - /** name of the Stan model */ char* name_ = nullptr; From ae1e15df6754245c428f671773ba1b371d02a9bb Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Tue, 4 Apr 2023 15:03:51 -0400 Subject: [PATCH 13/20] Formatting --- julia/src/model.jl | 104 ++++++++++++++++++++++------------ julia/test/model_tests.jl | 20 +++++-- python/bridgestan/__init__.py | 2 +- python/bridgestan/compile.py | 2 +- python/bridgestan/download.py | 6 +- python/setup.py | 2 +- python/test/conftest.py | 7 ++- python/test/test_download.py | 2 +- src/bridgestan.cpp | 9 +-- src/bridgestan.h | 10 ++-- src/bridgestanR.cpp | 16 +++--- src/bridgestanR.h | 11 ++-- 12 files changed, 118 insertions(+), 73 deletions(-) diff --git a/julia/src/model.jl b/julia/src/model.jl index eacb398c..fd6493b3 100644 --- a/julia/src/model.jl +++ b/julia/src/model.jl @@ -107,8 +107,12 @@ mutable struct StanRNG seed::UInt32 chain_id::UInt32 - function StanRNG(sm::StanModel, chain_id; seed=nothing) - seed = convert(UInt32, if isnothing(seed) sm.seed else seed end) + function StanRNG(sm::StanModel, chain_id; seed = nothing) + seed = convert(UInt32, if isnothing(seed) + sm.seed + else + seed + end) chain_id = convert(UInt32, chain_id) err = Ref{Cstring}() @@ -130,8 +134,9 @@ mutable struct StanRNG ccall( Libc.Libdl.dlsym(stanrng.lib, "bs_destruct_rng"), UInt32, - (Ptr{StanModelStruct},Ref{Cstring}), - stanrng.rng, C_NULL + (Ptr{StanModelStruct}, Ref{Cstring}), + stanrng.rng, + C_NULL, ) end @@ -289,8 +294,8 @@ function param_constrain!( out::Vector{Float64}; include_tp = false, include_gq = false, - chain_id::Union{Int, Nothing} = nothing, - rng::Union{StanRNG, Nothing} = nothing, + chain_id::Union{Int,Nothing} = nothing, + rng::Union{StanRNG,Nothing} = nothing, ) dims = param_num(sm; include_tp = include_tp, include_gq = include_gq) if length(out) != dims @@ -301,11 +306,11 @@ function param_constrain!( if chain_id === nothing && rng === nothing if include_gq - throw( - ArgumentError( - "Must provide either a chain_id or an RNG when including generated quantities", - ), - ) + throw( + ArgumentError( + "Must provide either a chain_id or an RNG when including generated quantities", + ), + ) else chain_id = 0 end @@ -316,32 +321,49 @@ function param_constrain!( if rng !== nothing rc = ccall( - Libc.Libdl.dlsym(sm.lib, "bs_param_constrain"), - Cint, - (Ptr{StanModelStruct}, Cint, Cint, Ref{Cdouble}, Ref{Cdouble}, Ptr{StanRNGStruct},Ref{Cstring}), - sm.stanmodel, - include_tp, - include_gq, - theta_unc, - out, - rng.rng, - err, - ) + Libc.Libdl.dlsym(sm.lib, "bs_param_constrain"), + Cint, + ( + Ptr{StanModelStruct}, + Cint, + Cint, + Ref{Cdouble}, + Ref{Cdouble}, + Ptr{StanRNGStruct}, + Ref{Cstring}, + ), + sm.stanmodel, + include_tp, + include_gq, + theta_unc, + out, + rng.rng, + err, + ) else chain_id = convert(UInt32, chain_id) rc = ccall( - Libc.Libdl.dlsym(sm.lib, "bs_param_constrain_seeded"), - Cint, - (Ptr{StanModelStruct}, Cint, Cint, Ref{Cdouble}, Ref{Cdouble}, Cuint, Cuint, Ref{Cstring}), - sm.stanmodel, - include_tp, - include_gq, - theta_unc, - out, - sm.seed, - chain_id, - err, - ) + Libc.Libdl.dlsym(sm.lib, "bs_param_constrain_seeded"), + Cint, + ( + Ptr{StanModelStruct}, + Cint, + Cint, + Ref{Cdouble}, + Ref{Cdouble}, + Cuint, + Cuint, + Ref{Cstring}, + ), + sm.stanmodel, + include_tp, + include_gq, + theta_unc, + out, + sm.seed, + chain_id, + err, + ) end if rc != 0 error(handle_error(sm.lib, err, "param_constrain")) @@ -372,11 +394,19 @@ function param_constrain( theta_unc::Vector{Float64}; include_tp = false, include_gq = false, - chain_id::Union{Int, Nothing} = nothing, - rng::Union{StanRNG, Nothing} = nothing, + chain_id::Union{Int,Nothing} = nothing, + rng::Union{StanRNG,Nothing} = nothing, ) out = zeros(param_num(sm, include_tp = include_tp, include_gq = include_gq)) - param_constrain!(sm, theta_unc, out; include_tp = include_tp, include_gq = include_gq, chain_id=chain_id, rng=rng) + param_constrain!( + sm, + theta_unc, + out; + include_tp = include_tp, + include_gq = include_gq, + chain_id = chain_id, + rng = rng, + ) end """ diff --git a/julia/test/model_tests.jl b/julia/test/model_tests.jl index 48248ac0..b786449c 100644 --- a/julia/test/model_tests.jl +++ b/julia/test/model_tests.jl @@ -176,17 +176,25 @@ end model2 = load_test_model("full", false) - rng = StanRNG(model2, 2; seed=1234) + rng = StanRNG(model2, 2; seed = 1234) @test 1 == length(BridgeStan.param_constrain(model2, a)) @test 2 == length(BridgeStan.param_constrain(model2, a; include_tp = true)) - @test 3 == length(BridgeStan.param_constrain(model2, a; include_gq = true, rng=rng)) + @test 3 == length(BridgeStan.param_constrain(model2, a; include_gq = true, rng = rng)) @test 4 == length( - BridgeStan.param_constrain(model2, a; include_tp = true, include_gq = true, rng=rng), + BridgeStan.param_constrain( + model2, + a; + include_tp = true, + include_gq = true, + rng = rng, + ), ) # reproducibility - @test isapprox(BridgeStan.param_constrain(model2, a; include_gq = true, chain_id=3), - BridgeStan.param_constrain(model2, a; include_gq = true, chain_id=3)) + @test isapprox( + BridgeStan.param_constrain(model2, a; include_gq = true, chain_id = 3), + BridgeStan.param_constrain(model2, a; include_gq = true, chain_id = 3), + ) # no seed or rng provided @test_throws ArgumentError BridgeStan.param_constrain(model2, a; include_gq = true) @@ -207,7 +215,7 @@ end model4, y; include_gq = true, - chain_id=1, + chain_id = 1, ) end diff --git a/python/bridgestan/__init__.py b/python/bridgestan/__init__.py index 7d15805f..24f4a9f5 100644 --- a/python/bridgestan/__init__.py +++ b/python/bridgestan/__init__.py @@ -1,5 +1,5 @@ +from .__version import __version__ from .compile import compile_model, set_bridgestan_path from .model import StanModel -from .__version import __version__ __all__ = ["StanModel", "set_bridgestan_path", "compile_model"] diff --git a/python/bridgestan/compile.py b/python/bridgestan/compile.py index a4bb1393..01c83457 100644 --- a/python/bridgestan/compile.py +++ b/python/bridgestan/compile.py @@ -4,9 +4,9 @@ from pathlib import Path from typing import List -from .util import validate_readable from .__version import __version__ from .download import CURRENT_BRIDGESTAN, get_bridgestan_src +from .util import validate_readable def verify_bridgestan_path(path: str) -> None: diff --git a/python/bridgestan/download.py b/python/bridgestan/download.py index 72c23d96..2fe8bb65 100644 --- a/python/bridgestan/download.py +++ b/python/bridgestan/download.py @@ -1,7 +1,7 @@ -from pathlib import Path -import urllib.request -import urllib.error import tarfile +import urllib.error +import urllib.request +from pathlib import Path from time import sleep from .__version import __version__ diff --git a/python/setup.py b/python/setup.py index 6441f7e3..f3f4621f 100644 --- a/python/setup.py +++ b/python/setup.py @@ -1,7 +1,7 @@ from ast import literal_eval -from setuptools import setup from pathlib import Path +from setuptools import setup VERSIONFILE = Path(__file__).parent / "bridgestan" / "__version.py" with open(VERSIONFILE, "rt") as f: diff --git a/python/test/conftest.py b/python/test/conftest.py index cd89932b..f6528c0b 100644 --- a/python/test/conftest.py +++ b/python/test/conftest.py @@ -17,6 +17,7 @@ types = ("default", "ad_hessian") + def pytest_addoption(parser): parser.addoption( "--run-type", @@ -29,8 +30,10 @@ def pytest_addoption(parser): def pytest_configure(config): for type in types[1:]: - config.addinivalue_line("markers", f"{type}: mark test as needing a seperate process in group {type}") - + config.addinivalue_line( + "markers", + f"{type}: mark test as needing a seperate process in group {type}", + ) def pytest_collection_modifyitems(config, items): diff --git a/python/test/test_download.py b/python/test/test_download.py index b646f4fd..81032674 100644 --- a/python/test/test_download.py +++ b/python/test/test_download.py @@ -1,7 +1,7 @@ import os +from unittest import mock import bridgestan as bs -from unittest import mock @mock.patch.dict(os.environ, {"BRIDGESTAN": ""}) diff --git a/src/bridgestan.cpp b/src/bridgestan.cpp index 0d525626..38600322 100644 --- a/src/bridgestan.cpp +++ b/src/bridgestan.cpp @@ -79,9 +79,9 @@ int bs_param_constrain(const bs_model* m, bool include_tp, bool include_gq, } int bs_param_constrain_seeded(const bs_model* m, bool include_tp, - bool include_gq, const double* theta_unc, - double* theta, unsigned int seed, unsigned int chain_id, - char** error_msg) { + bool include_gq, const double* theta_unc, + double* theta, unsigned int seed, + unsigned int chain_id, char** error_msg) { bs_rng rng(seed, chain_id); return bs_param_constrain(m, include_tp, include_gq, theta_unc, theta, &rng, error_msg); @@ -201,7 +201,8 @@ int bs_log_density_hessian(const bs_model* m, bool propto, bool jacobian, return -1; } -bs_rng* bs_construct_rng(unsigned int seed, unsigned int chain_id, char** error_msg) { +bs_rng* bs_construct_rng(unsigned int seed, unsigned int chain_id, + char** error_msg) { try { return new bs_rng(seed, chain_id); } catch (const std::exception& e) { diff --git a/src/bridgestan.h b/src/bridgestan.h index 08400b04..10219b04 100644 --- a/src/bridgestan.h +++ b/src/bridgestan.h @@ -87,8 +87,7 @@ const char* bs_model_info(const bs_model* m); * @param[in] include_gq `true` to include generated quantities * @return CSV-separated, indexed, parameter names */ -const char* bs_param_names(const bs_model* m, bool include_tp, - bool include_gq); +const char* bs_param_names(const bs_model* m, bool include_tp, bool include_gq); /** * Return a comma-separated sequence of unconstrained parameters. @@ -181,9 +180,10 @@ int bs_param_constrain(const bs_model* m, bool include_tp, bool include_gq, * @return code 0 if successful and code -1 if there is an exception * in the underlying Stan code */ -int bs_param_constrain_seeded(const bs_model* mr, bool include_tp, bool include_gq, - const double* theta_unc, double* theta, unsigned int seed, - unsigned int chain_id, char** error_msg); +int bs_param_constrain_seeded(const bs_model* mr, bool include_tp, + bool include_gq, const double* theta_unc, + double* theta, unsigned int seed, + unsigned int chain_id, char** error_msg); /** * Set the sequence of unconstrained parameters based on the diff --git a/src/bridgestanR.cpp b/src/bridgestanR.cpp index 69c2372b..2853b049 100644 --- a/src/bridgestanR.cpp +++ b/src/bridgestanR.cpp @@ -43,11 +43,13 @@ void bs_param_constrain_R(bs_model** model, int* include_tp, int* include_gq, *err_ptr = static_cast(*err_msg); } void bs_param_constrain_seeded_R(bs_model** model, int* include_tp, - int* include_gq, const double* theta_unc, - double* theta, int* seed, int* chain_id, int* return_code, - char** err_msg, void** err_ptr) { - *return_code = bs_param_constrain_seeded(*model, *include_tp, *include_gq, - theta_unc, theta, *seed, *chain_id, err_msg); + int* include_gq, const double* theta_unc, + double* theta, int* seed, int* chain_id, + int* return_code, char** err_msg, + void** err_ptr) { + *return_code + = bs_param_constrain_seeded(*model, *include_tp, *include_gq, theta_unc, + theta, *seed, *chain_id, err_msg); *err_ptr = static_cast(*err_msg); } void bs_param_unconstrain_R(bs_model** model, const double* theta, @@ -85,8 +87,8 @@ void bs_log_density_hessian_R(bs_model** model, int* propto, int* jacobian, val, grad, hess, err_msg); *err_ptr = static_cast(*err_msg); } -void bs_construct_rng_R(int* seed, int* chain_id, bs_rng** ptr_out, char** err_msg, - void** err_ptr) { +void bs_construct_rng_R(int* seed, int* chain_id, bs_rng** ptr_out, + char** err_msg, void** err_ptr) { *ptr_out = bs_construct_rng(*seed, *chain_id, err_msg); *err_ptr = static_cast(*err_msg); } diff --git a/src/bridgestanR.h b/src/bridgestanR.h index ec07826d..3d45b20c 100644 --- a/src/bridgestanR.h +++ b/src/bridgestanR.h @@ -44,9 +44,10 @@ void bs_param_constrain_R(bs_model** model, int* include_tp, int* include_gq, int* return_code, char** err_msg, void** err_ptr); void bs_param_constrain_seeded_R(bs_model** model, int* include_tp, - int* include_gq, const double* theta_unc, - double* theta, int* seed, int* chain_id, int* return_code, - char** err_msg, void** err_ptr); + int* include_gq, const double* theta_unc, + double* theta, int* seed, int* chain_id, + int* return_code, char** err_msg, + void** err_ptr); void bs_param_unconstrain_R(bs_model** model, const double* theta, double* theta_unc, int* return_code, char** err_msg, @@ -70,8 +71,8 @@ void bs_log_density_hessian_R(bs_model** model, int* propto, int* jacobian, double* grad, double* hess, int* return_code, char** err_msg, void** err_ptr); -void bs_construct_rng_R(int* seed, int* chain_id, bs_rng** ptr_out, char** err_msg, - void** err_ptr); +void bs_construct_rng_R(int* seed, int* chain_id, bs_rng** ptr_out, + char** err_msg, void** err_ptr); void bs_destruct_rng_R(bs_rng** rng, int* return_code); From c9b0b083596ffb9ffd6013479e6ee99da4c60a61 Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Thu, 6 Apr 2023 09:50:16 -0400 Subject: [PATCH 14/20] Update to reflect that destructor cannot fail --- R/R/bridgestan.R | 1 - julia/src/model.jl | 5 ++--- python/bridgestan/model.py | 6 +++--- src/bridgestan.cpp | 14 +------------- src/bridgestan.h | 4 +--- src/bridgestanR.cpp | 4 +--- src/bridgestanR.h | 2 +- 7 files changed, 9 insertions(+), 27 deletions(-) diff --git a/R/R/bridgestan.R b/R/R/bridgestan.R index a8bd000a..24036ae4 100755 --- a/R/R/bridgestan.R +++ b/R/R/bridgestan.R @@ -368,7 +368,6 @@ StanRNG <- R6::R6Class("StanRNG", finalize = function() { .C("bs_destruct_rng_R", as.raw(self$rng), - return_code = as.integer(0), PACKAGE = private$lib_name ) } diff --git a/julia/src/model.jl b/julia/src/model.jl index fd6493b3..1e6aea54 100644 --- a/julia/src/model.jl +++ b/julia/src/model.jl @@ -133,10 +133,9 @@ mutable struct StanRNG function f(stanrng) ccall( Libc.Libdl.dlsym(stanrng.lib, "bs_destruct_rng"), - UInt32, - (Ptr{StanModelStruct}, Ref{Cstring}), + Cvoid, + (Ptr{StanModelStruct}), stanrng.rng, - C_NULL, ) end diff --git a/python/bridgestan/model.py b/python/bridgestan/model.py index ab2659ed..2f5b93d2 100644 --- a/python/bridgestan/model.py +++ b/python/bridgestan/model.py @@ -642,12 +642,12 @@ def __init__(self, lib: ctypes.CDLL, seed: int, chain_id: int) -> None: raise RuntimeError("Failed to construct RNG.") self._destruct = self.stanlib.bs_destruct_rng - self._destruct.restype = ctypes.c_int - self._destruct.argtypes = [ctypes.c_void_p, star_star_char] + self._destruct.restype = None + self._destruct.argtypes = [ctypes.c_void_p] def __del__(self) -> None: """ Destroy the Stan model and free memory. """ if hasattr(self, "ptr") and hasattr(self, "_destruct"): - self._destruct(self.ptr, None) + self._destruct(self.ptr) diff --git a/src/bridgestan.cpp b/src/bridgestan.cpp index 38600322..9f30f2ad 100644 --- a/src/bridgestan.cpp +++ b/src/bridgestan.cpp @@ -223,16 +223,4 @@ bs_rng* bs_construct_rng(unsigned int seed, unsigned int chain_id, return nullptr; } -int bs_destruct_rng(bs_rng* mr, char** error_msg) { - try { - delete (mr); - return 0; - } catch (...) { - if (error_msg) { - std::stringstream error; - error << "destruct_rng() failed." << std::endl; - *error_msg = strdup(error.str().c_str()); - } - } - return -1; -} +void bs_destruct_rng(bs_rng* mr) { delete (mr); } diff --git a/src/bridgestan.h b/src/bridgestan.h index 10219b04..bfc20c0d 100644 --- a/src/bridgestan.h +++ b/src/bridgestan.h @@ -312,10 +312,8 @@ bs_rng* bs_construct_rng(unsigned int seed, unsigned int chain_id, * Destruct an RNG object. * * @param[in] rng pointer to RNG object - * @param[out] error_msg a pointer to a string that will be allocated if there - * is an error. This must later be freed by calling `bs_free_error_msg`. */ -int bs_destruct_rng(bs_rng* rng, char** error_msg); +void bs_destruct_rng(bs_rng* rng); #ifdef __cplusplus } diff --git a/src/bridgestanR.cpp b/src/bridgestanR.cpp index 2853b049..73ab69eb 100644 --- a/src/bridgestanR.cpp +++ b/src/bridgestanR.cpp @@ -92,6 +92,4 @@ void bs_construct_rng_R(int* seed, int* chain_id, bs_rng** ptr_out, *ptr_out = bs_construct_rng(*seed, *chain_id, err_msg); *err_ptr = static_cast(*err_msg); } -void bs_destruct_rng_R(bs_rng** rng, int* return_code) { - *return_code = bs_destruct_rng(*rng, nullptr); -} +void bs_destruct_rng_R(bs_rng** rng) { bs_destruct_rng(*rng); } diff --git a/src/bridgestanR.h b/src/bridgestanR.h index 3d45b20c..a226bbb3 100644 --- a/src/bridgestanR.h +++ b/src/bridgestanR.h @@ -74,7 +74,7 @@ void bs_log_density_hessian_R(bs_model** model, int* propto, int* jacobian, void bs_construct_rng_R(int* seed, int* chain_id, bs_rng** ptr_out, char** err_msg, void** err_ptr); -void bs_destruct_rng_R(bs_rng** rng, int* return_code); +void bs_destruct_rng_R(bs_rng** rng); #ifdef __cplusplus } From b084d6655e75c150327236186e2692e3894b6758 Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Thu, 6 Apr 2023 10:09:03 -0400 Subject: [PATCH 15/20] Typo fix --- julia/src/model.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/julia/src/model.jl b/julia/src/model.jl index 1e6aea54..7e26cd35 100644 --- a/julia/src/model.jl +++ b/julia/src/model.jl @@ -134,7 +134,7 @@ mutable struct StanRNG ccall( Libc.Libdl.dlsym(stanrng.lib, "bs_destruct_rng"), Cvoid, - (Ptr{StanModelStruct}), + (Ptr{StanModelStruct},), stanrng.rng, ) end From f856f9a7f4c4de406883d662ce53fb98cb1f9746 Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Tue, 11 Apr 2023 16:55:19 -0400 Subject: [PATCH 16/20] Simplify interface; remove concept of chain ID --- R/R/bridgestan.R | 67 ++++------- R/man/StanModel.Rd | 28 ++--- R/man/StanRNG.Rd | 4 +- R/man/handle_error.Rd | 2 +- R/tests/testthat/test-bridgestan.R | 8 +- docs/languages/julia.md | 48 ++++---- docs/languages/r.md | 21 ++-- julia/src/model.jl | 117 ++++++------------- julia/test/model_tests.jl | 19 +++- python/bridgestan/model.py | 175 ++++++++++++----------------- python/test/test_stanmodel.py | 8 +- src/bridgestan.cpp | 27 ++--- src/bridgestan.h | 69 ++++-------- src/bridgestanR.cpp | 14 +-- src/bridgestanR.h | 8 +- src/model_rng.cpp | 4 - src/model_rng.hpp | 8 +- 17 files changed, 240 insertions(+), 387 deletions(-) diff --git a/R/R/bridgestan.R b/R/R/bridgestan.R index 24036ae4..21ab62c4 100755 --- a/R/R/bridgestan.R +++ b/R/R/bridgestan.R @@ -11,7 +11,7 @@ StanModel <- R6::R6Class("StanModel", #' Create a Stan Model instance. #' @param lib A path to a compiled BridgeStan Shared Object file. #' @param data Either a JSON string literal or a path to a data file in JSON format ending in ".json". - #' @param seed Seed for the RNG in the model object. + #' @param seed Seed for the RNG used in constructing the model. #' @return A new StanModel. initialize = function(lib, data, seed) { if (.Platform$OS.type == "windows"){ @@ -140,42 +140,26 @@ StanModel <- R6::R6Class("StanModel", #' @param theta_unc The vector of unconstrained parameters. #' @param include_tp Whether to also output the transformed parameters of the model. #' @param include_gq Whether to also output the generated quantities of the model. - #' @param chain_id A chain ID used to offset a PRNG seeded with the model's seed - #' which should be unique between calls. One of `rng` or `seed` must be specified - #' if `include_gq` is `True`. - #' @param rng The random number generator to use. See `StanModel$new_rng()`. - #' One of `rng` or `seed` must be specified if `include_gq` is `True`. + #' @param rng The random number generator to use if `include_gq` is `TRUE`. See `StanModel$new_rng()`. #' @return The constrained parameters of the model. - param_constrain = function(theta_unc, include_tp = FALSE, include_gq = FALSE, chain_id, rng) { - if (missing(chain_id) && missing(rng)){ - if (include_gq) { - stop("Either chain_id or rng must be specified if include_gq is True.") - } else { - chain_id <- 0 + param_constrain = function(theta_unc, include_tp = FALSE, include_gq = FALSE, rng) { + if (missing(rng)) { + if (include_gq){ + stop("A rng must be provided if include_gq is True.") } - } - if (!missing(rng)){ - vars <- .C("bs_param_constrain_R", as.raw(private$model), - as.logical(include_tp), as.logical(include_gq), as.double(theta_unc), - theta = double(self$param_num(include_tp = include_tp, include_gq = include_gq)), - rng = as.raw(rng$rng), - return_code = as.integer(0), - err_msg = as.character(""), - err_ptr = raw(8), - PACKAGE = private$lib_name - ) + rng_ptr <- as.integer(0) } else { - vars <- .C("bs_param_constrain_seeded_R", as.raw(private$model), - as.logical(include_tp), as.logical(include_gq), as.double(theta_unc), - theta = double(self$param_num(include_tp = include_tp, include_gq = include_gq)), - seed = as.integer(private$seed), - chain_id = as.integer(chain_id), - return_code = as.integer(0), - err_msg = as.character(""), - err_ptr = raw(8), - PACKAGE = private$lib_name - ) + rng_ptr <- as.raw(rng$rng) } + vars <- .C("bs_param_constrain_R", as.raw(private$model), + as.logical(include_tp), as.logical(include_gq), as.double(theta_unc), + theta = double(self$param_num(include_tp = include_tp, include_gq = include_gq)), + rng = rng_ptr, + return_code = as.integer(0), + err_msg = as.character(""), + err_ptr = raw(8), + PACKAGE = private$lib_name + ) if (vars$return_code) { stop(handle_error(private$lib_name, vars$err_msg, vars$err_ptr, "param_constrain")) @@ -184,16 +168,10 @@ StanModel <- R6::R6Class("StanModel", }, #' @description #' Create a new persistent PRNG object for use in `param_constrain()`. - #' @param chain_id Identifier for a sequence in the RNG. This should be made - #' a distinct number for each PRNG created with the same seed - #' (for example, 1:N for N PRNGS). - #' @param seed The seed for the PRNG. If this is not specified, the model's seed is used. + #' @param seed The seed for the PRNG. #' @return A `StanRNG` object. - new_rng = function(chain_id, seed) { - if (missing(seed)){ - seed <- private$seed - } - StanRNG$new(private$lib_name, seed, chain_id) + new_rng = function(seed) { + StanRNG$new(private$lib_name, seed) }, #' @description #' Returns a vector of unconstrained parameters give the constrained parameters. @@ -340,14 +318,12 @@ StanRNG <- R6::R6Class("StanRNG", #' Create a StanRng #' @param lib_name The name of the Stan dynamic library. #' @param seed The seed for the RNG. - #' @param chain_id The chain ID for the RNG. #' @return A new StanRNG. - initialize = function(lib_name, seed, chain_id) { + initialize = function(lib_name, seed) { private$lib_name <- lib_name vars <- .C("bs_construct_rng_R", as.integer(seed), - as.integer(chain_id), ptr_out = raw(8), err_msg = as.character(""), err_ptr = raw(8), @@ -363,7 +339,6 @@ StanRNG <- R6::R6Class("StanRNG", rng = NA ), private = list( - lib = NA, lib_name = NA, finalize = function() { .C("bs_destruct_rng_R", diff --git a/R/man/StanModel.Rd b/R/man/StanModel.Rd index 89894be9..2b8822d0 100644 --- a/R/man/StanModel.Rd +++ b/R/man/StanModel.Rd @@ -20,6 +20,7 @@ as well as constraining and unconstraining transforms. \item \href{#method-StanModel-new}{\code{StanModel$new()}} \item \href{#method-StanModel-name}{\code{StanModel$name()}} \item \href{#method-StanModel-model_info}{\code{StanModel$model_info()}} +\item \href{#method-StanModel-model_version}{\code{StanModel$model_version()}} \item \href{#method-StanModel-param_names}{\code{StanModel$param_names()}} \item \href{#method-StanModel-param_unc_names}{\code{StanModel$param_unc_names()}} \item \href{#method-StanModel-param_num}{\code{StanModel$param_num()}} @@ -49,7 +50,7 @@ Create a Stan Model instance. \item{\code{data}}{Either a JSON string literal or a path to a data file in JSON format ending in ".json".} -\item{\code{seed}}{Seed for the RNG in the model object.} +\item{\code{seed}}{Seed for the RNG used in constructing the model.} } \if{html}{\out{}} } @@ -82,6 +83,15 @@ Get compile information about this Stan model. \subsection{Returns}{ A character vector of the Stan version and important flags. } +} +\if{html}{\out{
}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-StanModel-model_version}{}}} +\subsection{Method \code{model_version()}}{ +\subsection{Usage}{ +\if{html}{\out{
}}\preformatted{StanModel$model_version()}\if{html}{\out{
}} +} + } \if{html}{\out{
}} \if{html}{\out{}} @@ -177,7 +187,6 @@ See also \code{StanModel$param_unconstrain()}, the inverse of this function. theta_unc, include_tp = FALSE, include_gq = FALSE, - chain_id, rng )}\if{html}{\out{}} } @@ -191,12 +200,7 @@ See also \code{StanModel$param_unconstrain()}, the inverse of this function. \item{\code{include_gq}}{Whether to also output the generated quantities of the model.} -\item{\code{chain_id}}{A chain ID used to offset a PRNG seeded with the model's seed -which should be unique between calls. One of \code{rng} or \code{seed} must be specified -if \code{include_gq} is \code{True}.} - -\item{\code{rng}}{The random number generator to use. See \code{StanModel$new_rng()}. -One of \code{rng} or \code{seed} must be specified if \code{include_gq} is \code{True}.} +\item{\code{rng}}{The random number generator to use if \code{include_gq} is \code{TRUE}. See \code{StanModel$new_rng()}.} } \if{html}{\out{}} } @@ -210,17 +214,13 @@ The constrained parameters of the model. \subsection{Method \code{new_rng()}}{ Create a new persistent PRNG object for use in \code{param_constrain()}. \subsection{Usage}{ -\if{html}{\out{
}}\preformatted{StanModel$new_rng(chain_id, seed)}\if{html}{\out{
}} +\if{html}{\out{
}}\preformatted{StanModel$new_rng(seed)}\if{html}{\out{
}} } \subsection{Arguments}{ \if{html}{\out{
}} \describe{ -\item{\code{chain_id}}{Identifier for a sequence in the RNG. This should be made -a distinct number for each PRNG created with the same seed -(for example, 1:N for N PRNGS).} - -\item{\code{seed}}{The seed for the PRNG. If this is not specified, the model's seed is used.} +\item{\code{seed}}{The seed for the PRNG.} } \if{html}{\out{
}} } diff --git a/R/man/StanRNG.Rd b/R/man/StanRNG.Rd index f0af6940..12f42bc3 100644 --- a/R/man/StanRNG.Rd +++ b/R/man/StanRNG.Rd @@ -31,7 +31,7 @@ RNG object for use with \code{StanModel$param_constrain()} \subsection{Method \code{new()}}{ Create a StanRng \subsection{Usage}{ -\if{html}{\out{
}}\preformatted{StanRNG$new(lib_name, seed, chain_id)}\if{html}{\out{
}} +\if{html}{\out{
}}\preformatted{StanRNG$new(lib_name, seed)}\if{html}{\out{
}} } \subsection{Arguments}{ @@ -40,8 +40,6 @@ Create a StanRng \item{\code{lib_name}}{The name of the Stan dynamic library.} \item{\code{seed}}{The seed for the RNG.} - -\item{\code{chain_id}}{The chain ID for the RNG.} } \if{html}{\out{}} } diff --git a/R/man/handle_error.Rd b/R/man/handle_error.Rd index 9c434bb2..25b611b4 100644 --- a/R/man/handle_error.Rd +++ b/R/man/handle_error.Rd @@ -4,7 +4,7 @@ \alias{handle_error} \title{Get and free the error message stored at the C++ pointer} \usage{ -handle_error(function_name, err_msg, err_ptr, lib_name) +handle_error(lib_name, err_msg, err_ptr, function_name) } \description{ Get and free the error message stored at the C++ pointer diff --git a/R/tests/testthat/test-bridgestan.R b/R/tests/testthat/test-bridgestan.R index 97cfa683..955e9eab 100644 --- a/R/tests/testthat/test-bridgestan.R +++ b/R/tests/testthat/test-bridgestan.R @@ -115,11 +115,11 @@ test_that("param_constrain handles rng arguments", { expect_equal(4, length(full$param_constrain(c(1.2), include_tp=TRUE, include_gq=TRUE, rng=rng))) # check reproducibility - expect_equal(full$param_constrain(c(1.2), include_gq=TRUE, chain_id=3), - full$param_constrain(c(1.2), include_gq=TRUE, chain_id=3)) + expect_equal(full$param_constrain(c(1.2), include_gq=TRUE, rng=full$new_rng(456)), + full$param_constrain(c(1.2), include_gq=TRUE, rng=full$new_rng(456))) # require at least one present - expect_error(full$param_constrain(c(1.2), include_gq=TRUE), "chain_id or rng must be specified") + expect_error(full$param_constrain(c(1.2), include_gq=TRUE), "rng must be provided") }) @@ -142,5 +142,5 @@ test_that("param_constrain propagates errors", { m2 <- load_model("throw_gq",include_data=FALSE) m2$param_constrain(c(1.2)) # no error - expect_error(m2$param_constrain(c(1.2), include_gq=TRUE, chain_id=1), "find this text: gqfails") + expect_error(m2$param_constrain(c(1.2), include_gq=TRUE, rng=m2$new_rng(1234)), "find this text: gqfails") }) diff --git a/docs/languages/julia.md b/docs/languages/julia.md index e31b2b29..e8a6d205 100644 --- a/docs/languages/julia.md +++ b/docs/languages/julia.md @@ -152,7 +152,7 @@ Return the log density of the specified unconstrained parameters. This calculation drops constant terms that do not depend on the parameters if `propto` is `true` and includes change of variables terms for constrained parameters if `jacobian` is `true`. -source
+source
# **`BridgeStan.log_density_gradient`** — *Function*. @@ -170,7 +170,7 @@ This calculation drops constant terms that do not depend on the parameters if `p This allocates new memory for the gradient output each call. See `log_density_gradient!` for a version which allows re-using existing memory. -source
+source
# **`BridgeStan.log_density_hessian`** — *Function*. @@ -188,7 +188,7 @@ This calculation drops constant terms that do not depend on the parameters if `p This allocates new memory for the gradient and Hessian output each call. See `log_density_gradient!` for a version which allows re-using existing memory. -source
+source
# **`BridgeStan.param_constrain`** — *Function*. @@ -196,19 +196,19 @@ This allocates new memory for the gradient and Hessian output each call. See `lo ```julia -param_constrain(sm, theta_unc, out; include_tp=false, include_gq=false, chain_id=nothing, rng=nothing) +param_constrain(sm, theta_unc, out; include_tp=false, include_gq=false, rng=nothing) ``` Returns a vector constrained parameters given unconstrained parameters. Additionally (if `include_tp` and `include_gq` are set, respectively) returns transformed parameters and generated quantities. -If `include_gq` is set, then either `chain_id` or `rng` must be provided. `chain_id` specifies an offset in a PRNG seeded with the model's seed which should be unique between calls. See `StanRNG` for details on how to construct persistent RNGs. +If `include_gq` is `true`, then `rng` must be provided. See `StanRNG` for details on how to construct RNGs. This allocates new memory for the output each call. See `param_constrain!` for a version which allows re-using existing memory. This is the inverse of `param_unconstrain`. -source
+source
# **`BridgeStan.param_unconstrain`** — *Function*. @@ -228,7 +228,7 @@ This allocates new memory for the output each call. See `param_unconstrain!` for This is the inverse of `param_constrain`. -source
+source
# **`BridgeStan.param_unconstrain_json`** — *Function*. @@ -246,7 +246,7 @@ The JSON is expected to be in the [JSON Format for CmdStan](https://mc-stan.org/ This allocates new memory for the output each call. See `param_unconstrain_json!` for a version which allows re-using existing memory. -source
+source
# **`BridgeStan.name`** — *Function*. @@ -260,7 +260,7 @@ name(sm) Return the name of the model `sm` -source
+source
# **`BridgeStan.model_info`** — *Function*. @@ -276,7 +276,7 @@ Return information about the model `sm`. This includes the Stan version and important compiler flags. -source
+source
# **`BridgeStan.param_num`** — *Function*. @@ -292,7 +292,7 @@ Return the number of (constrained) parameters in the model. This is the total of all the sizes of items declared in the `parameters` block of the model. If `include_tp` or `include_gq` are true, items declared in the `transformed parameters` and `generate quantities` blocks are included, respectively. -source
+source
# **`BridgeStan.param_unc_num`** — *Function*. @@ -308,7 +308,7 @@ Return the number of unconstrained parameters in the model. This function is mainly different from `param_num` when variables are declared with constraints. For example, `simplex[5]` has a constrained size of 5, but an unconstrained size of 4. -source
+source
# **`BridgeStan.param_names`** — *Function*. @@ -326,7 +326,7 @@ For containers, indexes are separated by periods (.). For example, the scalar `a` has indexed name `"a"`, the vector entry `a[1]` has indexed name `"a.1"` and the matrix entry `a[2, 3]` has indexed names `"a.2.3"`. Parameter order of the output is column major and more generally last-index major for containers. -source
+source
# **`BridgeStan.param_unc_names`** — *Function*. @@ -342,7 +342,7 @@ Return the indexed names of the unconstrained parameters. For example, a scalar unconstrained parameter `b` has indexed name `b` and a vector entry `b[3]` has indexed name `b.3`. -source
+source
# **`BridgeStan.log_density_gradient!`** — *Function*. @@ -360,7 +360,7 @@ This calculation drops constant terms that do not depend on the parameters if `p The gradient is stored in the vector `out`, and a reference is returned. See `log_density_gradient` for a version which allocates fresh memory. -source
+source
# **`BridgeStan.log_density_hessian!`** — *Function*. @@ -378,7 +378,7 @@ This calculation drops constant terms that do not depend on the parameters if `p The gradient is stored in the vector `out_grad` and the Hessian is stored in `out_hess` and references are returned. See `log_density_hessian` for a version which allocates fresh memory. -source
+source
# **`BridgeStan.param_constrain!`** — *Function*. @@ -386,19 +386,19 @@ The gradient is stored in the vector `out_grad` and the Hessian is stored in `ou ```julia -param_constrain!(sm, theta_unc, out; include_tp=false, include_gq=false, chain_id=nothing, rng=nothing) +param_constrain!(sm, theta_unc, out; include_tp=false, include_gq=false, rng=nothing) ``` Returns a vector constrained parameters given unconstrained parameters. Additionally (if `include_tp` and `include_gq` are set, respectively) returns transformed parameters and generated quantities. -If `include_gq` is set, then either `chain_id` or `rng` must be provided. `chain_id` specifies an offset in a PRNG seeded with the model's seed which should be unique between calls. See `StanRNG` for details on how to construct persistent RNGs. +If `include_gq` is `true`, then `rng` must be provided. See `StanRNG` for details on how to construct RNGs. The result is stored in the vector `out`, and a reference is returned. See `param_constrain` for a version which allocates fresh memory. This is the inverse of `param_unconstrain!`. -source
+source
# **`BridgeStan.param_unconstrain!`** — *Function*. @@ -418,7 +418,7 @@ The result is stored in the vector `out`, and a reference is returned. See `para This is the inverse of `param_constrain!`. -source
+source
# **`BridgeStan.param_unconstrain_json!`** — *Function*. @@ -436,7 +436,7 @@ The JSON is expected to be in the [JSON Format for CmdStan](https://mc-stan.org/ The result is stored in the vector `out`, and a reference is returned. See `param_unconstrain_json` for a version which allocates fresh memory. -source
+source
# **`BridgeStan.StanRNG`** — *Type*. @@ -444,17 +444,17 @@ The result is stored in the vector `out`, and a reference is returned. See `para ```julia -StanRNG(sm::StanModel, chain_id; seed=nothing) +StanRNG(sm::StanModel, seed) ``` -Construct a StanRNG instance from a `StanModel` instance and a chain ID. This ID serves as an offset in the PRNG stream and should be unique for each PRNG created with the same seed. A seed can be supplied to use in place of the model's seed, which is used by default. +Construct a StanRNG instance from a `StanModel` instance and a seed. This can be used in the `param_constrain` and `param_constrain!` methods when using the generated quantities block. This object is not thread-safe, one should be created per thread. -source
+source
diff --git a/docs/languages/r.md b/docs/languages/r.md index 20c5199e..01b4b045 100644 --- a/docs/languages/r.md +++ b/docs/languages/r.md @@ -72,7 +72,7 @@ _Arguments_ - `data` Either a JSON string literal or a path to a data file in JSON format ending in ".json". - - `rng_seed` Seed for the RNG in the model object. + - `rng_seed` Seed for the RNG used in constructing the model. _Returns_ @@ -210,14 +210,13 @@ _Returns_ **Method** `param_constrain()`: -Returns a vector of constrained parameters given the unconstrained parameters. -parameters See also `StanModel$param_unconstrain()`, the inverse +Returns a vector of constrained parameters given the unconstrained parameters. See also `StanModel$param_unconstrain()`, the inverse of this function. _Usage_ ```R -StanModel$param_constrain(theta_unc, include_tp = FALSE, include_gq = FALSE, chain_id, rng) +StanModel$param_constrain(theta_unc, include_tp = FALSE, include_gq = FALSE, rng) ``` @@ -231,11 +230,8 @@ _Arguments_ - `include_gq` Whether to also output the generated quantities of the model. - - `chain_id` A chain ID used to offset a PRNG seeded with the model's seed - which should be unique between calls. One of `rng` or `chain_id` must be specified if `include_gq` is `True` - - - `rng` StanRNG to use in the model object. See `StanModel$new_rng()`. - One of `rng` or `chain_id` must be specified if `include_gq` is `True` + - `rng` The random number generator to use if `include_gq` is + `TRUE`. See `StanModel$new_rng()`. _Returns_ @@ -251,16 +247,13 @@ Create a new persistent PRNG object for use in `param_constrain()`. _Usage_ ```R -StanModel$new_rng(chain_id, seed) +StanModel$new_rng(seed) ``` _Arguments_ - - `chain_id` Identifier for a sequence in the RNG. This should be made a distinct number for each PRNG created with the sane seed (for example, 1:N for N PRNGS). - - - `seed` The seed for the PRNG. If this is not specified, the model's seed is used. - + - `seed` The seed for the PRNG. _Returns_ diff --git a/julia/src/model.jl b/julia/src/model.jl index 7e26cd35..e9188a5e 100644 --- a/julia/src/model.jl +++ b/julia/src/model.jl @@ -89,12 +89,9 @@ mutable struct StanModel end """ - StanRNG(sm::StanModel, chain_id; seed=nothing) + StanRNG(sm::StanModel, seed) -Construct a StanRNG instance from a `StanModel` instance and a chain ID. -This ID serves as an offset in the PRNG stream and should be unique for each -PRNG created with the same seed. -A seed can be supplied to use in place of the model's seed, which is used by default. +Construct a StanRNG instance from a `StanModel` instance and a seed. This can be used in the `param_constrain` and `param_constrain!` methods when using the generated quantities block. @@ -105,30 +102,23 @@ mutable struct StanRNG lib::Ptr{Nothing} rng::Ptr{StanRNGStruct} seed::UInt32 - chain_id::UInt32 - function StanRNG(sm::StanModel, chain_id; seed = nothing) - seed = convert(UInt32, if isnothing(seed) - sm.seed - else - seed - end) - chain_id = convert(UInt32, chain_id) + function StanRNG(sm::StanModel, seed) + seed = convert(UInt32, seed) err = Ref{Cstring}() rng = ccall( Libc.Libdl.dlsym(sm.lib, "bs_construct_rng"), Ptr{StanModelStruct}, - (UInt32, UInt32, Ref{Cstring}), + (UInt32, Ref{Cstring}), seed, - chain_id, err, ) if rng == C_NULL error(_handle_error(sm.lib, err, "bs_construct_rng")) end - stanrng = new(sm.lib, rng, seed, chain_id) + stanrng = new(sm.lib, rng, seed) function f(stanrng) ccall( @@ -271,16 +261,14 @@ function param_unc_names(sm::StanModel) end """ - param_constrain!(sm, theta_unc, out; include_tp=false, include_gq=false, chain_id=nothing, rng=nothing) + param_constrain!(sm, theta_unc, out; include_tp=false, include_gq=false, rng=nothing) Returns a vector constrained parameters given unconstrained parameters. Additionally (if `include_tp` and `include_gq` are set, respectively) returns transformed parameters and generated quantities. -If `include_gq` is set, then either `chain_id` or `rng` must be provided. -`chain_id` specifies an offset in a PRNG seeded with the model's -seed which should be unique between calls. -See `StanRNG` for details on how to construct persistent RNGs. +If `include_gq` is `true`, then `rng` must be provided. +See `StanRNG` for details on how to construct RNGs. The result is stored in the vector `out`, and a reference is returned. See `param_constrain` for a version which allocates fresh memory. @@ -293,7 +281,6 @@ function param_constrain!( out::Vector{Float64}; include_tp = false, include_gq = false, - chain_id::Union{Int,Nothing} = nothing, rng::Union{StanRNG,Nothing} = nothing, ) dims = param_num(sm; include_tp = include_tp, include_gq = include_gq) @@ -303,67 +290,37 @@ function param_constrain!( ) end - if chain_id === nothing && rng === nothing + if rng === nothing if include_gq - throw( - ArgumentError( - "Must provide either a chain_id or an RNG when including generated quantities", - ), - ) - else - chain_id = 0 + throw(ArgumentError("Must provide an RNG when including generated quantities")) end + rng_ptr = C_NULL + else + rng_ptr = rng.rng end - err = Ref{Cstring}() - if rng !== nothing - rc = ccall( - Libc.Libdl.dlsym(sm.lib, "bs_param_constrain"), + rc = ccall( + Libc.Libdl.dlsym(sm.lib, "bs_param_constrain"), + Cint, + ( + Ptr{StanModelStruct}, Cint, - ( - Ptr{StanModelStruct}, - Cint, - Cint, - Ref{Cdouble}, - Ref{Cdouble}, - Ptr{StanRNGStruct}, - Ref{Cstring}, - ), - sm.stanmodel, - include_tp, - include_gq, - theta_unc, - out, - rng.rng, - err, - ) - else - chain_id = convert(UInt32, chain_id) - rc = ccall( - Libc.Libdl.dlsym(sm.lib, "bs_param_constrain_seeded"), Cint, - ( - Ptr{StanModelStruct}, - Cint, - Cint, - Ref{Cdouble}, - Ref{Cdouble}, - Cuint, - Cuint, - Ref{Cstring}, - ), - sm.stanmodel, - include_tp, - include_gq, - theta_unc, - out, - sm.seed, - chain_id, - err, - ) - end + Ref{Cdouble}, + Ref{Cdouble}, + Ptr{StanRNGStruct}, + Ref{Cstring}, + ), + sm.stanmodel, + include_tp, + include_gq, + theta_unc, + out, + rng_ptr, + err, + ) if rc != 0 error(handle_error(sm.lib, err, "param_constrain")) end @@ -371,16 +328,14 @@ function param_constrain!( end """ - param_constrain(sm, theta_unc, out; include_tp=false, include_gq=false, chain_id=nothing, rng=nothing) + param_constrain(sm, theta_unc, out; include_tp=false, include_gq=false, rng=nothing) Returns a vector constrained parameters given unconstrained parameters. Additionally (if `include_tp` and `include_gq` are set, respectively) returns transformed parameters and generated quantities. -If `include_gq` is set, then either `chain_id` or `rng` must be provided. -`chain_id` specifies an offset in a PRNG seeded with the model's -seed which should be unique between calls. -See `StanRNG` for details on how to construct persistent RNGs. +If `include_gq` is `true`, then `rng` must be provided. +See `StanRNG` for details on how to construct RNGs. This allocates new memory for the output each call. See `param_constrain!` for a version which allows @@ -393,7 +348,6 @@ function param_constrain( theta_unc::Vector{Float64}; include_tp = false, include_gq = false, - chain_id::Union{Int,Nothing} = nothing, rng::Union{StanRNG,Nothing} = nothing, ) out = zeros(param_num(sm, include_tp = include_tp, include_gq = include_gq)) @@ -403,7 +357,6 @@ function param_constrain( out; include_tp = include_tp, include_gq = include_gq, - chain_id = chain_id, rng = rng, ) end diff --git a/julia/test/model_tests.jl b/julia/test/model_tests.jl index b786449c..8b6da210 100644 --- a/julia/test/model_tests.jl +++ b/julia/test/model_tests.jl @@ -176,7 +176,7 @@ end model2 = load_test_model("full", false) - rng = StanRNG(model2, 2; seed = 1234) + rng = StanRNG(model2, 1234) @test 1 == length(BridgeStan.param_constrain(model2, a)) @test 2 == length(BridgeStan.param_constrain(model2, a; include_tp = true)) @test 3 == length(BridgeStan.param_constrain(model2, a; include_gq = true, rng = rng)) @@ -192,8 +192,18 @@ end # reproducibility @test isapprox( - BridgeStan.param_constrain(model2, a; include_gq = true, chain_id = 3), - BridgeStan.param_constrain(model2, a; include_gq = true, chain_id = 3), + BridgeStan.param_constrain( + model2, + a; + include_gq = true, + rng = StanRNG(model2, 45678), + ), + BridgeStan.param_constrain( + model2, + a; + include_gq = true, + rng = StanRNG(model2, 45678), + ), ) # no seed or rng provided @@ -210,12 +220,13 @@ end ) model4 = load_test_model("throw_gq", false) + rng_model4 = StanRNG(model4, 1234) BridgeStan.param_constrain(model4, y) @test_throw_string "find this text: gqfails" BridgeStan.param_constrain( model4, y; include_gq = true, - chain_id = 1, + rng = rng_model4, ) end diff --git a/python/bridgestan/model.py b/python/bridgestan/model.py index 2f5b93d2..d95bb541 100644 --- a/python/bridgestan/model.py +++ b/python/bridgestan/model.py @@ -126,19 +126,6 @@ def __init__( star_star_char, ] - self._param_constrain_seeded = self.stanlib.bs_param_constrain_seeded - self._param_constrain_seeded.restype = ctypes.c_int - self._param_constrain_seeded.argtypes = [ - ctypes.c_void_p, - ctypes.c_int, - ctypes.c_int, - double_array, - double_array, - ctypes.c_uint, - ctypes.c_uint, - star_star_char, - ] - self._param_unconstrain = self.stanlib.bs_param_unconstrain self._param_unconstrain.restype = ctypes.c_int self._param_unconstrain.argtypes = [ @@ -210,7 +197,7 @@ def from_stan_file( """ Construct a StanModel instance from a ``.stan`` file, compiling if necessary. - This is equivalent to calling :func:`bridgestan.compile_model` and then the + This is equivalent to calling :func:`bridgestan.compile_model`` and then the constructor of this class. :param stan_file: A path to a Stan model file. @@ -222,7 +209,7 @@ def from_stan_file( threading for the compiled model. If the same flags are defined in ``make/local``, the versions passed here will take precedent. :param seed: A pseudo random number generator seed. - :raises FileNotFoundError or PermissionError: If `stan_file` does not exist + :raises FileNotFoundError or PermissionError: If ``stan_file`` does not exist or is not readable. :raises ValueError: If BridgeStan cannot be located. :raises RuntimeError: If compilation fails. @@ -274,8 +261,8 @@ def param_num(self, *, include_tp: bool = False, include_gq: bool = False) -> in Return the number of parameters, including transformed parameters and/or generated quantities as indicated. - :param include_tp: `True` to include the transformed parameters. - :param include_gq: `True` to include the generated quantities. + :param include_tp: ``True`` to include the transformed parameters. + :param include_gq: ``True`` to include the generated quantities. :return: The number of parameters. """ return self._param_num(self.model_rng, int(include_tp), int(include_gq)) @@ -295,14 +282,14 @@ def param_names( Return the indexed names of the parameters, including transformed parameters and/or generated quantities as indicated. For containers, indexes are separated by periods (`.`). - For example, the scalar `a` has - indexed name `a`, the vector entry `a[1]` has indexed name `a.1` - and the matrix entry `a[2, 3]` has indexed name `a.2.3`. + For example, the scalar ``a`` has + indexed name ``a``, the vector entry ``a[1]`` has indexed name ``a.1`` + and the matrix entry ``a[2, 3]`` has indexed name ``a.2.3``. Parameter order of the output is column major and more generally last-index major for containers. - :param include_tp: `True` to include transformed parameters. - :param include_gq: `True` to include generated quantities. + :param include_tp: ``True`` to include transformed parameters. + :param include_gq: ``True`` to include generated quantities. :return: The indexed names of the parameters. """ return ( @@ -314,8 +301,8 @@ def param_names( def param_unc_names(self) -> List[str]: """ Return the indexed names of the unconstrained parameters. - For example, a scalar unconstrained parameter `b` has indexed - name `b` and a vector entry `b[3]` has indexed name `b.3`. + For example, a scalar unconstrained parameter ``b`` has indexed + name ``b`` and a vector entry ``b[3]`` has indexed name ``b.3``. :return: The indexed names of the unconstrained parameters. """ @@ -328,7 +315,6 @@ def param_constrain( include_tp: bool = False, include_gq: bool = False, out: Optional[FloatArray] = None, - chain_id: Optional[int] = None, rng: Optional["StanRNG"] = None, ) -> FloatArray: """ @@ -336,37 +322,32 @@ def param_constrain( unconstrained parameters as an array, optionally including the transformed parameters and/or generated quantitities as specified. Including generated quantities uses the PRNG and may update its state. - Setting `out` avoids allocation of a new array for the return value. + Setting ``out`` avoids allocation of a new array for the return value. :param theta_unc: Unconstrained parameter array. - :param include_tp: `True` to include transformed parameters. - :param include_gq: `True` to include generated quantities. + :param include_tp: ``True`` to include transformed parameters. + :param include_gq: ``True`` to include generated quantities. :param out: A location into which the result is stored. If - provided, it must have shape `(D, )`, where `D` is the number of - constrained parameters. If not provided or `None`, a freshly + provided, it must have shape ``(D, )``, where ``D`` is the number of + constrained parameters. If not provided or ``None``, a freshly allocated array is returned. - :param chain_id: A chain ID used to offset a PRNG seeded with the model's - seed which should be unique between calls. One of ``rng`` or - ``chain_id`` must be specified if ``include_gq`` is ``True``. :param rng: A ``StanRNG`` object to use for generating random - numbers, see :meth:`~StanModel.new_rng`. - One of ``rng`` or ``chain_id`` must be specified if ``include_gq`` - is ``True``. + numbers, see :meth:`~StanModel.new_rng``. Must be specified + if ``include_gq`` is ``True``. :return: The constrained parameter array. :raises ValueError: If ``out`` is specified and is not the same shape as the return. - :raises ValueError: If neither ``rng`` nor ``chain_id`` is specified - and ``include_gq`` is ``True``. + :raises ValueError: If ``rng`` is ``None`` and ``include_gq`` is ``True``. :raises RuntimeError: If the C++ Stan model throws an exception. """ - if chain_id is None and rng is None: + if rng is None: if include_gq: raise ValueError( - "Error: must specify rng or chain_id when including generated quantities" + "Error: must specify rng when including generated quantities" ) - else: - # neither specified, but not doing gq, so use a fixed chain id - chain_id = 0 + rng_ptr = None + else: + rng_ptr = rng.ptr dims = self.param_num(include_tp=include_tp, include_gq=include_gq) if out is None: @@ -375,61 +356,46 @@ def param_constrain( raise ValueError( "Error: out must be same size as number of constrained parameters" ) + err = ctypes.pointer(ctypes.c_char_p()) - if rng is not None: - rc = self._param_constrain( - self.model_rng, - int(include_tp), - int(include_gq), - theta_unc, - out, - rng.ptr, - err, - ) - else: - rc = self._param_constrain_seeded( - self.model_rng, - int(include_tp), - int(include_gq), - theta_unc, - out, - self.seed, - chain_id, - err, - ) + rc = self._param_constrain( + self.model_rng, + int(include_tp), + int(include_gq), + theta_unc, + out, + rng_ptr, + err, + ) if rc: raise self._handle_error(err.contents, "param_constrain") return out - def new_rng(self, chain_id: int, *, seed=None) -> "StanRNG": + def new_rng(self, seed) -> "StanRNG": """ - Return a new PRNG for use in :meth:`~StanModel.param_constrain`. + Return a new PRNG for use in :meth:`~StanModel.param_constrain``. - :param chain_id: Identifier for a sequence in the RNG. This should be made - a distinct number for each PRNG created with the same seed - (for example, 1:N for N PRNGS). - :param seed: A seed for the PRNG. If not specified, the - model's seed is used. - :return: A new PRNG for the model. + :param seed: A seed for the PRNG. + :return: A new PRNG wrapper. """ - return StanRNG(self.stanlib, seed or self.seed, chain_id) + return StanRNG(self.stanlib, seed) def param_unconstrain( self, theta: FloatArray, *, out: Optional[FloatArray] = None ) -> FloatArray: """ Return the unconstrained parameters derived from the specified - constrained parameters. Setting `out` avoids allocation of a + constrained parameters. Setting ``out`` avoids allocation of a new array for the return value. :param theta: Constrained parameter array. :param out: A location into which the result is stored. If - provided, it must have shape `(D, )`, where `D` is the number of - unconstrained parameters. If not provided or `None`, a freshly + provided, it must have shape ``(D, )``, where ``D`` is the number of + unconstrained parameters. If not provided or ``None``, a freshly allocated array is returned. - :raises ValueError: If `out` is specified and is not the same + :raises ValueError: If ``out`` is specified and is not the same shape as the return. :raises RuntimeError: If the C++ Stan model throws an exception. """ @@ -456,11 +422,11 @@ def param_unconstrain_json( :param theta_json: The JSON encoded constrained parameters. :param out: A location into which the result is stored. If - provided, it must have shape `(D, )`, where `D` is the number of - unconstrained parameters. If not provided or `None`, a freshly + provided, it must have shape ``(D, )``, where ``D`` is the number of + unconstrained parameters. If not provided or ``None``, a freshly allocated array is returned. :return: The unconstrained parameter array. - :raises ValueError: If `out` is specified and is not the same + :raises ValueError: If ``out`` is specified and is not the same shape as the return value. :raises RuntimeError: If the C++ Stan model throws an exception. """ @@ -488,12 +454,12 @@ def log_density( """ Return the log density of the specified unconstrained parameters, dropping constant terms that do not depend on the - parameters if `propto` is `True` and including change of - variables terms for constrained parameters if `jacobian` is `True`. + parameters if ``propto`` is ``True`` and including change of + variables terms for constrained parameters if ``jacobian`` is ``True``. :param theta_unc: Unconstrained parameter array. - :param propto: `True` if constant terms should be dropped from the log density. - :param jacobian: `True` if change-of-variables terms for + :param propto: ``True`` if constant terms should be dropped from the log density. + :param jacobian: ``True`` if change-of-variables terms for constrained parameters should be included in the log density. :return: The log density. :raises RuntimeError: If the C++ Stan model throws an exception. @@ -518,20 +484,20 @@ def log_density_gradient( """ Return a tuple of the log density and gradient of the specified unconstrained parameters, dropping constant terms that do not depend - on the parameters if `propto` is `True` and including change of - variables terms for constrained parameters if `jacobian` - is `True`. + on the parameters if ``propto`` is ``True`` and including change of + variables terms for constrained parameters if ``jacobian`` + is ``True``. :param theta_unc: Unconstrained parameter array. - :param propto: `True` if constant terms should be dropped from the log density. - :param jacobian: `True` if change-of-variables terms for + :param propto: ``True`` if constant terms should be dropped from the log density. + :param jacobian: ``True`` if change-of-variables terms for constrained parameters should be included in the log density. :param out: A location into which the gradient is stored. If - provided, it must have shape `(D, )` where `D` is the number + provided, it must have shape `(D, )` where ``D`` is the number of parameters. If not provided, a freshly allocated array is returned. :return: A tuple consisting of log density and gradient. - :raises ValueError: If `out` is specified and is not the same + :raises ValueError: If ``out`` is specified and is not the same shape as the gradient. :raises RuntimeError: If the C++ Stan model throws an exception. """ @@ -561,25 +527,25 @@ def log_density_hessian( """ Return a tuple of the log density, gradient, and Hessian of the specified unconstrained parameters, dropping constant terms that do - not depend on the parameters if `propto` is `True` and including + not depend on the parameters if ``propto`` is ``True`` and including change of variables terms for constrained parameters if - `jacobian` is `True`. + ``jacobian`` is ``True``. :param theta_unc: Unconstrained parameter array. - :param propto: `True` if constant terms should be dropped from the log density. - :param jacobian: `True` if change-of-variables terms for + :param propto: ``True`` if constant terms should be dropped from the log density. + :param jacobian: ``True`` if change-of-variables terms for constrained parameters should be included in the log density. :param out_grad: A location into which the gradient is stored. If - provided, it must have shape `(D, )` where `D` is the number + provided, it must have shape `(D, )` where ``D`` is the number of parameters. If not provided, a freshly allocated array is returned. :param out_hess: A location into which the Hessian is stored. If - provided, it must have shape `(D, D)`, where `D` is the + provided, it must have shape `(D, D)`, where ``D`` is the number of parameters. If not provided, a freshly allocated array is returned. :return: A tuple consisting of the log density, gradient, and Hessian. - :raises ValueError: If `out_grad` is specified and is not the - same shape as the gradient or if `out_hess` is specified and it + :raises ValueError: If ``out_grad`` is specified and is not the + same shape as the gradient or if ``out_hess`` is specified and it is not the same shape as the Hessian. :raises RuntimeError: If the C++ Stan model throws an exception. """ @@ -630,13 +596,18 @@ def _handle_error(self, err: ctypes.c_char_p, method: str) -> Exception: class StanRNG: - def __init__(self, lib: ctypes.CDLL, seed: int, chain_id: int) -> None: + def __init__(self, lib: ctypes.CDLL, seed: int) -> None: + """ + Construct a Stan random number generator. + This should not be called directly. Instead, use + :meth:`StanModel.new_rng`. + """ self.stanlib = lib construct = self.stanlib.bs_construct_rng construct.restype = ctypes.c_void_p - construct.argtypes = [ctypes.c_uint, ctypes.c_uint, star_star_char] - self.ptr = construct(seed, chain_id, None) + construct.argtypes = [ctypes.c_uint, star_star_char] + self.ptr = construct(seed, None) if not self.ptr: raise RuntimeError("Failed to construct RNG.") diff --git a/python/test/test_stanmodel.py b/python/test/test_stanmodel.py index 4d319ac1..971a0fb1 100644 --- a/python/test/test_stanmodel.py +++ b/python/test/test_stanmodel.py @@ -208,7 +208,7 @@ def test_param_constrain(): full_so = str(STAN_FOLDER / "full" / "full_model.so") bridge2 = bs.StanModel(full_so) - rng = bridge2.new_rng(chain_id=1, seed=1234) + rng = bridge2.new_rng(seed=1234) np.testing.assert_equal(1, bridge2.param_constrain(a).size) np.testing.assert_equal(2, bridge2.param_constrain(a, include_tp=True).size) @@ -221,8 +221,8 @@ def test_param_constrain(): # reproducibility test np.testing.assert_equal( - bridge2.param_constrain(a, include_gq=True, chain_id=4), - bridge2.param_constrain(a, include_gq=True, chain_id=4), + bridge2.param_constrain(a, include_gq=True, rng=bridge2.new_rng(seed=4567)), + bridge2.param_constrain(a, include_gq=True, rng=bridge2.new_rng(seed=4567)), ) # test error if neither seed or rng is provided @@ -251,7 +251,7 @@ def test_param_constrain(): bridge3 = bs.StanModel(throw_gq_so) bridge3.param_constrain(y, include_gq=False) with pytest.raises(RuntimeError, match="find this text: gqfails"): - bridge3.param_constrain(y, include_gq=True, chain_id=1) + bridge3.param_constrain(y, include_gq=True, rng=bridge3.new_rng(seed=1)) def test_param_unconstrain(): diff --git a/src/bridgestan.cpp b/src/bridgestan.cpp index 9f30f2ad..cbb31f67 100644 --- a/src/bridgestan.cpp +++ b/src/bridgestan.cpp @@ -59,7 +59,18 @@ int bs_param_constrain(const bs_model* m, bool include_tp, bool include_gq, const double* theta_unc, double* theta, bs_rng* rng, char** error_msg) { try { - m->param_constrain(include_tp, include_gq, theta_unc, theta, rng->rng_); + if (rng == nullptr) { + // If RNG is not provided (e.g., we are not using include_gq), use a dummy + // RNG. + // SAFETY: this can be static because we know the rng is never advanced. + static boost::ecuyer1988 dummy_rng(0); + + if (include_gq) + throw std::invalid_argument("include_gq=true but rng=nullptr"); + + m->param_constrain(include_tp, include_gq, theta_unc, theta, dummy_rng); + } else + m->param_constrain(include_tp, include_gq, theta_unc, theta, rng->rng_); return 0; } catch (const std::exception& e) { if (error_msg) { @@ -78,15 +89,6 @@ int bs_param_constrain(const bs_model* m, bool include_tp, bool include_gq, return 1; } -int bs_param_constrain_seeded(const bs_model* m, bool include_tp, - bool include_gq, const double* theta_unc, - double* theta, unsigned int seed, - unsigned int chain_id, char** error_msg) { - bs_rng rng(seed, chain_id); - return bs_param_constrain(m, include_tp, include_gq, theta_unc, theta, &rng, - error_msg); -} - int bs_param_unconstrain(const bs_model* m, const double* theta, double* theta_unc, char** error_msg) { try { @@ -201,10 +203,9 @@ int bs_log_density_hessian(const bs_model* m, bool propto, bool jacobian, return -1; } -bs_rng* bs_construct_rng(unsigned int seed, unsigned int chain_id, - char** error_msg) { +bs_rng* bs_construct_rng(unsigned int seed, char** error_msg) { try { - return new bs_rng(seed, chain_id); + return new bs_rng(seed); } catch (const std::exception& e) { if (error_msg) { std::stringstream error; diff --git a/src/bridgestan.h b/src/bridgestan.h index bfc20c0d..a6f091ff 100644 --- a/src/bridgestan.h +++ b/src/bridgestan.h @@ -34,7 +34,7 @@ bs_model* bs_construct(const char* data_file, unsigned int seed, /** * Destroy the model. * - * @param[in] m pointer to model and RNG structure + * @param[in] m pointer to model structure */ void bs_destruct(bs_model* m); @@ -49,7 +49,7 @@ void bs_free_error_msg(char* error_msg); * Return the name of the specified model as a C-style string. * * The returned string should not be modified; it is freed when the - * model and RNG wrapper is destroyed. + * model wrapper is destroyed. * * @param[in] m pointer to model and RNG structure * @return name of model @@ -60,9 +60,9 @@ const char* bs_name(const bs_model* m); * Return information about the compiled model as a C-style string. * * The returned string should not be modified; it is freed when the - * model and RNG wrapper is destroyed. + * model wrapper is destroyed. * - * @param[in] m pointer to model and RNG structure + * @param[in] m pointer to model structure * @return Information about the model including Stan version, Stan defines, and * compiler flags. */ @@ -80,9 +80,9 @@ const char* bs_model_info(const bs_model* m); * 3]` as `b.2.3`. The numbering follows Stan and is indexed from 1. * * The returned string should not be modified; it is freed when the - * model and RNG wrapper is destroyed. + * model wrapper is destroyed. * - * @param[in] m pointer to model and RNG structure + * @param[in] m pointer to model structure * @param[in] include_tp `true` to include transformed parameters * @param[in] include_gq `true` to include generated quantities * @return CSV-separated, indexed, parameter names @@ -101,9 +101,9 @@ const char* bs_param_names(const bs_model* m, bool include_tp, bool include_gq); * 3]` as `b.2.3`. The numbering follows Stan and is indexed from 1. * * The returned string should not be modified; it is freed when the - * model and RNG wrapper is destroyed. + * model wrapper is destroyed. * - * @param[in] m pointer to model and RNG structure + * @param[in] m pointer to model structure * @return CSV-separated, indexed, unconstrained parameter names */ const char* bs_param_unc_names(const bs_model* m); @@ -113,7 +113,7 @@ const char* bs_param_unc_names(const bs_model* m); * number of transformed parameters and/or generated quantities. * For example, a 2 x 3 matrix counts as 6 scalar parameters. * - * @param[in] m pointer to model and RNG structure + * @param[in] m pointer to model structure * @param[in] include_tp `true` to include transformed parameters * @param[in] include_gq `true` to include generated quantities * @return number of parameters @@ -126,7 +126,7 @@ int bs_param_num(const bs_model* m, bool include_tp, bool include_gq); * parameters if the unconstrained space has fewer dimensions than * the constrained (e.g., for simplexes or correlation matrices). * - * @param[in] m pointer to model and RNG structure + * @param[in] m pointer to model structure * @return number of unconstrained parameters */ int bs_param_unc_num(const bs_model* m); @@ -139,13 +139,14 @@ int bs_param_unc_num(const bs_model* m); * in the Stan program, with multivariate parameters given in * last-index-major order. * - * @param[in] m pointer to model and RNG structure + * @param[in] m pointer to model structure * @param[in] include_tp `true` to include transformed parameters * @param[in] include_gq `true` to include generated quantities * @param[in] theta_unc sequence of unconstrained parameters * @param[out] theta sequence of constrained parameters * @param[in] rng pointer to pseudorandom number generator, should be created - * by `bs_construct_rng` + * by `bs_construct_rng`. This is only required when `include_gq` is `true`, + * otherwise it can be null. * @param[out] error_msg a pointer to a string that will be allocated if there * is an error. This must later be freed by calling `bs_free_error_msg`. * @return code 0 if successful and code -1 if there is an exception @@ -155,36 +156,6 @@ int bs_param_constrain(const bs_model* m, bool include_tp, bool include_gq, const double* theta_unc, double* theta, bs_rng* rng, char** error_msg); -/** - * Set the sequence of constrained parameters based on the specified - * unconstrained parameters, including transformed parameters and/or - * generated quantities as specified, and return a return code of 0 - * for success and -1 for failure. Parameter order is as declared - * in the Stan program, with multivariate parameters given in - * last-index-major order. - * - * This version accepts a chain_id which is used to create a PRNG - * offset from the model's seed which lives only for the duration - * of this call. - * - * @param[in] mr pointer to model and RNG structure - * @param[in] include_tp `true` to include transformed parameters - * @param[in] include_gq `true` to include generated quantities - * @param[in] theta_unc sequence of unconstrained parameters - * @param[out] theta sequence of constrained parameters - * @param[in] chain_id offset for pseudorandom number generator which will be - * created and destroyed during this call. seeded with model seed. See - * `bs_param_constrain` for an option with a persistent RNG. - * @param[out] error_msg a pointer to a string that will be allocated if there - * is an error. This must later be freed by calling `bs_free_error_msg`. - * @return code 0 if successful and code -1 if there is an exception - * in the underlying Stan code - */ -int bs_param_constrain_seeded(const bs_model* mr, bool include_tp, - bool include_gq, const double* theta_unc, - double* theta, unsigned int seed, - unsigned int chain_id, char** error_msg); - /** * Set the sequence of unconstrained parameters based on the * specified constrained parameters, and return a return code of 0 @@ -192,7 +163,7 @@ int bs_param_constrain_seeded(const bs_model* mr, bool include_tp, * in the Stan program, with multivariate parameters given in * last-index-major order. * - * @param[in] m pointer to model and RNG structure + * @param[in] m pointer to model structure * @param[in] theta sequence of constrained parameters * @param[out] theta_unc sequence of unconstrained parameters * @param[out] error_msg a pointer to a string that will be allocated if there @@ -211,7 +182,7 @@ int bs_param_unconstrain(const bs_model* m, const double* theta, * in last-index-major order. The JSON schema assumed is fully * defined in the *CmdStan Reference Manual*. * - * @param[in] m pointer to model and RNG structure + * @param[in] m pointer to model structure * @param[in] json JSON-encoded constrained parameters * @param[out] theta_unc sequence of unconstrained parameters * @param[out] error_msg a pointer to a string that will be allocated if there @@ -229,7 +200,7 @@ int bs_param_unconstrain_json(const bs_model* m, const char* json, * and return a return code of 0 for success and -1 if there is an * exception executing the Stan program. * - * @param[in] m pointer to model and RNG structure + * @param[in] m pointer to model structure * @param[in] propto `true` to discard constant terms * @param[in] jacobian `true` to include change-of-variables terms * @param[in] theta_unc unconstrained parameters @@ -252,7 +223,7 @@ int bs_log_density(const bs_model* m, bool propto, bool jacobian, * * The gradients are computed using automatic differentiation. * - * @param[in] m pointer to model and RNG structure + * @param[in] m pointer to model structure * @param[in] propto `true` to discard constant terms * @param[in] jacobian `true` to include change-of-variables terms * @param[in] theta_unc unconstrained parameters @@ -279,7 +250,7 @@ int bs_log_density_gradient(const bs_model* m, bool propto, bool jacobian, * The gradients are computed using automatic differentiation. the * Hessians are * - * @param[in] m pointer to model and RNG structure + * @param[in] m pointer to model structure * @param[in] propto `true` to discard constant terms * @param[in] jacobian `true` to include change-of-variables terms * @param[in] theta_unc unconstrained parameters @@ -301,12 +272,10 @@ int bs_log_density_hessian(const bs_model* m, bool propto, bool jacobian, * destructed for each thread. * * @param[in] seed seed for the RNG - * @param[in] chain_id identifier for the current sequence of PRNG draws * @param[out] error_msg a pointer to a string that will be allocated if there * is an error. This must later be freed by calling `bs_free_error_msg`. */ -bs_rng* bs_construct_rng(unsigned int seed, unsigned int chain_id, - char** error_msg); +bs_rng* bs_construct_rng(unsigned int seed, char** error_msg); /** * Destruct an RNG object. diff --git a/src/bridgestanR.cpp b/src/bridgestanR.cpp index 73ab69eb..e79df490 100644 --- a/src/bridgestanR.cpp +++ b/src/bridgestanR.cpp @@ -42,16 +42,6 @@ void bs_param_constrain_R(bs_model** model, int* include_tp, int* include_gq, theta, *rng, err_msg); *err_ptr = static_cast(*err_msg); } -void bs_param_constrain_seeded_R(bs_model** model, int* include_tp, - int* include_gq, const double* theta_unc, - double* theta, int* seed, int* chain_id, - int* return_code, char** err_msg, - void** err_ptr) { - *return_code - = bs_param_constrain_seeded(*model, *include_tp, *include_gq, theta_unc, - theta, *seed, *chain_id, err_msg); - *err_ptr = static_cast(*err_msg); -} void bs_param_unconstrain_R(bs_model** model, const double* theta, double* theta_unc, int* return_code, char** err_msg, void** err_ptr) { @@ -87,9 +77,9 @@ void bs_log_density_hessian_R(bs_model** model, int* propto, int* jacobian, val, grad, hess, err_msg); *err_ptr = static_cast(*err_msg); } -void bs_construct_rng_R(int* seed, int* chain_id, bs_rng** ptr_out, +void bs_construct_rng_R(int* seed, bs_rng** ptr_out, char** err_msg, void** err_ptr) { - *ptr_out = bs_construct_rng(*seed, *chain_id, err_msg); + *ptr_out = bs_construct_rng(*seed, err_msg); *err_ptr = static_cast(*err_msg); } void bs_destruct_rng_R(bs_rng** rng) { bs_destruct_rng(*rng); } diff --git a/src/bridgestanR.h b/src/bridgestanR.h index a226bbb3..73f8fe3d 100644 --- a/src/bridgestanR.h +++ b/src/bridgestanR.h @@ -43,12 +43,6 @@ void bs_param_constrain_R(bs_model** model, int* include_tp, int* include_gq, const double* theta_unc, double* theta, bs_rng** rng, int* return_code, char** err_msg, void** err_ptr); -void bs_param_constrain_seeded_R(bs_model** model, int* include_tp, - int* include_gq, const double* theta_unc, - double* theta, int* seed, int* chain_id, - int* return_code, char** err_msg, - void** err_ptr); - void bs_param_unconstrain_R(bs_model** model, const double* theta, double* theta_unc, int* return_code, char** err_msg, void** err_ptr); @@ -71,7 +65,7 @@ void bs_log_density_hessian_R(bs_model** model, int* propto, int* jacobian, double* grad, double* hess, int* return_code, char** err_msg, void** err_ptr); -void bs_construct_rng_R(int* seed, int* chain_id, bs_rng** ptr_out, +void bs_construct_rng_R(int* seed, bs_rng** ptr_out, char** err_msg, void** err_ptr); void bs_destruct_rng_R(bs_rng** rng); diff --git a/src/model_rng.cpp b/src/model_rng.cpp index 3ef28186..48f20ec1 100644 --- a/src/model_rng.cpp +++ b/src/model_rng.cpp @@ -323,7 +323,3 @@ void bs_model::log_density_hessian(bool propto, bool jacobian, Eigen::VectorXd::Map(grad, N) = grad_vec; Eigen::MatrixXd::Map(hessian, N, N) = hess_mat; } - -bs_rng::bs_rng(unsigned int seed, unsigned int chain_id) { - rng_ = stan::services::util::create_rng(seed, chain_id + 1); -} diff --git a/src/model_rng.hpp b/src/model_rng.hpp index 482cf2f0..b1d5ead7 100644 --- a/src/model_rng.hpp +++ b/src/model_rng.hpp @@ -20,8 +20,6 @@ class bs_model { * * @param[in] data_file file from which to read data or "" if none * @param[in] seed pseudorandom number generator seed - * @param[in] chain_id number of gaps to skip in the pseudorandom - * number generator for concurrent computations */ bs_model(const char* data_file, unsigned int seed); @@ -230,7 +228,11 @@ class bs_model { */ class bs_rng { public: - bs_rng(unsigned int seed, unsigned int chain_id); + bs_rng(unsigned int seed) : rng_(seed) { + // discard first value as workaround for + // https://github.com/stan-dev/stan/issues/3167 + rng_.discard(1); + } boost::ecuyer1988 rng_; }; From c9c6d6da93df7cd1cd368d29d228eb74543ca24d Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Thu, 13 Apr 2023 14:56:07 -0400 Subject: [PATCH 17/20] Doc changes and naming consistency --- R/R/bridgestan.R | 10 +++---- docs/languages/julia.md | 40 ++++++++++++++-------------- docs/languages/r.md | 2 +- julia/src/model.jl | 17 ++++++------ python/bridgestan/model.py | 53 ++++++++++++++++---------------------- src/bridgestan.cpp | 2 +- src/bridgestan.h | 10 ++++--- 7 files changed, 64 insertions(+), 70 deletions(-) diff --git a/R/R/bridgestan.R b/R/R/bridgestan.R index 21ab62c4..d7d61da8 100755 --- a/R/R/bridgestan.R +++ b/R/R/bridgestan.R @@ -10,7 +10,7 @@ StanModel <- R6::R6Class("StanModel", #' @description #' Create a Stan Model instance. #' @param lib A path to a compiled BridgeStan Shared Object file. - #' @param data Either a JSON string literal or a path to a data file in JSON format ending in ".json". + #' @param data Either a JSON string literal,a path to a data file in JSON format ending in ".json", or the empty string. #' @param seed Seed for the RNG used in constructing the model. #' @return A new StanModel. initialize = function(lib, data, seed) { @@ -149,7 +149,7 @@ StanModel <- R6::R6Class("StanModel", } rng_ptr <- as.integer(0) } else { - rng_ptr <- as.raw(rng$rng) + rng_ptr <- as.raw(rng$ptr) } vars <- .C("bs_param_constrain_R", as.raw(private$model), as.logical(include_tp), as.logical(include_gq), as.double(theta_unc), @@ -333,16 +333,16 @@ StanRNG <- R6::R6Class("StanRNG", if (all(vars$ptr_out == 0)) { stop(handle_error("construct_rng", vars$err_msg, vars$err_ptr, private$lib_name)) } else { - self$rng <- vars$ptr_out + self$ptr <- vars$ptr_out } }, - rng = NA + ptr = NA ), private = list( lib_name = NA, finalize = function() { .C("bs_destruct_rng_R", - as.raw(self$rng), + as.raw(self$ptr), PACKAGE = private$lib_name ) } diff --git a/docs/languages/julia.md b/docs/languages/julia.md index e8a6d205..f466a68d 100644 --- a/docs/languages/julia.md +++ b/docs/languages/julia.md @@ -121,7 +121,7 @@ StanModel(lib, datafile="", seed=204) A StanModel instance encapsulates a Stan model instantiated with data. -Construct a Stan model from the supplied library file path and data. Data should either be a string containing a JSON string literal or a path to a data file ending in `.json`. If seed is supplied, it is used to initialize the RNG used by the model's constructor. +Construct a Stan model from the supplied library file path and data. Data should either be a string containing a JSON string literal, a path to a data file ending in `.json`, or the empty string. If seed is supplied, it is used to initialize the RNG used by the model's constructor. ``` StanModel(;stan_file, data="", seed=204) @@ -136,7 +136,7 @@ StanModel(;stan_file, stanc_args=[], make_args=[], data="", seed=204) Construct a `StanModel` instance from a `.stan` file. Compilation occurs if no shared object file exists for the supplied Stan file or if a shared object file exists and the Stan file has changed since last compilation. This is equivalent to calling `compile_model` and then the original constructor of `StanModel`. -source
+source
# **`BridgeStan.log_density`** — *Function*. @@ -152,7 +152,7 @@ Return the log density of the specified unconstrained parameters. This calculation drops constant terms that do not depend on the parameters if `propto` is `true` and includes change of variables terms for constrained parameters if `jacobian` is `true`. -source
+source
# **`BridgeStan.log_density_gradient`** — *Function*. @@ -170,7 +170,7 @@ This calculation drops constant terms that do not depend on the parameters if `p This allocates new memory for the gradient output each call. See `log_density_gradient!` for a version which allows re-using existing memory. -source
+source
# **`BridgeStan.log_density_hessian`** — *Function*. @@ -188,7 +188,7 @@ This calculation drops constant terms that do not depend on the parameters if `p This allocates new memory for the gradient and Hessian output each call. See `log_density_gradient!` for a version which allows re-using existing memory. -source
+source
# **`BridgeStan.param_constrain`** — *Function*. @@ -208,7 +208,7 @@ This allocates new memory for the output each call. See `param_constrain!` for a This is the inverse of `param_unconstrain`. -source
+source
# **`BridgeStan.param_unconstrain`** — *Function*. @@ -228,7 +228,7 @@ This allocates new memory for the output each call. See `param_unconstrain!` for This is the inverse of `param_constrain`. -source
+source
# **`BridgeStan.param_unconstrain_json`** — *Function*. @@ -246,7 +246,7 @@ The JSON is expected to be in the [JSON Format for CmdStan](https://mc-stan.org/ This allocates new memory for the output each call. See `param_unconstrain_json!` for a version which allows re-using existing memory. -source
+source
# **`BridgeStan.name`** — *Function*. @@ -260,7 +260,7 @@ name(sm) Return the name of the model `sm` -source
+source
# **`BridgeStan.model_info`** — *Function*. @@ -276,7 +276,7 @@ Return information about the model `sm`. This includes the Stan version and important compiler flags. -source
+source
# **`BridgeStan.param_num`** — *Function*. @@ -292,7 +292,7 @@ Return the number of (constrained) parameters in the model. This is the total of all the sizes of items declared in the `parameters` block of the model. If `include_tp` or `include_gq` are true, items declared in the `transformed parameters` and `generate quantities` blocks are included, respectively. -source
+source
# **`BridgeStan.param_unc_num`** — *Function*. @@ -308,7 +308,7 @@ Return the number of unconstrained parameters in the model. This function is mainly different from `param_num` when variables are declared with constraints. For example, `simplex[5]` has a constrained size of 5, but an unconstrained size of 4. -source
+source
# **`BridgeStan.param_names`** — *Function*. @@ -326,7 +326,7 @@ For containers, indexes are separated by periods (.). For example, the scalar `a` has indexed name `"a"`, the vector entry `a[1]` has indexed name `"a.1"` and the matrix entry `a[2, 3]` has indexed names `"a.2.3"`. Parameter order of the output is column major and more generally last-index major for containers. -source
+source
# **`BridgeStan.param_unc_names`** — *Function*. @@ -342,7 +342,7 @@ Return the indexed names of the unconstrained parameters. For example, a scalar unconstrained parameter `b` has indexed name `b` and a vector entry `b[3]` has indexed name `b.3`. -source
+source
# **`BridgeStan.log_density_gradient!`** — *Function*. @@ -360,7 +360,7 @@ This calculation drops constant terms that do not depend on the parameters if `p The gradient is stored in the vector `out`, and a reference is returned. See `log_density_gradient` for a version which allocates fresh memory. -source
+source
# **`BridgeStan.log_density_hessian!`** — *Function*. @@ -378,7 +378,7 @@ This calculation drops constant terms that do not depend on the parameters if `p The gradient is stored in the vector `out_grad` and the Hessian is stored in `out_hess` and references are returned. See `log_density_hessian` for a version which allocates fresh memory. -source
+source
# **`BridgeStan.param_constrain!`** — *Function*. @@ -398,7 +398,7 @@ The result is stored in the vector `out`, and a reference is returned. See `para This is the inverse of `param_unconstrain!`. -source
+source
# **`BridgeStan.param_unconstrain!`** — *Function*. @@ -418,7 +418,7 @@ The result is stored in the vector `out`, and a reference is returned. See `para This is the inverse of `param_constrain!`. -source
+source
# **`BridgeStan.param_unconstrain_json!`** — *Function*. @@ -436,7 +436,7 @@ The JSON is expected to be in the [JSON Format for CmdStan](https://mc-stan.org/ The result is stored in the vector `out`, and a reference is returned. See `param_unconstrain_json` for a version which allocates fresh memory. -source
+source
# **`BridgeStan.StanRNG`** — *Type*. @@ -454,7 +454,7 @@ This can be used in the `param_constrain` and `param_constrain!` methods when us This object is not thread-safe, one should be created per thread. -source
+source
diff --git a/docs/languages/r.md b/docs/languages/r.md index 01b4b045..0046d749 100644 --- a/docs/languages/r.md +++ b/docs/languages/r.md @@ -70,7 +70,7 @@ _Arguments_ - `lib` A path to a compiled BridgeStan Shared Object file. - - `data` Either a JSON string literal or a path to a data file in JSON format ending in ".json". + - `data` Either a JSON string literal, a path to a data file in JSON format ending in ".json", or the empty string. - `rng_seed` Seed for the RNG used in constructing the model. diff --git a/julia/src/model.jl b/julia/src/model.jl index e9188a5e..981473c8 100644 --- a/julia/src/model.jl +++ b/julia/src/model.jl @@ -18,7 +18,8 @@ end A StanModel instance encapsulates a Stan model instantiated with data. Construct a Stan model from the supplied library file path and data. Data -should either be a string containing a JSON string literal or a path to a data file ending in `.json`. +should either be a string containing a JSON string literal, a path to a data file ending in `.json`, +or the empty string. If seed is supplied, it is used to initialize the RNG used by the model's constructor. StanModel(;stan_file, data="", seed=204) @@ -100,32 +101,32 @@ This object is not thread-safe, one should be created per thread. """ mutable struct StanRNG lib::Ptr{Nothing} - rng::Ptr{StanRNGStruct} + ptr::Ptr{StanRNGStruct} seed::UInt32 function StanRNG(sm::StanModel, seed) seed = convert(UInt32, seed) err = Ref{Cstring}() - rng = ccall( + ptr = ccall( Libc.Libdl.dlsym(sm.lib, "bs_construct_rng"), Ptr{StanModelStruct}, (UInt32, Ref{Cstring}), seed, err, ) - if rng == C_NULL - error(_handle_error(sm.lib, err, "bs_construct_rng")) + if ptr == C_NULL + error(handle_error(sm.lib, err, "bs_construct_rng")) end - stanrng = new(sm.lib, rng, seed) + stanrng = new(sm.lib, ptr, seed) function f(stanrng) ccall( Libc.Libdl.dlsym(stanrng.lib, "bs_destruct_rng"), Cvoid, (Ptr{StanModelStruct},), - stanrng.rng, + stanrng.ptr, ) end @@ -296,7 +297,7 @@ function param_constrain!( end rng_ptr = C_NULL else - rng_ptr = rng.rng + rng_ptr = rng.ptr end err = Ref{Cstring}() diff --git a/python/bridgestan/model.py b/python/bridgestan/model.py index d95bb541..1df674bc 100644 --- a/python/bridgestan/model.py +++ b/python/bridgestan/model.py @@ -20,17 +20,6 @@ class StanModel: A StanModel instance encapsulates a Stan model instantiated with data and provides methods to access parameter names, transforms, log densities, gradients, and Hessians. - - The constructor method instantiates a Stan model and sets constant - return values. The constructor arguments are - - :param model_lib: A path to a compiled shared object. - :param model_data: Either a JSON string literal or a - path to a data file in JSON format ending in ``.json``. - :param seed: A pseudo random number generator seed. - :raises FileNotFoundError or PermissionError: If ``model_lib`` is not readable or - ``model_data`` is specified and not a path to a readable file. - :raises RuntimeError: If there is an error instantiating the Stan model. """ def __init__( @@ -41,13 +30,15 @@ def __init__( seed: int = 1234, ) -> None: """ - Construct a StanModel object for a Stan model and data given + Construct a StanModel object for a compiled Stan model and data given constructor arguments. :param model_lib: A system path to compiled shared object. - :param model_data: Either a JSON string literal or a - system path to a data file in JSON format ending in ``.json``. - :param seed: A pseudo random number generator seed. + :param model_data: Either a JSON string literal, a + system path to a data file in JSON format ending in ``.json``, + or the empty string. + :param seed: A pseudo random number generator seed, used for RNG functions + in the ``transformed data`` block. :raises FileNotFoundError or PermissionError: If ``model_lib`` is not readable or ``model_data`` is specified and not a path to a readable file. :raises RuntimeError: If there is an error instantiating the @@ -74,9 +65,9 @@ def __init__( self._free_error.argtypes = [ctypes.c_char_p] err = ctypes.pointer(ctypes.c_char_p()) - self.model_rng = self._construct(str.encode(self.data_path), self.seed, err) + self.model = self._construct(str.encode(self.data_path), self.seed, err) - if not self.model_rng: + if not self.model: raise self._handle_error(err.contents, "bs_construct") if self.model_version() != __version_info__: @@ -221,8 +212,8 @@ def __del__(self) -> None: """ Destroy the Stan model and free memory. """ - if hasattr(self, "model_rng") and hasattr(self, "_destruct"): - self._destruct(self.model_rng) + if hasattr(self, "model") and hasattr(self, "_destruct"): + self._destruct(self.model) def __repr__(self) -> str: data = f"{self.data_path!r}, " if self.data_path else "" @@ -234,7 +225,7 @@ def name(self) -> str: :return: The name of Stan model. """ - return self._name(self.model_rng).decode("utf-8") + return self._name(self.model).decode("utf-8") def model_info(self) -> str: """ @@ -244,7 +235,7 @@ def model_info(self) -> str: :return: Information about the compiled Stan model. """ - return self._model_info(self.model_rng).decode("utf-8") + return self._model_info(self.model).decode("utf-8") def model_version(self) -> Tuple[int, int, int]: """ @@ -265,7 +256,7 @@ def param_num(self, *, include_tp: bool = False, include_gq: bool = False) -> in :param include_gq: ``True`` to include the generated quantities. :return: The number of parameters. """ - return self._param_num(self.model_rng, int(include_tp), int(include_gq)) + return self._param_num(self.model, int(include_tp), int(include_gq)) def param_unc_num(self) -> int: """ @@ -273,7 +264,7 @@ def param_unc_num(self) -> int: :return: The number of unconstrained parameters. """ - return self._param_unc_num(self.model_rng) + return self._param_unc_num(self.model) def param_names( self, *, include_tp: bool = False, include_gq: bool = False @@ -293,7 +284,7 @@ def param_names( :return: The indexed names of the parameters. """ return ( - self._param_names(self.model_rng, int(include_tp), int(include_gq)) + self._param_names(self.model, int(include_tp), int(include_gq)) .decode("utf-8") .split(",") ) @@ -306,7 +297,7 @@ def param_unc_names(self) -> List[str]: :return: The indexed names of the unconstrained parameters. """ - return self._param_unc_names(self.model_rng).decode("utf-8").split(",") + return self._param_unc_names(self.model).decode("utf-8").split(",") def param_constrain( self, @@ -360,7 +351,7 @@ def param_constrain( err = ctypes.pointer(ctypes.c_char_p()) rc = self._param_constrain( - self.model_rng, + self.model, int(include_tp), int(include_gq), theta_unc, @@ -407,7 +398,7 @@ def param_unconstrain( f"out size = {out.size} != unconstrained params size = {dims}" ) err = ctypes.pointer(ctypes.c_char_p()) - rc = self._param_unconstrain(self.model_rng, theta, out, err) + rc = self._param_unconstrain(self.model, theta, out, err) if rc: raise self._handle_error(err.contents, "param_unconstrain") return out @@ -439,7 +430,7 @@ def param_unconstrain_json( ) chars = theta_json.encode("UTF-8") err = ctypes.pointer(ctypes.c_char_p()) - rc = self._param_unconstrain_json(self.model_rng, chars, out, err) + rc = self._param_unconstrain_json(self.model, chars, out, err) if rc: raise self._handle_error(err.contents, "param_unconstrain_json") return out @@ -467,7 +458,7 @@ def log_density( lp = ctypes.pointer(ctypes.c_double()) err = ctypes.pointer(ctypes.c_char_p()) rc = self._log_density( - self.model_rng, int(propto), int(jacobian), theta_unc, lp, err + self.model, int(propto), int(jacobian), theta_unc, lp, err ) if rc: raise self._handle_error(err.contents, "log_density") @@ -509,7 +500,7 @@ def log_density_gradient( lp = ctypes.pointer(ctypes.c_double()) err = ctypes.pointer(ctypes.c_char_p()) rc = self._log_density_gradient( - self.model_rng, int(propto), int(jacobian), theta_unc, lp, out, err + self.model, int(propto), int(jacobian), theta_unc, lp, out, err ) if rc: raise self._handle_error(err.contents, "log_density_gradient") @@ -564,7 +555,7 @@ def log_density_hessian( lp = ctypes.pointer(ctypes.c_double()) err = ctypes.pointer(ctypes.c_char_p()) rc = self._log_density_hessian( - self.model_rng, + self.model, int(propto), int(jacobian), theta_unc, diff --git a/src/bridgestan.cpp b/src/bridgestan.cpp index 8cee4977..a361df4d 100644 --- a/src/bridgestan.cpp +++ b/src/bridgestan.cpp @@ -225,4 +225,4 @@ bs_rng* bs_construct_rng(unsigned int seed, char** error_msg) { return nullptr; } -void bs_destruct_rng(bs_rng* mr) { delete (mr); } +void bs_destruct_rng(bs_rng* rng) { delete (rng); } diff --git a/src/bridgestan.h b/src/bridgestan.h index d4362799..a9903461 100644 --- a/src/bridgestan.h +++ b/src/bridgestan.h @@ -15,15 +15,17 @@ extern int bs_minor_version; extern int bs_patch_version; /** - * Construct an instance of a model and pseudorandom number - * generator (PRNG) wrapper. Data must be encoded in JSON as - * indicated in the *CmdStan Reference Manual*. + * Construct an instance of a model wrapper. + * Data must be encoded in JSON as indicated + * in the *CmdStan Reference Manual*. * * @param[in] data C-style string. This is either a * path to JSON-encoded data file (must end with ".json"), * a JSON string literal, or nullptr. An empty string or null * pointer are both interpreted as no data. - * @param[in] seed seed for PRNG + * @param[in] seed seed for PRNG used during model construction. + * This PRNG is used for RNG functions in the `transformed data` + * block of the model, and then discarded. * @param[out] error_msg a pointer to a string that will be allocated if there * is an error. This must later be freed by calling `bs_free_error_msg`. * @return pointer to constructed model or `nullptr` if construction From d22191f7318f332c87d29630d472e9aa07ce80dd Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Thu, 13 Apr 2023 15:03:51 -0400 Subject: [PATCH 18/20] Consistently name construct/destruct functions --- R/R/bridgestan.R | 8 ++++---- c-example/example.c | 4 ++-- julia/src/model.jl | 12 ++++++------ julia/test/model_tests.jl | 34 ++++++++++++++++++++++++++++++++++ python/bridgestan/model.py | 10 +++++----- src/bridgestan.cpp | 9 +++++---- src/bridgestan.h | 11 ++++++----- src/bridgestanR.cpp | 14 +++++++------- src/bridgestanR.h | 10 +++++----- 9 files changed, 74 insertions(+), 38 deletions(-) diff --git a/R/R/bridgestan.R b/R/R/bridgestan.R index d7d61da8..6f0bc8a3 100755 --- a/R/R/bridgestan.R +++ b/R/R/bridgestan.R @@ -32,7 +32,7 @@ StanModel <- R6::R6Class("StanModel", } dyn.load(private$lib, PACKAGE = private$lib_name) - ret <- .C("bs_construct_R", + ret <- .C("bs_model_construct_R", as.character(data), as.integer(seed), ptr_out = raw(8), err_msg = as.character(""), @@ -288,7 +288,7 @@ StanModel <- R6::R6Class("StanModel", model = NA, seed = NA, finalize = function() { - .C("bs_destruct_R", + .C("bs_model_destruct_R", as.raw(private$model), PACKAGE = private$lib_name ) @@ -322,7 +322,7 @@ StanRNG <- R6::R6Class("StanRNG", initialize = function(lib_name, seed) { private$lib_name <- lib_name - vars <- .C("bs_construct_rng_R", + vars <- .C("bs_rng_construct_R", as.integer(seed), ptr_out = raw(8), err_msg = as.character(""), @@ -341,7 +341,7 @@ StanRNG <- R6::R6Class("StanRNG", private = list( lib_name = NA, finalize = function() { - .C("bs_destruct_rng_R", + .C("bs_rng_destruct_R", as.raw(self$ptr), PACKAGE = private$lib_name ) diff --git a/c-example/example.c b/c-example/example.c index 352aeeca..fa6d4856 100644 --- a/c-example/example.c +++ b/c-example/example.c @@ -14,7 +14,7 @@ int main(int argc, char** argv) { // this could potentially error, and we may get information back about why. char* err; - bs_model* model = bs_construct(data, 123, &err); + bs_model* model = bs_model_construct(data, 123, &err); if (!model) { if (err) { printf("Error: %s", err); @@ -26,6 +26,6 @@ int main(int argc, char** argv) { printf("This model's name is %s.\n", bs_name(model)); printf("It has %d parameters.\n", bs_param_num(model, 0, 0)); - bs_destruct(model); + bs_model_destruct(model); return 0; } diff --git a/julia/src/model.jl b/julia/src/model.jl index 981473c8..c966abc2 100644 --- a/julia/src/model.jl +++ b/julia/src/model.jl @@ -63,7 +63,7 @@ mutable struct StanModel err = Ref{Cstring}() stanmodel = ccall( - Libc.Libdl.dlsym(lib, "bs_construct"), + Libc.Libdl.dlsym(lib, "bs_model_construct"), Ptr{StanModelStruct}, (Cstring, UInt32, Ref{Cstring}), data, @@ -71,14 +71,14 @@ mutable struct StanModel err, ) if stanmodel == C_NULL - error(handle_error(lib, err, "bs_construct")) + error(handle_error(lib, err, "bs_model_construct")) end sm = new(lib, stanmodel, data, seed) function f(sm) ccall( - Libc.Libdl.dlsym(sm.lib, "bs_destruct"), + Libc.Libdl.dlsym(sm.lib, "bs_model_destruct"), Cvoid, (Ptr{StanModelStruct},), sm.stanmodel, @@ -109,21 +109,21 @@ mutable struct StanRNG err = Ref{Cstring}() ptr = ccall( - Libc.Libdl.dlsym(sm.lib, "bs_construct_rng"), + Libc.Libdl.dlsym(sm.lib, "bs_rng_construct"), Ptr{StanModelStruct}, (UInt32, Ref{Cstring}), seed, err, ) if ptr == C_NULL - error(handle_error(sm.lib, err, "bs_construct_rng")) + error(handle_error(sm.lib, err, "bs_rng_construct")) end stanrng = new(sm.lib, ptr, seed) function f(stanrng) ccall( - Libc.Libdl.dlsym(stanrng.lib, "bs_destruct_rng"), + Libc.Libdl.dlsym(stanrng.lib, "bs_rng_destruct"), Cvoid, (Ptr{StanModelStruct},), stanrng.ptr, diff --git a/julia/test/model_tests.jl b/julia/test/model_tests.jl index 8b6da210..4fb92670 100644 --- a/julia/test/model_tests.jl +++ b/julia/test/model_tests.jl @@ -542,6 +542,40 @@ end end +@testset "threaded model: multi" begin + # Multivariate Gaussian + # make test_models/multi/multi_model.so + + function gaussian(x) + return -0.5 * x' * x + end + + function grad_gaussian(x) + return -x + end + + model = load_test_model("multi") + nt = Threads.nthreads() + + R = 1000 + ld = Vector{Bool}(undef, R) + g = Vector{Bool}(undef, R) + + @sync for it = 1:nt + Threads.@spawn for r = it:nt:R + x = randn(BridgeStan.param_num(model)) + (lp, grad) = BridgeStan.log_density_gradient(model, x) + + ld[r] = isapprox(lp, gaussian(x)) + g[r] = isapprox(grad, grad_gaussian(x)) + end + end + + @test all(ld) + @test all(g) +end + + @testset "gaussian" begin # Guassian with positive constrained standard deviation # make test_models/gaussian/gaussian_model.so diff --git a/python/bridgestan/model.py b/python/bridgestan/model.py index 1df674bc..9e36da2d 100644 --- a/python/bridgestan/model.py +++ b/python/bridgestan/model.py @@ -52,7 +52,7 @@ def __init__( self.data_path = model_data or "" self.seed = seed - self._construct = self.stanlib.bs_construct + self._construct = self.stanlib.bs_model_construct self._construct.restype = ctypes.c_void_p self._construct.argtypes = [ ctypes.c_char_p, @@ -68,7 +68,7 @@ def __init__( self.model = self._construct(str.encode(self.data_path), self.seed, err) if not self.model: - raise self._handle_error(err.contents, "bs_construct") + raise self._handle_error(err.contents, "bs_model_construct") if self.model_version() != __version_info__: warnings.warn( @@ -171,7 +171,7 @@ def __init__( star_star_char, ] - self._destruct = self.stanlib.bs_destruct + self._destruct = self.stanlib.bs_model_destruct self._destruct.restype = None self._destruct.argtypes = [ctypes.c_void_p] @@ -595,7 +595,7 @@ def __init__(self, lib: ctypes.CDLL, seed: int) -> None: """ self.stanlib = lib - construct = self.stanlib.bs_construct_rng + construct = self.stanlib.bs_rng_construct construct.restype = ctypes.c_void_p construct.argtypes = [ctypes.c_uint, star_star_char] self.ptr = construct(seed, None) @@ -603,7 +603,7 @@ def __init__(self, lib: ctypes.CDLL, seed: int) -> None: if not self.ptr: raise RuntimeError("Failed to construct RNG.") - self._destruct = self.stanlib.bs_destruct_rng + self._destruct = self.stanlib.bs_rng_destruct self._destruct.restype = None self._destruct.argtypes = [ctypes.c_void_p] diff --git a/src/bridgestan.cpp b/src/bridgestan.cpp index a361df4d..39e670d6 100644 --- a/src/bridgestan.cpp +++ b/src/bridgestan.cpp @@ -9,7 +9,8 @@ int bs_patch_version = BRIDGESTAN_PATCH; #include -bs_model* bs_construct(const char* data, unsigned int seed, char** error_msg) { +bs_model* bs_model_construct(const char* data, unsigned int seed, + char** error_msg) { try { return new bs_model(data, seed); } catch (const std::exception& e) { @@ -33,7 +34,7 @@ bs_model* bs_construct(const char* data, unsigned int seed, char** error_msg) { return nullptr; } -void bs_destruct(bs_model* m) { delete (m); } +void bs_model_destruct(bs_model* m) { delete (m); } void bs_free_error_msg(char* error_msg) { free(error_msg); } @@ -204,7 +205,7 @@ int bs_log_density_hessian(const bs_model* m, bool propto, bool jacobian, return -1; } -bs_rng* bs_construct_rng(unsigned int seed, char** error_msg) { +bs_rng* bs_rng_construct(unsigned int seed, char** error_msg) { try { return new bs_rng(seed); } catch (const std::exception& e) { @@ -225,4 +226,4 @@ bs_rng* bs_construct_rng(unsigned int seed, char** error_msg) { return nullptr; } -void bs_destruct_rng(bs_rng* rng) { delete (rng); } +void bs_rng_destruct(bs_rng* rng) { delete (rng); } diff --git a/src/bridgestan.h b/src/bridgestan.h index a9903461..6a48e415 100644 --- a/src/bridgestan.h +++ b/src/bridgestan.h @@ -31,14 +31,15 @@ extern int bs_patch_version; * @return pointer to constructed model or `nullptr` if construction * fails */ -bs_model* bs_construct(const char* data, unsigned int seed, char** error_msg); +bs_model* bs_model_construct(const char* data, unsigned int seed, + char** error_msg); /** * Destroy the model. * * @param[in] m pointer to model structure */ -void bs_destruct(bs_model* m); +void bs_model_destruct(bs_model* m); /** * Free the error messages created by other methods. @@ -147,7 +148,7 @@ int bs_param_unc_num(const bs_model* m); * @param[in] theta_unc sequence of unconstrained parameters * @param[out] theta sequence of constrained parameters * @param[in] rng pointer to pseudorandom number generator, should be created - * by `bs_construct_rng`. This is only required when `include_gq` is `true`, + * by `bs_rng_construct`. This is only required when `include_gq` is `true`, * otherwise it can be null. * @param[out] error_msg a pointer to a string that will be allocated if there * is an error. This must later be freed by calling `bs_free_error_msg`. @@ -277,14 +278,14 @@ int bs_log_density_hessian(const bs_model* m, bool propto, bool jacobian, * @param[out] error_msg a pointer to a string that will be allocated if there * is an error. This must later be freed by calling `bs_free_error_msg`. */ -bs_rng* bs_construct_rng(unsigned int seed, char** error_msg); +bs_rng* bs_rng_construct(unsigned int seed, char** error_msg); /** * Destruct an RNG object. * * @param[in] rng pointer to RNG object */ -void bs_destruct_rng(bs_rng* rng); +void bs_rng_destruct(bs_rng* rng); #ifdef __cplusplus } diff --git a/src/bridgestanR.cpp b/src/bridgestanR.cpp index 51c9d091..9cfacf01 100644 --- a/src/bridgestanR.cpp +++ b/src/bridgestanR.cpp @@ -1,9 +1,9 @@ #include "bridgestanR.h" #include "bridgestan.h" -void bs_construct_R(char** data, int* rng, bs_model** ptr_out, char** err_msg, - void** err_ptr) { - *ptr_out = bs_construct(*data, *rng, err_msg); +void bs_model_construct_R(char** data, int* rng, bs_model** ptr_out, + char** err_msg, void** err_ptr) { + *ptr_out = bs_model_construct(*data, *rng, err_msg); *err_ptr = static_cast(*err_msg); } void bs_free_error_msg_R(void** err_msg) { @@ -14,7 +14,7 @@ void bs_version_R(int* major, int* minor, int* patch) { *minor = bs_minor_version; *patch = bs_patch_version; } -void bs_destruct_R(bs_model** model) { bs_destruct(*model); } +void bs_model_destruct_R(bs_model** model) { bs_model_destruct(*model); } void bs_name_R(bs_model** model, char const** name_out) { *name_out = bs_name(*model); } @@ -77,9 +77,9 @@ void bs_log_density_hessian_R(bs_model** model, int* propto, int* jacobian, val, grad, hess, err_msg); *err_ptr = static_cast(*err_msg); } -void bs_construct_rng_R(int* seed, bs_rng** ptr_out, char** err_msg, +void bs_rng_construct_R(int* seed, bs_rng** ptr_out, char** err_msg, void** err_ptr) { - *ptr_out = bs_construct_rng(*seed, err_msg); + *ptr_out = bs_rng_construct(*seed, err_msg); *err_ptr = static_cast(*err_msg); } -void bs_destruct_rng_R(bs_rng** rng) { bs_destruct_rng(*rng); } +void bs_rng_destruct_R(bs_rng** rng) { bs_rng_destruct(*rng); } diff --git a/src/bridgestanR.h b/src/bridgestanR.h index 79476fa9..c6a2e7f4 100644 --- a/src/bridgestanR.h +++ b/src/bridgestanR.h @@ -12,12 +12,12 @@ typedef int bool; // Shim to convert to R interface requirement of void with pointer args // All calls directly delegated to versions without _R suffix -void bs_construct_R(char** data, int* rng, bs_model** ptr_out, char** err_msg, - void** err_ptr); +void bs_model_construct_R(char** data, int* rng, bs_model** ptr_out, + char** err_msg, void** err_ptr); void bs_version_R(int* major, int* minor, int* patch); -void bs_destruct_R(bs_model** model); +void bs_model_destruct_R(bs_model** model); /** * Free error message allocated in C++. Because R performs copies @@ -65,10 +65,10 @@ void bs_log_density_hessian_R(bs_model** model, int* propto, int* jacobian, double* grad, double* hess, int* return_code, char** err_msg, void** err_ptr); -void bs_construct_rng_R(int* seed, bs_rng** ptr_out, char** err_msg, +void bs_rng_construct_R(int* seed, bs_rng** ptr_out, char** err_msg, void** err_ptr); -void bs_destruct_rng_R(bs_rng** rng); +void bs_rng_destruct_R(bs_rng** rng); #ifdef __cplusplus } From 70f20609a8decbdd0b6c62231f84c6241aab7192 Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Thu, 13 Apr 2023 15:37:44 -0400 Subject: [PATCH 19/20] Test thread-safety of param_constrain in Julia --- julia/test/model_tests.jl | 54 +++++++++++++++++++-------------------- 1 file changed, 26 insertions(+), 28 deletions(-) diff --git a/julia/test/model_tests.jl b/julia/test/model_tests.jl index 4fb92670..9c70dd41 100644 --- a/julia/test/model_tests.jl +++ b/julia/test/model_tests.jl @@ -509,8 +509,6 @@ end @testset "threaded model: multi" begin - # Multivariate Gaussian - # make test_models/multi/multi_model.so function gaussian(x) return -0.5 * x' * x @@ -527,8 +525,8 @@ end ld = Vector{Bool}(undef, R) g = Vector{Bool}(undef, R) - @sync for it = 1:nt - Threads.@spawn for r = it:nt:R + @Threads.threads for it = 1:nt + for r = it:nt:R x = randn(BridgeStan.param_num(model)) (lp, grad) = BridgeStan.log_density_gradient(model, x) @@ -542,37 +540,37 @@ end end -@testset "threaded model: multi" begin - # Multivariate Gaussian - # make test_models/multi/multi_model.so - - function gaussian(x) - return -0.5 * x' * x - end +@testset "threaded model: full" begin - function grad_gaussian(x) - return -x - end - - model = load_test_model("multi") + model = load_test_model("full", false) nt = Threads.nthreads() + seeds = rand(UInt32, nt) - R = 1000 - ld = Vector{Bool}(undef, R) - g = Vector{Bool}(undef, R) - - @sync for it = 1:nt - Threads.@spawn for r = it:nt:R - x = randn(BridgeStan.param_num(model)) - (lp, grad) = BridgeStan.log_density_gradient(model, x) + x = [0.5] # bernoulli parameter - ld[r] = isapprox(lp, gaussian(x)) - g[r] = isapprox(grad, grad_gaussian(x)) + R = 1000 + out_size = BridgeStan.param_num(model;include_tp=false, include_gq=true) + + # to test the thread safety of our RNGs, we do two runs + # the first we do in parallel + gq1 = zeros(Float64, R, out_size) + @Threads.threads for it = 1:nt + rng = StanRNG(model,seeds[it]) # RNG is created per-thread + for r = it:nt:R + BridgeStan.param_constrain!(model, x, gq1[r,:]; include_gq=true, rng=rng) + end + end + # the second we do sequentially. + gq2 = zeros(Float64, R, out_size) + for it = 1:nt + rng = StanRNG(model,seeds[it]) + for r = it:nt:R + BridgeStan.param_constrain!(model, x, gq2[r,:]; include_gq=true, rng=rng) end end - @test all(ld) - @test all(g) + # these should be the same if the param_constrain is thread safe + @test gq1 == gq2 end From 852856343201fc6ba3a3addab8005c08ffe86e08 Mon Sep 17 00:00:00 2001 From: "Edward A. Roualdes" Date: Thu, 13 Apr 2023 13:22:33 -0700 Subject: [PATCH 20/20] rework julia test threaded model: full --- julia/test/model_tests.jl | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/julia/test/model_tests.jl b/julia/test/model_tests.jl index 9c70dd41..c1c4d0ed 100644 --- a/julia/test/model_tests.jl +++ b/julia/test/model_tests.jl @@ -553,19 +553,20 @@ end # to test the thread safety of our RNGs, we do two runs # the first we do in parallel - gq1 = zeros(Float64, R, out_size) + gq1 = zeros(Float64, out_size, R) @Threads.threads for it = 1:nt - rng = StanRNG(model,seeds[it]) # RNG is created per-thread - for r = it:nt:R - BridgeStan.param_constrain!(model, x, gq1[r,:]; include_gq=true, rng=rng) - end + rng = StanRNG(model,seeds[it]) # RNG is created per-thread + for r = it:nt:R + gq1[:, r] = BridgeStan.param_constrain(model, x; include_gq=true, rng=rng) end - # the second we do sequentially. - gq2 = zeros(Float64, R, out_size) + end + + # the second we do sequentially + gq2 = zeros(Float64, out_size, R) for it = 1:nt rng = StanRNG(model,seeds[it]) for r = it:nt:R - BridgeStan.param_constrain!(model, x, gq2[r,:]; include_gq=true, rng=rng) + gq2[:, r] = BridgeStan.param_constrain(model, x; include_gq=true, rng=rng) end end