diff --git a/R/cpp_exports.R b/R/cpp_exports.R index d8c66ba..680f932 100644 --- a/R/cpp_exports.R +++ b/R/cpp_exports.R @@ -40,17 +40,17 @@ hessian_impl <- function(model_ptr, upars, jacobian = TRUE) { .Call(`hessian_`, model_ptr, upars, jacobian) } -unconstrain_variables_impl <- function(model_ptr, variables) { - .Call(`unconstrain_variables_`, model_ptr, variables) +unconstrain_variables_impl <- function(variables, lb, ub) { + .Call(`unconstrain_variables_`, variables, lb, ub) } -unconstrain_draws_impl <- function(model_ptr, draws, match_format = TRUE) { +unconstrain_draws_impl <- function(draws, lb, ub, match_format = TRUE) { draws_matrix <- posterior::as_draws_matrix(draws) par_cols <- grep("^par", colnames(draws_matrix)) if (length(par_cols) == 0) { stop("No parameter columns found in draws object", call. = FALSE) } - unconstrained_variables <- .Call(`unconstrain_draws_`, model_ptr, draws_matrix[, par_cols]) + unconstrained_variables <- .Call(`unconstrain_draws_`, draws_matrix[, par_cols], lb, ub) draws_matrix[, par_cols] <- unconstrained_variables if (match_format) { match_draws_format(draws, draws_matrix) @@ -59,8 +59,8 @@ unconstrain_draws_impl <- function(model_ptr, draws, match_format = TRUE) { } } -constrain_variables_impl <- function(model_ptr, upars) { - .Call(`constrain_variables_`, model_ptr, upars) +constrain_variables_impl <- function(upars, lb, ub) { + .Call(`constrain_variables_`, upars, lb, ub) } lub_constrain <- function(x, lb, ub) { diff --git a/R/model_methods.R b/R/model_methods.R index 78781f0..ae184eb 100644 --- a/R/model_methods.R +++ b/R/model_methods.R @@ -120,8 +120,7 @@ setMethod("hessian", "StanBase", #' @aliases unconstrain_variables,StanBase,StanBase-method setMethod("unconstrain_variables", "StanBase", function(stan_object, variables) { - unconstrain_variables_impl(stan_object@model_methods$model_pointer, - variables) + unconstrain_variables_impl(variables, stan_object@lower_bounds, stan_object@upper_bounds) } ) @@ -132,7 +131,7 @@ setMethod("unconstrain_draws", "StanBase", if (is.null(draws)) { draws <- stan_object@draws } - unconstrain_draws_impl(stan_object@model_methods$model_pointer, draws) + unconstrain_draws_impl(draws, stan_object@lower_bounds, stan_object@upper_bounds) } ) @@ -142,9 +141,9 @@ setMethod("unconstrain_draws", "StanOptimize", function(stan_object, draws) { if (is.null(draws)) { variables <- stan_object@estimates - unconstrain_draws_impl(stan_object@model_methods$model_pointer, stan_object@estimates, match_format = FALSE) + unconstrain_draws_impl(stan_object@estimates, stan_object@lower_bounds, stan_object@upper_bounds, match_format = FALSE) } else { - unconstrain_draws_impl(stan_object@model_methods$model_pointer, draws) + unconstrain_draws_impl(draws, stan_object@lower_bounds, stan_object@upper_bounds) } } ) @@ -153,8 +152,7 @@ setMethod("unconstrain_draws", "StanOptimize", #' @aliases constrain_variables,StanBase,StanBase-method setMethod("constrain_variables", "StanBase", function(stan_object, unconstrained_variables) { - constrain_variables_impl(stan_object@model_methods$model_pointer, - unconstrained_variables) + constrain_variables_impl(unconstrained_variables, stan_object@lower_bounds, stan_object@upper_bounds) } ) diff --git a/inst/include/estimator/estimator_ext_header.hpp b/inst/include/estimator/estimator_ext_header.hpp index ae29a96..d80867c 100644 --- a/inst/include/estimator/estimator_ext_header.hpp +++ b/inst/include/estimator/estimator_ext_header.hpp @@ -41,14 +41,9 @@ double r_function(const T& v, const Eigen::Map& lower_bounds, const Eigen::Map& upper_bounds, std::ostream* pstream__) { - if (jacobian__) { - double lp = 0; - auto v_cons = stan::math::lub_constrain(v, lower_bounds, upper_bounds, lp); - return Rcpp::as(internal::ll_fun(v_cons)) + lp; - } else { - auto v_cons = stan::math::lub_constrain(v, lower_bounds, upper_bounds); - return Rcpp::as(internal::ll_fun(v_cons)); - } + double lp = 0; + auto v_cons = stan::math::lub_constrain(v, lower_bounds, upper_bounds, lp); + return Rcpp::as(internal::ll_fun(v_cons)) + lp; } template * = nullptr> @@ -64,19 +59,24 @@ stan::math::var r_function(const T& v, auto funwrap = [&](const auto& x) { return r_function(x, finite_diff, no_bounds, bounds_types, lower_bounds, upper_bounds, pstream__); }; - stan::arena_t> arena_v = v; + stan::math::var lp(0); + stan::arena_t> arena_v; stan::arena_t arena_grad; + double rtn; if (finite_diff) { + arena_v = v; arena_grad = fdiff(funwrap, arena_v.val()); + rtn = funwrap(v.val()); } else { - Eigen::VectorXd v_cons = stan::math::lub_constrain(arena_v.val(), lower_bounds, upper_bounds); - arena_grad = Rcpp::as(internal::grad_fun(v_cons)); + arena_v = stan::math::lub_constrain(v, lower_bounds, upper_bounds, lp); + arena_grad = Rcpp::as(internal::grad_fun(arena_v.val())); + rtn = Rcpp::as(internal::ll_fun(arena_v.val())); } return make_callback_var( - funwrap(v.val()), + rtn, [arena_v, arena_grad](auto& vi) mutable { arena_v.adj() += vi.adj() * arena_grad; - }); + }) + lp; } #ifdef USING_R diff --git a/src/call_stan.cpp b/src/StanEstimators.cpp similarity index 81% rename from src/call_stan.cpp rename to src/StanEstimators.cpp index e8a06b2..b15a18a 100644 --- a/src/call_stan.cpp +++ b/src/StanEstimators.cpp @@ -1,6 +1,7 @@ #include #include #include +#include #include #include @@ -241,40 +242,34 @@ RcppExport SEXP hessian_(SEXP ext_model_ptr, SEXP upars_, SEXP jacobian_) { END_RCPP } -RcppExport SEXP unconstrain_variables_(SEXP ext_model_ptr, SEXP variables_) { +RcppExport SEXP unconstrain_variables_(SEXP variables_, SEXP lb_, SEXP ub_) { BEGIN_RCPP - Rcpp::XPtr ptr(ext_model_ptr); Eigen::Map variables = Rcpp::as>(variables_); - Eigen::VectorXd unconstrained_variables; - ptr->unconstrain_array(variables, unconstrained_variables, &Rcpp::Rcout); - return Rcpp::wrap(unconstrained_variables); + Eigen::Map lb = Rcpp::as>(lb_); + Eigen::Map ub = Rcpp::as>(ub_); + return Rcpp::wrap(stan::math::lub_free(variables, lb, ub)); END_RCPP } -RcppExport SEXP unconstrain_draws_(SEXP ext_model_ptr, SEXP draws_matrix_) { +RcppExport SEXP unconstrain_draws_(SEXP draws_matrix_, SEXP lb_, SEXP ub_) { BEGIN_RCPP - Rcpp::XPtr ptr(ext_model_ptr); Eigen::Map variables = Rcpp::as>(draws_matrix_); - Eigen::MatrixXd unconstrained_draws(variables.cols(), variables.rows()); - for (int i = 0; i < variables.rows(); i++) { - Eigen::VectorXd unconstrained_variables; - ptr->unconstrain_array(variables.transpose().col(i), unconstrained_variables, &Rcpp::Rcout); - unconstrained_draws.col(i) = unconstrained_variables; + Eigen::Map lb = Rcpp::as>(lb_); + Eigen::Map ub = Rcpp::as>(ub_); + Eigen::MatrixXd unconstrained_draws(variables.rows(), variables.cols()); + for (int i = 0; i < variables.cols(); i++) { + unconstrained_draws.col(i) = stan::math::lub_free(variables.col(i), lb[i], ub[i]); } - return Rcpp::wrap(unconstrained_draws.transpose()); + return Rcpp::wrap(unconstrained_draws); END_RCPP } -RcppExport SEXP constrain_variables_(SEXP ext_model_ptr, SEXP upars_) { +RcppExport SEXP constrain_variables_(SEXP variables_, SEXP lb_, SEXP ub_) { BEGIN_RCPP - Rcpp::XPtr ptr(ext_model_ptr); - std::vector upars = Rcpp::as>(upars_); - std::vector params_i; - std::vector pars_constrained; - // RNG only used for *_rng calls in generated_quantities, which we don't use - stan::rng_t dummy_rng(0); - ptr->write_array(dummy_rng, upars, params_i, pars_constrained, false, false); - return Rcpp::wrap(pars_constrained); + Eigen::Map variables = Rcpp::as>(variables_); + Eigen::Map lb = Rcpp::as>(lb_); + Eigen::Map ub = Rcpp::as>(ub_); + return Rcpp::wrap(stan::math::lub_constrain(variables, lb, ub)); END_RCPP } @@ -297,3 +292,28 @@ RcppExport SEXP lub_free_(SEXP y_, SEXP lb_, SEXP ub_) { return Rcpp::wrap(stan::math::lub_free(y, lb, ub)); END_RCPP } + +RcppExport SEXP stan_versions_() { + BEGIN_RCPP + std::string math_version = + stan::math::MAJOR_VERSION + "." + + stan::math::MINOR_VERSION + "." + + stan::math::PATCH_VERSION; + + std::string stan_version = + stan::MAJOR_VERSION + "." + + stan::MINOR_VERSION + "." + + stan::PATCH_VERSION; + + std::string cmdstan_version = + cmdstan::MAJOR_VERSION + "." + + cmdstan::MINOR_VERSION + "." + + cmdstan::PATCH_VERSION; + + return Rcpp::List::create( + Rcpp::Named("Math") = math_version, + Rcpp::Named("Stan") = stan_version, + Rcpp::Named("CmdStan") = cmdstan_version + ); + END_RCPP +} diff --git a/src/init.cpp b/src/init.cpp index 4af6399..88f3c25 100644 --- a/src/init.cpp +++ b/src/init.cpp @@ -18,9 +18,9 @@ SEXP make_model_pointer_(SEXP data_json_string_, SEXP seed_); SEXP log_prob_(SEXP ext_model_ptr, SEXP upars_, SEXP jacobian_); SEXP grad_log_prob_(SEXP ext_model_ptr, SEXP upars_, SEXP jacobian_); SEXP hessian_(SEXP ext_model_ptr, SEXP upars_, SEXP jacobian_); -SEXP unconstrain_variables_(SEXP ext_model_ptr, SEXP cons_json_string_); -SEXP unconstrain_draws_(SEXP ext_model_ptr, SEXP draws_matrix_); -SEXP constrain_variables_(SEXP ext_model_ptr, SEXP upars_); +SEXP unconstrain_variables_(SEXP variables_, SEXP lb_, SEXP ub_); +SEXP unconstrain_draws_(SEXP draws_matrix_, SEXP lb_, SEXP ub_); +SEXP constrain_variables_(SEXP variables_, SEXP lb_, SEXP ub_); SEXP lub_constrain_(SEXP y_, SEXP lb_, SEXP ub_); SEXP lub_free_(SEXP y_, SEXP lb_, SEXP ub_); @@ -39,9 +39,9 @@ static const R_CallMethodDef CallEntries[] = { CALLDEF(log_prob_, 3), CALLDEF(grad_log_prob_, 3), CALLDEF(hessian_, 3), - CALLDEF(unconstrain_variables_, 2), - CALLDEF(unconstrain_draws_, 2), - CALLDEF(constrain_variables_, 2), + CALLDEF(unconstrain_variables_, 3), + CALLDEF(unconstrain_draws_, 3), + CALLDEF(constrain_variables_, 3), CALLDEF(lub_constrain_, 3), CALLDEF(lub_free_, 3), {NULL, NULL, 0} diff --git a/src/stan_versions.cpp b/src/stan_versions.cpp deleted file mode 100644 index 9596659..0000000 --- a/src/stan_versions.cpp +++ /dev/null @@ -1,29 +0,0 @@ -#include -#include -#include -#include - -RcppExport SEXP stan_versions_() { - BEGIN_RCPP - std::string math_version = - stan::math::MAJOR_VERSION + "." + - stan::math::MINOR_VERSION + "." + - stan::math::PATCH_VERSION; - - std::string stan_version = - stan::MAJOR_VERSION + "." + - stan::MINOR_VERSION + "." + - stan::PATCH_VERSION; - - std::string cmdstan_version = - cmdstan::MAJOR_VERSION + "." + - cmdstan::MINOR_VERSION + "." + - cmdstan::PATCH_VERSION; - - return Rcpp::List::create( - Rcpp::Named("Math") = math_version, - Rcpp::Named("Stan") = stan_version, - Rcpp::Named("CmdStan") = cmdstan_version - ); - END_RCPP -}