Getting Rcpp to understand stan::math::var

Developers
1

I have been trying to make it possible to pass Stan’s reverse mode autodiff type back and forth between R and C++ by stuffing the value and adjoint into the real and imaginary parts of a complex number, which is a hack but the complex number type is native to R and C++. I was trying to follow @coatless’s example in the attached file (4.1 KB) but failed miserably when I Rcpp::sourceCpp it. That file comments out some stuff that does not compile leaving “one” compiler error, which I do not understand but maybe someone like @Krzysztof_Sakrejda does?

In file included from Rcpp_var.cpp:10:
In file included from /home/ben/r-devel/library/Rcpp/include/RcppCommon.h:166:
/home/ben/r-devel/library/Rcpp/include/Rcpp/internal/export.h:155:4: error: no matching function for call to 'export_indexing__dispatch'
                        export_indexing__dispatch<T,value_type>(
                        ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/home/ben/r-devel/library/Rcpp/include/Rcpp/internal/Exporter.h:64:35: note: in instantiation of function template specialization
      'Rcpp::internal::export_indexing<Eigen::Matrix<stan::math::var, -1, 1, 0, -1, 1>, stan::math::var>' requested here
                ::Rcpp::internal::export_indexing<T,value_type>( object, result ) ;
                                  ^
/home/ben/r-devel/library/Rcpp/include/Rcpp/as.h:89:29: note: in instantiation of member function 'Rcpp::traits::IndexingExporter<Eigen::Matrix<stan::math::var, -1, 1, 0, -1, 1>,
      stan::math::var>::get' requested here
            return exporter.get();
                            ^
/home/ben/r-devel/library/Rcpp/include/Rcpp/as.h:152:26: note: in instantiation of function template specialization 'Rcpp::internal::as<Eigen::Matrix<stan::math::var, -1, 1, 0, -1, 1> >'
      requested here
        return internal::as<T>(x, typename traits::r_type_traits<T>::r_category());
                         ^
Rcpp_var.cpp:138:36: note: in instantiation of function template specialization 'Rcpp::as<Eigen::Matrix<stan::math::var, -1, 1, 0, -1, 1> >' requested here
  stan::math::vector_v v_x = Rcpp::as<stan::math::vector_v>(x);
                                   ^
/home/ben/r-devel/library/Rcpp/include/Rcpp/internal/export.h:131:8: note: candidate function template not viable: no known conversion from
      'typename ::Rcpp::traits::r_type_traits<var>::r_category' (aka 'Rcpp::traits::r_type_generic_tag') to '::Rcpp::traits::r_type_primitive_tag' for 3rd argument
                void export_indexing__dispatch( SEXP x, T& res, ::Rcpp::traits::r_type_primitive_tag ) {
                     ^
/home/ben/r-devel/library/Rcpp/include/Rcpp/internal/export.h:140:8: note: candidate function template not viable: no known conversion from
      'typename ::Rcpp::traits::r_type_traits<var>::r_category' (aka 'Rcpp::traits::r_type_generic_tag') to '::Rcpp::traits::r_type_string_tag' for 3rd argument
                void export_indexing__dispatch( SEXP x, T& res, ::Rcpp::traits::r_type_string_tag ) {
                     ^
2

It would be cool to get this working, I’ll check it out. Could you give one or two examples of what you’d like usage on the R side to look like?

I think that compiler message is telling us that Rcpp::as is peeling off only part of the Eigen type and then doesn’t know what to do with the integers that make up the rest. I’d have to look at implementation to figure out why. Probably will make sense to try it for scalars first and get the behavior right

3

There is a failing test at the end that tries to round trip something from R to C++ and back. But suppose you evaluated normal_lpdf(x, mu, 1) in C++ when mu is a var or vector_v or matrix_v and then you returned mu to R. It might look something like

[1,] -0.1871065-0.5674697i
...
[N,] -0.4945891+2.0543956i

where the real parts hold mu.val() and the imaginary parts hold mu.adj(). We could probably define some convenience functions on the R side like val() and adj() that would wrap the call to Re() and Im().

And then you could go the other way by constructing a complex vector or matrix in R and passing it to a C++ function that would make an Eigen thing of the same shape with var as the scalar type. Except it does not seem straightforward to create a var whose adjoint is anything but 1.0 so I had to do stuff like stan::math::var(real_part / imaginary_part) * imaginary_part;.

4

No idea about Rcpp, but wouldn’t it be easier to use std::pair<double, double> because it’s just a bare container?

5

A pair would come back to R as a list of size two, rather than something numeric. So, possible but it isn’t as easy to work with, for example, a matrix whose elements are each a list of two numbers.

6

I see. Using a vector might be clearer, but it would be suboptimal as it wastes a bunch of memory compared to pair or complex. Now if you could do Eigen::Matrix<double, 2, 1> with hardcoding, I think that’s a tight data structure.

7

At the level we’re working at here we could make a stan::math::var equivalent transparently to a re-classed ‘complex’ number, because in R s3 classes are just tags. It would (with a lot of shenanigans to avoid loosing the class) give us pretty clean representation in R.

> o = as.complex(3.3)
> class(o)
[1] "complex"
> class(o) = 'var'
> str(o)
 'var' cplx 3.3+0i
> class(o)
[1] "var"
val = function(x) Re(x)
adj = function(x) Im(x)
> matrix(data = rep(o,3), nrow=3, ncol=1)
       [,1]
[1,] 3.3+0i
[2,] 3.3+0i
[3,] 3.3+0i

At the last step when you made a matrix you lost the class tag so that would require a wrapper

Anyway, the R side seems ok in this, and the C++ side could stick with the native Stan types to avoid confusion.

I had to sit on this for a while to figure out how it would make sense. Do you think it’s a feature that would get much use? I’m sure I could use it and it would make it plausible to make a real R wrapper for the stan math library but I’m not sure how much use it would get outside of our circles.

8

That was basically my thinking.

9

Great, I think the problems with the C++ as/wrap are just syntax but I’ll have to work through how to implement this anyway. I’ll give it a shot this weekend.