From 4039e191f909a066151f30855a7b634584b3479e Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Sun, 14 Jan 2024 19:41:17 +0200 Subject: [PATCH 1/5] Implement Eigen-to-R wrapping and drop RcppEigen suggest --- DESCRIPTION | 5 ++-- R/fit.R | 5 ++-- R/utils.R | 4 +-- inst/include/hessian.cpp | 4 +-- inst/include/rcpp_eigen_interop.hpp | 42 ++++++++++++++++++++++++++++ man/cmdstanr-package.Rd | 30 ++++++++++++++++++++ man/fit-method-init_model_methods.Rd | 4 +-- 7 files changed, 80 insertions(+), 14 deletions(-) create mode 100644 inst/include/rcpp_eigen_interop.hpp diff --git a/DESCRIPTION b/DESCRIPTION index 053dda7ac..3d55cb486 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -29,7 +29,7 @@ URL: https://mc-stan.org/cmdstanr/, https://discourse.mc-stan.org BugReports: https://github.com/stan-dev/cmdstanr/issues Encoding: UTF-8 LazyData: true -RoxygenNote: 7.2.3 +RoxygenNote: 7.3.0 Roxygen: list(markdown = TRUE, r6 = FALSE) SystemRequirements: CmdStan (https://mc-stan.org/users/interfaces/cmdstan) Depends: @@ -50,6 +50,5 @@ Suggests: loo (>= 2.0.0), rmarkdown, testthat (>= 2.1.0), - Rcpp, - RcppEigen + Rcpp VignetteBuilder: knitr diff --git a/R/fit.R b/R/fit.R index 9dc1c3240..554a79b32 100644 --- a/R/fit.R +++ b/R/fit.R @@ -312,8 +312,8 @@ CmdStanFit$set("public", name = "init", value = init) #' @description The `$init_model_methods()` method compiles and initializes the #' `log_prob`, `grad_log_prob`, `constrain_variables`, `unconstrain_variables` #' and `unconstrain_draws` functions. These are then available as methods of -#' the fitted model object. This requires the additional `Rcpp` and -#' `RcppEigen` packages, which are not required for fitting models using +#' the fitted model object. This requires the additional `Rcpp` package, +#' which are not required for fitting models using #' CmdStanR. #' #' Note: there may be many compiler warnings emitted during compilation but @@ -339,7 +339,6 @@ init_model_methods <- function(seed = 0, verbose = FALSE, hessian = FALSE) { call. = FALSE) } require_suggested_package("Rcpp") - require_suggested_package("RcppEigen") if (length(private$model_methods_env_$hpp_code_) == 0) { stop("Model methods cannot be used with a pre-compiled Stan executable, ", "the model must be compiled again", call. = FALSE) diff --git a/R/utils.R b/R/utils.R index 26d57fec4..0fe42864c 100644 --- a/R/utils.R +++ b/R/utils.R @@ -928,8 +928,7 @@ compile_functions <- function(env, verbose = FALSE, global = FALSE) { mod_stan_funs <- paste(c( env$hpp_code[1:(funs[1] - 1)], "#include ", - "#include ", - "// [[Rcpp::depends(RcppEigen)]]", + "#include ", stan_funs), collapse = "\n") if (global) { @@ -980,7 +979,6 @@ expose_stan_functions <- function(function_env, global = FALSE, verbose = FALSE) call. = FALSE) } require_suggested_package("Rcpp") - require_suggested_package("RcppEigen") if (function_env$compiled) { if (!global) { message("Functions already compiled, nothing to do!") diff --git a/inst/include/hessian.cpp b/inst/include/hessian.cpp index 003c722e0..3883c130a 100644 --- a/inst/include/hessian.cpp +++ b/inst/include/hessian.cpp @@ -1,8 +1,6 @@ -#include +#include #include -// [[Rcpp::depends(RcppEigen)]] - template struct hessian_wrapper { const M& model; diff --git a/inst/include/rcpp_eigen_interop.hpp b/inst/include/rcpp_eigen_interop.hpp new file mode 100644 index 000000000..7ea165d4d --- /dev/null +++ b/inst/include/rcpp_eigen_interop.hpp @@ -0,0 +1,42 @@ +#ifndef CMDSTANR_RCPP_EIGEN_INTEROP_HPP +#define CMDSTANR_RCPP_EIGEN_INTEROP_HPP + +#include +#include + +namespace Rcpp { + namespace traits { + template + class Exporter> { + private: + Rcpp::NumericVector vec_x; + public: + Exporter(SEXP x) : vec_x(x) { } + Eigen::Matrix get() { + if (R == 1) { + return Eigen::Map(vec_x.begin(), vec_x.size()); + } else if (C == 1) { + return Eigen::Map(vec_x.begin(), vec_x.size()); + } else { + Rcpp::Dimension vec_dims(vec_x.attr("dim")); + return Eigen::Map>(vec_x.begin(), vec_dims[0], vec_dims[1]); + } + } + }; + } // traits + + // Rcpp is hardcoded to assume that Eigen types will be wrapped using the + // eigen_wrap function in the RcppEigen package + namespace RcppEigen { + template + SEXP eigen_wrap(const T& x) { + Rcpp::NumericVector vec_rtn(Rcpp::wrap(stan::math::to_array_1d(x))); + if (x.cols() > 1) { + vec_rtn.attr("dim") = Rcpp::Dimension(x.rows(), x.cols()); + } + return vec_rtn; + } + } // RcppEigen +} // Rcpp + +#endif diff --git a/man/cmdstanr-package.Rd b/man/cmdstanr-package.Rd index 930105ee3..c98b0b3a7 100644 --- a/man/cmdstanr-package.Rd +++ b/man/cmdstanr-package.Rd @@ -198,4 +198,34 @@ The Stan and CmdStan documentation: \item Stan documentation: \href{https://mc-stan.org/users/documentation/}{mc-stan.org/users/documentation} \item CmdStan User’s Guide: \href{https://mc-stan.org/docs/cmdstan-guide/}{mc-stan.org/docs/cmdstan-guide} } + +Useful links: +\itemize{ + \item \url{https://mc-stan.org/cmdstanr/} + \item \url{https://discourse.mc-stan.org} + \item Report bugs at \url{https://github.com/stan-dev/cmdstanr/issues} +} + +} +\author{ +\strong{Maintainer}: Jonah Gabry \email{jsg2201@columbia.edu} + +Authors: +\itemize{ + \item Rok Češnovar \email{rok.cesnovar@fri.uni-lj.si} + \item Andrew Johnson (\href{https://orcid.org/0000-0001-7000-8065}{ORCID}) +} + +Other contributors: +\itemize{ + \item Ben Bales [contributor] + \item Mitzi Morris [contributor] + \item Mikhail Popov [contributor] + \item Mike Lawrence [contributor] + \item William Michael Landau \email{will.landau@gmail.com} (\href{https://orcid.org/0000-0003-1878-3253}{ORCID}) [contributor] + \item Jacob Socolar [contributor] + \item Martin Modrák [contributor] + \item Steve Bronder [contributor] +} + } diff --git a/man/fit-method-init_model_methods.Rd b/man/fit-method-init_model_methods.Rd index 0fd0c0752..1a96cd595 100644 --- a/man/fit-method-init_model_methods.Rd +++ b/man/fit-method-init_model_methods.Rd @@ -19,8 +19,8 @@ init_model_methods(seed = 0, verbose = FALSE, hessian = FALSE) The \verb{$init_model_methods()} method compiles and initializes the \code{log_prob}, \code{grad_log_prob}, \code{constrain_variables}, \code{unconstrain_variables} and \code{unconstrain_draws} functions. These are then available as methods of -the fitted model object. This requires the additional \code{Rcpp} and -\code{RcppEigen} packages, which are not required for fitting models using +the fitted model object. This requires the additional \code{Rcpp} package, +which are not required for fitting models using CmdStanR. Note: there may be many compiler warnings emitted during compilation but From 2c8b1bfcf6cb33c9974721018f966b3be6f52081 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Sun, 14 Jan 2024 20:57:03 +0200 Subject: [PATCH 2/5] Fix test failures --- R/utils.R | 2 +- inst/include/rcpp_eigen_interop.hpp | 40 +++++++++++++++++++---------- 2 files changed, 28 insertions(+), 14 deletions(-) diff --git a/R/utils.R b/R/utils.R index 0fe42864c..54c78e3ef 100644 --- a/R/utils.R +++ b/R/utils.R @@ -928,7 +928,7 @@ compile_functions <- function(env, verbose = FALSE, global = FALSE) { mod_stan_funs <- paste(c( env$hpp_code[1:(funs[1] - 1)], "#include ", - "#include ", + "#include ", stan_funs), collapse = "\n") if (global) { diff --git a/inst/include/rcpp_eigen_interop.hpp b/inst/include/rcpp_eigen_interop.hpp index 7ea165d4d..a3fc793eb 100644 --- a/inst/include/rcpp_eigen_interop.hpp +++ b/inst/include/rcpp_eigen_interop.hpp @@ -8,20 +8,33 @@ namespace Rcpp { namespace traits { template class Exporter> { - private: - Rcpp::NumericVector vec_x; - public: - Exporter(SEXP x) : vec_x(x) { } - Eigen::Matrix get() { - if (R == 1) { - return Eigen::Map(vec_x.begin(), vec_x.size()); - } else if (C == 1) { - return Eigen::Map(vec_x.begin(), vec_x.size()); - } else { - Rcpp::Dimension vec_dims(vec_x.attr("dim")); - return Eigen::Map>(vec_x.begin(), vec_dims[0], vec_dims[1]); + private: + SEXP object; + public: + Exporter(SEXP x) : object(x){} + ~Exporter(){} + + Eigen::Matrix get() { + if (R == 1) { + Eigen::Matrix result(::Rf_length(object)); + ::Rcpp::internal::export_indexing, T>(object, result); + return result ; + } else if (C == 1) { + Eigen::Matrix result(::Rf_length(object)); + ::Rcpp::internal::export_indexing, T>(object, result); + return result ; + } else { + Shield dims( ::Rf_getAttrib( object, R_DimSymbol ) ); + if( Rf_isNull(dims) || ::Rf_length(dims) != 2 ){ + throw ::Rcpp::not_a_matrix(); } + int* dims_ = INTEGER(dims); + Eigen::Matrix result(dims_[0],dims_[1]); + T* result_data = result.data(); + ::Rcpp::internal::export_indexing( object, result_data); + return result ; } + } }; } // traits @@ -30,7 +43,8 @@ namespace Rcpp { namespace RcppEigen { template SEXP eigen_wrap(const T& x) { - Rcpp::NumericVector vec_rtn(Rcpp::wrap(stan::math::to_array_1d(x))); + const static int RTYPE = ::Rcpp::traits::r_sexptype_traits>::rtype; + Rcpp::Vector vec_rtn(Rcpp::wrap(stan::math::to_array_1d(x))); if (x.cols() > 1) { vec_rtn.attr("dim") = Rcpp::Dimension(x.rows(), x.cols()); } From 345c1be4c3418fe2a1e8d3ffb9298f3e24f4ef7a Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Sun, 14 Jan 2024 21:03:34 +0200 Subject: [PATCH 3/5] Tweaks --- inst/include/rcpp_eigen_interop.hpp | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/inst/include/rcpp_eigen_interop.hpp b/inst/include/rcpp_eigen_interop.hpp index a3fc793eb..18343e080 100644 --- a/inst/include/rcpp_eigen_interop.hpp +++ b/inst/include/rcpp_eigen_interop.hpp @@ -14,25 +14,26 @@ namespace Rcpp { Exporter(SEXP x) : object(x){} ~Exporter(){} + // Adapted from Rcpp/internal/Exporter.h Eigen::Matrix get() { if (R == 1) { - Eigen::Matrix result(::Rf_length(object)); - ::Rcpp::internal::export_indexing, T>(object, result); - return result ; + Eigen::Matrix result(Rf_length(object)); + Rcpp::internal::export_indexing, T>(object, result); + return result; } else if (C == 1) { - Eigen::Matrix result(::Rf_length(object)); - ::Rcpp::internal::export_indexing, T>(object, result); - return result ; + Eigen::Matrix result(Rf_length(object)); + Rcpp::internal::export_indexing, T>(object, result); + return result; } else { Shield dims( ::Rf_getAttrib( object, R_DimSymbol ) ); - if( Rf_isNull(dims) || ::Rf_length(dims) != 2 ){ + if( Rf_isNull(dims) || Rf_length(dims) != 2 ){ throw ::Rcpp::not_a_matrix(); } int* dims_ = INTEGER(dims); Eigen::Matrix result(dims_[0],dims_[1]); T* result_data = result.data(); - ::Rcpp::internal::export_indexing( object, result_data); - return result ; + Rcpp::internal::export_indexing(object, result_data); + return result; } } }; From 429b648b77d6fb15ea61016c75f9db7bc92dcbc3 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Sun, 14 Jan 2024 21:08:48 +0200 Subject: [PATCH 4/5] Formatting --- inst/include/rcpp_eigen_interop.hpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/inst/include/rcpp_eigen_interop.hpp b/inst/include/rcpp_eigen_interop.hpp index 18343e080..9e6c4b54a 100644 --- a/inst/include/rcpp_eigen_interop.hpp +++ b/inst/include/rcpp_eigen_interop.hpp @@ -25,12 +25,12 @@ namespace Rcpp { Rcpp::internal::export_indexing, T>(object, result); return result; } else { - Shield dims( ::Rf_getAttrib( object, R_DimSymbol ) ); - if( Rf_isNull(dims) || Rf_length(dims) != 2 ){ - throw ::Rcpp::not_a_matrix(); + Shield dims(Rf_getAttrib(object, R_DimSymbol)); + if (Rf_isNull(dims) || Rf_length(dims) != 2) { + throw Rcpp::not_a_matrix(); } int* dims_ = INTEGER(dims); - Eigen::Matrix result(dims_[0],dims_[1]); + Eigen::Matrix result(dims_[0], dims_[1]); T* result_data = result.data(); Rcpp::internal::export_indexing(object, result_data); return result; @@ -44,7 +44,7 @@ namespace Rcpp { namespace RcppEigen { template SEXP eigen_wrap(const T& x) { - const static int RTYPE = ::Rcpp::traits::r_sexptype_traits>::rtype; + const static int RTYPE = Rcpp::traits::r_sexptype_traits>::rtype; Rcpp::Vector vec_rtn(Rcpp::wrap(stan::math::to_array_1d(x))); if (x.cols() > 1) { vec_rtn.attr("dim") = Rcpp::Dimension(x.rows(), x.cols()); From b33dd0294f459b303d9d9776ce2da62e2a4c73f1 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Sun, 14 Jan 2024 22:11:56 +0200 Subject: [PATCH 5/5] Fix handling of 1-row matrices --- inst/include/rcpp_eigen_interop.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/inst/include/rcpp_eigen_interop.hpp b/inst/include/rcpp_eigen_interop.hpp index 9e6c4b54a..f916fff98 100644 --- a/inst/include/rcpp_eigen_interop.hpp +++ b/inst/include/rcpp_eigen_interop.hpp @@ -46,7 +46,7 @@ namespace Rcpp { SEXP eigen_wrap(const T& x) { const static int RTYPE = Rcpp::traits::r_sexptype_traits>::rtype; Rcpp::Vector vec_rtn(Rcpp::wrap(stan::math::to_array_1d(x))); - if (x.cols() > 1) { + if (!stan::is_eigen_col_vector::value) { vec_rtn.attr("dim") = Rcpp::Dimension(x.rows(), x.cols()); } return vec_rtn;