CODEWITHSUNDEEP
c++rrcpprcppeigen
Rafa
Rafa

Reputation: 169

Error when using cross() method in RcppEigen

I am trying to implement a simple function that I can call from R to calculate the cross-product of a pair of 3D vectors. To be clear, I mean this cross-product.

I am trying to use for that the cross() method from Eigen library through RcppEigen. My failed attempt currently looks like this:

#include <Rcpp.h>
#include <RcppParallel.h>
#include <RcppEigen.h>
#include <Eigen/Eigen>
#include <Eigen/Geometry>
using namespace Rcpp;
using namespace RcppParallel;
using Eigen::Map;
using Eigen::VectorXd;
using Eigen::Vector3d;
using Eigen::MatrixXd;
// [[Rcpp::plugins("cpp11")]]
// [[Rcpp::export]]
NumericVector vectorCrossProduct3D(NumericVector x, NumericVector y) {
    const Map<VectorXd> xEigen(as<Map<VectorXd> >(x));
    const Map<VectorXd> yEigen(as<Map<VectorXd> >(y));
    VectorXd crossProduct = xEigen.cross(yEigen);
    NumericVector crossProductR = wrap(xEigen);
    return(crossProductR);
}

This is currently resulting in errors of the following type when trying to compile:

error: static_assert failed due to requirement 'Matrix<double, -1, 1, 0, -1, 1>::IsVectorAtCompileTime && Matrix<double, -1, 1, 0, -1, 1>::SizeAtCompileTime == 4' "THIS_METHOD_IS_ONLY_FOR_VECTORS_OF_A_SPECIFIC_SIZE"
      EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(Matrix, 4)

I have also tried declaring the output of the cross() method as MatrixXd, but that also failed. I also tried calling cross3() instead, but also failed.

I feel like I'm missing something very basic here, but I'm new to RcppEigen and can't figure it out... I would greatly appreciate any insights into what is going on

UPDATE

So it seems to work by declaring new variables of type Vector3d and applying cross to those, like this:

// [[Rcpp::plugins("cpp11")]]
// [[Rcpp::export]]
NumericVector vectorCrossProduct3D(NumericVector x, NumericVector y) {
    Map<VectorXd> xEigen(as<Map<VectorXd> >(x));
    Map<VectorXd> yEigen(as<Map<VectorXd> >(y));
    Vector3d xEigen3 = xEigen;
    Vector3d yEigen3 = yEigen;
    Vector3d crossProduct = xEigen3.cross(yEigen3);
    NumericVector crossProductR = wrap(crossProduct);
    return(crossProductR);
}

I guess the compiler prevents from applying cross to vectors of unknown size. However, when I benchmarked it against an implementation of cross product in plain R, it is around twice slower, which does not seem right:

vectorCrossProduct3D_R <- function(u, v) {
    w <- c(u[2]*v[3] - u[3]*v[2],
           u[3]*v[1] - u[1]*v[3],
           u[1]*v[2] - u[2]*v[1])
    return(w)
}

library(microbenchmark)

microbenchmark(vectorCrossProduct3D(1:3, 4:6), vectorCrossProduct3D_R(1:3, 4:6))

Upvotes: 0

Views: 106

Answers (1)

Allan Cameron
Allan Cameron

Reputation: 174468

It's not clear to me why you are using Eigen. Your R code demonstrates that this is a simple arithmetic calculation, and you can simply convert it to C++. Your R version does essentially the same thing as pracma::cross, but without the type and length checks that are necessary in production code.

vectorCrossProduct_R <- function(u, v) {
  w <- c(u[2]*v[3] - u[3]*v[2],
         u[3]*v[1] - u[1]*v[3],
         u[1]*v[2] - u[2]*v[1])
  return(w)
}

Does the same thing as:

Rcpp::cppFunction('
NumericVector vectorCrossProduct_cpp(NumericVector u, NumericVector v) {
  if(v.length() != 3 || u.length() != 3) {
    stop("Input vectors must be length 3");
  }
  NumericVector out = {u[1]*v[2] - u[2]*v[1],
                       u[2]*v[0] - u[0]*v[2],
                       u[0]*v[1] - u[1]*v[0]};
  return out;
}
')

Testing, we have:

vectorCrossProduct_cpp(1:3, 4:6)
#> [1] -3  6 -3
vectorCrossProduct_R(1:3, 4:6)
#> [1] -3  6 -3

However, benchmarking will show little meaningful difference between these two methods. R is relatively slow when iterating through a large vector or list, unless it that iteration is vectorized in the underlying C code. Creation of a single length-three vector by 3 simple calculations is going to be fast.

If you have lots of 3-vectors stored in an existing data structure, then you are likely to gain a significant benefit from iterating through it in C++, but for occasional calls the R version is fine.

Upvotes: 1

Related Questions