Merge pull request #899 from stan-dev/drop-rcppeigen-dep

Drop RcppEigen dependency, implement basic Eigen -> C++ interop

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

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)

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 <rcpp_tuple_interop.hpp>",
-    "#include <RcppEigen.h>",
-    "// [[Rcpp::depends(RcppEigen)]]",
+    "#include <rcpp_eigen_interop.hpp>",
     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!")

inst/include/hessian.cpp

@@ -1,8 +1,6 @@
-#include <RcppEigen.h>
+#include <rcpp_eigen_interop.hpp>
 #include <stan/model/hessian.hpp>
 
-// [[Rcpp::depends(RcppEigen)]]
-
 template <class M>
 struct hessian_wrapper {
   const M& model;

inst/include/rcpp_eigen_interop.hpp

@@ -0,0 +1,57 @@
+#ifndef CMDSTANR_RCPP_EIGEN_INTEROP_HPP
+#define CMDSTANR_RCPP_EIGEN_INTEROP_HPP
+
+#include <stan/math/prim/fun/to_array_1d.hpp>
+#include <Rcpp.h>
+
+namespace Rcpp {
+  namespace traits {
+    template <typename T, int R, int C>
+    class Exporter<Eigen::Matrix<T, R, C>> {
+    private:
+        SEXP object;
+    public:
+      Exporter(SEXP x) : object(x){}
+      ~Exporter(){}
+
+      // Adapted from Rcpp/internal/Exporter.h
+      Eigen::Matrix<T, R, C> get() {
+        if (R == 1) {
+          Eigen::Matrix<T, 1, C> result(Rf_length(object));
+          Rcpp::internal::export_indexing<Eigen::Matrix<T, 1, C>, T>(object, result);
+          return result;
+        } else if (C == 1) {
+          Eigen::Matrix<T, R, 1> result(Rf_length(object));
+          Rcpp::internal::export_indexing<Eigen::Matrix<T, R, 1>, T>(object, result);
+          return result;
+        } else {
+          Shield<SEXP> 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<T, R, C> result(dims_[0], dims_[1]);
+          T* result_data = result.data();
+          Rcpp::internal::export_indexing<T*, T>(object, result_data);
+          return result;
+        }
+      }
+    };
+  } // traits
+
+  // Rcpp is hardcoded to assume that Eigen types will be wrapped using the
+  //   eigen_wrap function in the RcppEigen package
+  namespace RcppEigen {
+    template <typename T>
+    SEXP eigen_wrap(const T& x) {
+      const static int RTYPE = Rcpp::traits::r_sexptype_traits<stan::scalar_type_t<T>>::rtype;
+      Rcpp::Vector<RTYPE> vec_rtn(Rcpp::wrap(stan::math::to_array_1d(x)));
+      if (!stan::is_eigen_col_vector<T>::value) {
+        vec_rtn.attr("dim") = Rcpp::Dimension(x.rows(), x.cols());
+      }
+      return vec_rtn;
+    }
+  } // RcppEigen
+} // Rcpp
+
+#endif