Fix constrain methods and lp adjustments

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) {

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)
   }
 )
 

inst/include/estimator/estimator_ext_header.hpp

@@ -41,14 +41,9 @@ double r_function(const T& v,
                   const Eigen::Map<Eigen::VectorXd>& lower_bounds,
                   const Eigen::Map<Eigen::VectorXd>& 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<double>(internal::ll_fun(v_cons)) + lp;
-  } else {
-    auto v_cons = stan::math::lub_constrain(v, lower_bounds, upper_bounds);
-    return Rcpp::as<double>(internal::ll_fun(v_cons));
-  }
+  double lp = 0;
+  auto v_cons = stan::math::lub_constrain<jacobian__>(v, lower_bounds, upper_bounds, lp);
+  return Rcpp::as<double>(internal::ll_fun(v_cons)) + lp;
 }
 
 template <bool jacobian__, typename T, stan::require_st_var<T>* = nullptr>
@@ -64,19 +59,24 @@ stan::math::var r_function(const T& v,
   auto funwrap = [&](const auto& x) {
     return r_function<jacobian__>(x, finite_diff, no_bounds, bounds_types, lower_bounds, upper_bounds, pstream__);
   };
-  stan::arena_t<Eigen::Matrix<stan::math::var, -1, 1>> arena_v = v;
+  stan::math::var lp(0);
+  stan::arena_t<Eigen::Matrix<stan::math::var, -1, 1>> arena_v;
   stan::arena_t<Eigen::VectorXd> 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<Eigen::VectorXd>(internal::grad_fun(v_cons));
+    arena_v = stan::math::lub_constrain<jacobian__>(v, lower_bounds, upper_bounds, lp);
+    arena_grad = Rcpp::as<Eigen::VectorXd>(internal::grad_fun(arena_v.val()));
+    rtn = Rcpp::as<double>(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

src/call_stan.cpp → src/StanEstimators.cpp

@@ -1,6 +1,7 @@
 #include <Rcpp.h>
 #include <cstdint>
 #include <cmdstan/command.hpp>
+#include <cmdstan/version.hpp>
 #include <estimator/estimator_ext_header.hpp>
 #include <estimator/estimator.hpp>
 
@@ -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<stan::model::model_base> ptr(ext_model_ptr);
   Eigen::Map<Eigen::VectorXd> variables = Rcpp::as<Eigen::Map<Eigen::VectorXd>>(variables_);
-  Eigen::VectorXd unconstrained_variables;
-  ptr->unconstrain_array(variables, unconstrained_variables, &Rcpp::Rcout);
-  return Rcpp::wrap(unconstrained_variables);
+  Eigen::Map<Eigen::VectorXd> lb = Rcpp::as<Eigen::Map<Eigen::VectorXd>>(lb_);
+  Eigen::Map<Eigen::VectorXd> ub = Rcpp::as<Eigen::Map<Eigen::VectorXd>>(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<stan::model::model_base> ptr(ext_model_ptr);
   Eigen::Map<Eigen::MatrixXd> variables = Rcpp::as<Eigen::Map<Eigen::MatrixXd>>(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<Eigen::VectorXd> lb = Rcpp::as<Eigen::Map<Eigen::VectorXd>>(lb_);
+  Eigen::Map<Eigen::VectorXd> ub = Rcpp::as<Eigen::Map<Eigen::VectorXd>>(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<stan::model::model_base> ptr(ext_model_ptr);
-  std::vector<double> upars = Rcpp::as<std::vector<double>>(upars_);
-  std::vector<int> params_i;
-  std::vector<double> 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<Eigen::VectorXd> variables = Rcpp::as<Eigen::Map<Eigen::VectorXd>>(variables_);
+  Eigen::Map<Eigen::VectorXd> lb = Rcpp::as<Eigen::Map<Eigen::VectorXd>>(lb_);
+  Eigen::Map<Eigen::VectorXd> ub = Rcpp::as<Eigen::Map<Eigen::VectorXd>>(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
+}

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}