adcomp icon indicating copy to clipboard operation
adcomp copied to clipboard

Namespace declarations to enable Eigen/Tensor in TMB models

Open jeffeaton opened this issue 3 years ago • 5 comments

This PR makes a few small header edits to enable including the Eigen/Tensor module in TMB previously discussed on the TMB mailing list here: https://groups.google.com/g/tmb-users/c/lIsIJ0njMOA

It make the following changes:

  • Explicitly reference tmbutils::array in tmbutils/density.hpp and tmbutils/R_inla.hpp.

    This avoids a namespace clash between Eigen:: and tmbutils:: when Eigen/Tensor headers are included:

    > TMB/include/tmbutils/density.hpp:106:3: error: 
    >    reference to 'array' is ambiguous
    
  • Add #define R_NO_REMAP before #include <Rinternals.h> in TMB.hpp. This avoids a compiler error when including <TMB.hpp> is included before <Tensor> (compiler Apple clang version 13.0.0):

 >  unsupported/Eigen/CXX11/src/Tensor/TensorTraits.h:122:8: error: 
 >    explicit specialization of non-template struct 'Rf_eval'
  • Specify Rf_ prefix in a few locations in tmb_core.hpp and dynamic_data.hpp where it was omitted (Rf_findVar, Rf_install, Rf_setAttrib, Rf_mkString, Rf_ScalarLogical). This was necessitated by the addition of #define R_NO_REMAP. I think these were just oversights because the Rf_ prefix is used in most calls to those functions in those files.

  • Update inst/include/unsupported/ headers to Eigen 3.4.0 for consistency with updated headers in inst/include/Eigen.

I have aimed to be relatively light touch with edits because my understanding of the internals is not very deep. But let me know if there is anything more that I can do to help or clarify.

Thanks for your consideration, and, as ever, many thanks for the brilliant software.

Thanks, Jeff

jeffeaton avatar Feb 13 '22 13:02 jeffeaton

I see in the test failures that adding #define R_NO_REMAP has knock on effects for glmmTMB:

  • https://github.com/glmmTMB/glmmTMB/blob/master/glmmTMB/src/glmmTMB.cpp#L170
  • https://github.com/glmmTMB/glmmTMB/blob/master/glmmTMB/src/glmmTMB.cpp#L171
  • https://github.com/glmmTMB/glmmTMB/blob/master/glmmTMB/src/glmmTMB.cpp#L176
  • https://github.com/glmmTMB/glmmTMB/blob/master/glmmTMB/src/glmmTMB.cpp#L177
  • https://github.com/glmmTMB/glmmTMB/blob/master/glmmTMB/src/distrib.h#L182
  • https://github.com/glmmTMB/glmmTMB/blob/master/glmmTMB/src/distrib.h#L276
  • https://github.com/glmmTMB/glmmTMB/blob/master/glmmTMB/src/distrib.h#L298

and maybe elsewhere. Curious your thoughts about making PRs to add Rf_ prefix there also, or another solution than adding #define R_NO_REMAP here.

jeffeaton avatar Feb 14 '22 08:02 jeffeaton

@jeffeaton I agree that getting Eigen tensors to work in TMB seems like a relevant request. Obviously, I can't merge this PR because it breaks too much stuff, but this is a good starting point - thanks! I'll experiment a bit with this myself for now.

Curious your thoughts about making PRs to add Rf_ prefix there also

It's my impression that error and warning are so widely used that it would be a pain to change. Another solution might be to set R_NO_REMAP and then add a few remaps manually:

#define error Rf_error
#define warning Rf_warning

kaskr avatar Mar 04 '22 16:03 kaskr

Many thanks! Let me know if anything further I can do to help.

jeffeaton avatar Mar 04 '22 17:03 jeffeaton

@jeffeaton

I updated Eigen/unsupported to version 3.4 and can now do (without further changes):

tensor.cpp

#include <TMB.hpp>
#undef eval // Workaround
#include <unsupported/Eigen/CXX11/Tensor>

template<class Type>
Type objective_function<Type>::operator() ()
{
  PARAMETER(p);
  Eigen::Tensor<Type, 2> a(2, 3);
  a.setValues({{p,p+1,p+2},{p+3,p+4,p+5}});
  Rcout << a << "\n\n";
  Eigen::Tensor<Type, 0> b = a.sum();  
  return b(0);
}

tensor.R

library(TMB)
## Compile and load the model
compile("tensor.cpp")
dyn.load(dynlib("tensor"))
## Data and parameters
data <- list()
parameters <- list(p=0)
## Make a function object
obj <- MakeADFun(data, parameters, DLL="tensor")

kaskr avatar Mar 29 '22 09:03 kaskr

Thanks so much for the easy solution on this.

Just leaving a note that to get the example above working on Windows (with Rtools 4.2), I had to add two more #undef statements:

tensor.cpp

#include <TMB.hpp>

#undef eval // Workaround
#undef Realloc
#undef Free

#include <unsupported/Eigen/CXX11/Tensor>

template<class Type>
Type objective_function<Type>::operator() ()
{
  PARAMETER(p);
  Eigen::Tensor<Type, 2> a(2, 3);
  a.setValues({{p,p+1,p+2},{p+3,p+4,p+5}});
  Rcout << a << "\n\n";
  Eigen::Tensor<Type, 0> b = a.sum();  
  return b(0);
}

The resolution was derived from this SO post: https://stackoverflow.com/questions/11588765/using-rcpp-with-windows-specific-includes

jeffeaton avatar Sep 28 '22 16:09 jeffeaton