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/>
_______________________________________________
Rcpp-devel mailing list
Rcpp-devel@lists.r-forge.r-project.org
https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
