CODEWITHSUNDEEP
c++eigenrcpp
jds
jds

Reputation: 614

RcppEigen Templated Function to Fill with Unit Normals

I have the following code:

#include <RcppEigen.h>
using namespace Rcpp;
using Eigen::MatrixXd;
using Eigen::VectorXd;
using Eigen::Lower;
using Eigen::Map;

// fills passed dense objects with unit normal random variables
template <typename Derived>
void fillUnitNormal(Eigen::DenseBase<Derived>& Z){
  int m = Z.rows();
  int n = Z.cols();
  NumericVector r(m*n);
  r = rnorm(m*n, 0, 1); // using vectorization from Rcpp sugar
  Map<VectorXd> rvec(as<Map<VectorXd> >(r));
  Map<MatrixXd> rmat(rvec.data(), m, n);
  Z = rmat;
}

That had been working well for me for some time. However, I realized that if Z is a VectorXd object then the function would fail. What would be the proper way to fill each element of an Eigen object inheriting from class Eigen::DenseBase with a normal(0,1) draw?

Upvotes: 2

Views: 75

Answers (1)

Ralf Stubner
Ralf Stubner

Reputation: 26843

One way would be to just std::copy the random values into Z. Since Eigen does not support std::begin(), I decided to use the raw pointer provided by .data(). However, that is not available at the Eigen::DenseBase level. Two levels up in the hierarchy at Eigen::PlainObjectBase it works, though: 

// [[Rcpp::depends(RcppEigen)]]
#include <RcppEigen.h>

// fills passed dense objects with unit normal random variables
template <typename T>
void fillUnitNormal(Eigen::PlainObjectBase<T>& Z){
  int m = Z.rows();
  int n = Z.cols();
  Rcpp::NumericVector r(m*n);
  r = Rcpp::rnorm(m*n, 0, 1); // using vectorization from Rcpp sugar
  std::copy(std::begin(r), std::end(r), Z.data());
}


// [[Rcpp::export]]
Rcpp::List test(int n) {
  Eigen::MatrixXd mat(n, n);
  Eigen::VectorXd vec(n);
  fillUnitNormal(mat);
  fillUnitNormal(vec);
  // gives compile time error: fillUnitNormal(Rcpp::NumericVector::create(n));
  return Rcpp::List::create(mat, vec);
}

/*** R
test(5)
*/

Upvotes: 2

Related Questions