[Rcpp-devel] Expose cpp class with built-in boost odeint to R

Simon Woodward Simon.Woodward at dairynz.co.nz
Fri Sep 14 05:17:12 CEST 2018


I am developing a simulation model in C++ and I want to expose it to R for testing. The model is in a C++ class and uses the boost::odeint library for numerical integration. I am also using an unordered_map to expose the model variables to the user.
I am having some trouble exposing the class functionality to C++, either as separate // [[Rcpp::export]] methods or using RCPP_MODULE. The following code successfully instantiates the model and returns the unordered_map but crashes Rcpp if I try to access the advance(t,dt) method.

# sample code embedded in an R script
library(Rcpp)

Sys.setenv("PKG_CXXFLAGS"='-I"C:/boost/boost_1_66_0"')

sourceCpp(code='
          #include <unordered_map>
          #include <string>
          #define BOOST_CONFIG_SUPPRESS_OUTDATED_MESSAGE
          #include <boost/numeric/odeint.hpp>
          #include <Rcpp.h>

          // Enable C++11 via this plugin (Rcpp 0.10.3 or later)
          // [[Rcpp::plugins(cpp11)]]

          using namespace std;
          using namespace boost::numeric::odeint;
          using namespace Rcpp;

          class spring{

          public:

          // unordered_map gives user efficient access to variables by name
          // https://stackoverflow.com/questions/8483985/obtaining-list-of-keys-and-values-from-unordered-map
          // will be more efficient if we intelligently reserve storage for the map at construction time
          // https://stackoverflow.com/questions/15135900/how-to-pass-a-mapdouble-t-from-c-to-r-with-rcpp
          unordered_map< string , double > variable;

          private:

          // declare state_type and system_time
          typedef boost::array< double , 3 > state_type;
          double system_time;
          const int system_variables = 12; // number of variables available to user

          // declare boost::odeint stepper
          typedef runge_kutta4< state_type > stepper_type;
          stepper_type stepper;

          // state variables
          double time;
          double xd;
          double x;

          // constants
          const double k = 0.12;
          const double a = 1.0;
          const double w = 1.0;
          const double g = 9.81;
          const double mass = 0.03;
          const double xdic = 0.0;
          const double xic = 0.0;

          // other variables
          double xdd = -999;

          public:

          // constructor
          spring(){

          // reserve buckets to minimise storage and avoid rehashing
          // note: only need variables the user needs to set or get
          variable.reserve( system_variables );

          }

          void initialise_model( double a_system_time ){

          // initialise system time
          system_time = a_system_time;

          // initial calculations

          // initial conditions
          time = 0.0;
          xd = xdic;
          x = xic;

          }

          void pull_variables_from_model(){

          // generic system property
          variable["system_time"] = system_time;

          // state variables
          variable["time"] = time;
          variable["xd"] = xd;
          variable["x"] = x;

          // constants
          variable["k"] = k;
          variable["a"] = a;
          variable["w"] = w;
          variable["g"] = g;
          variable["mass"] = mass;
          variable["xdic"] = xdic;
          variable["xic"] = xic;

          // other variables
          variable["xdd"] = xdd;

          }

          void push_variables_to_model(){

          // generic system property
          system_time = variable["system_time"];

          // state variables
          time = variable["time"];
          xd = variable["xd"];
          x = variable["x"];

          // constants (cant change these)
          //                 k = variable["k"];
          //                 a = variable["a"];
          //                 w = variable["w"];
          //                 g = variable["g"];
          //                 mass = variable["mass"];
          //                 xdic = variable["xdic"];
          //                 xic = variable["xic"];

          // other variables
          xdd = variable["xdd"];

          }

          private:

          state_type get_state(){

          state_type a_state;

          // return current state
          a_state[0] = time;
          a_state[1] = xd;
          a_state[2] = x;

          return( a_state );

          }

          void set_state( state_type a_state ){

          // set state
          time = a_state[0];
          xd = a_state[1];
          x = a_state[2];

          }

          public:

          void calculate_rate(){

          // calculations
          xdd = ( mass * g - k * xd - a * x ) / mass;

          }

          public:

          // called by boost::odeint::integrate()
          void operator()( const state_type &a_state , state_type &a_rate, double a_time ){

          // set state
          system_time = a_time;
          time = a_state[0];
          xd = a_state[1];
          x = a_state[2];

          // calculate rate
          calculate_rate();

          // return rate
          a_rate[0] = 1.0;
          a_rate[1] = xdd;
          a_rate[2] = xd;

          }

          void advance( const double end_time , const double time_step ){

          double a_time;
          state_type a_state;

          a_time = system_time;
          a_state = get_state();

          // https://stackoverflow.com/questions/10976078/using-boostnumericodeint-inside-the-class
          integrate_const( stepper , *this , a_state, a_time , end_time , time_step );

          system_time = end_time;
          set_state( a_state );
          calculate_rate();

          }

          }; // end class spring_type

          spring my_spring;

          // https://stackoverflow.com/questions/34181135/converting-an-stdmap-to-an-rcpplist

          // [[Rcpp::export]]
          void initialise_model( double t ){
          my_spring.initialise_model( t );
          }

         // [[Rcpp::export]]
          void pull_variables_from_model(){
          my_spring.pull_variables_from_model();
          }

          // [[Rcpp::export]]
          NumericVector get_spring_variables(){
          NumericVector my_vector;
          my_vector = my_spring.variable; // coerces unordered_map to named vector
          return my_vector;
          }

          // whenever I try to expose the advance method I get an error
          // // [[Rcpp::export]]
          // void advance( double t , double dt ){
          // my_spring.advance( t , dt );
          // }

          int main(){
          // do nothing
          }
          ') # end of sourceCpp

# instantiate and query the model
initialise_model(0)
pull_variables_from_model()
x <- get_spring_variables()
print(x) # this all works.

# I now want to do this
advance( 1 , 0.1)
pull_variables_from_model()
x <- get_spring_variables() # is this necessary?
print(x)

# end of R script

Thank you

Simon Woodward
Senior Scientist (Mathematical Modelling)
DairyNZ
Cnr Ruakura & Morrinsville Roads | Newstead | Private Bag 3221| Hamilton 3240 | NEW ZEALAND
Ph +64 7 858 3750 | DDI +64 7 858 3787 | Fax +64 7 858 3751
Web www.dairynz.co.nz<http://www.dairynz.co.nz/> | www.GoDairy.co.nz<http://www.godairy.co.nz/> | www.getfresh.co.nz<http://www.getfresh.co.nz/>


-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20180914/45071be9/attachment-0001.html>


More information about the Rcpp-devel mailing list