Define a function within another function in C++

Question

To solve differential equations, I use the Rcpp package in R and the boost package in C++. I created the eq function in C++, which transform the R function mod_cpp into a "C++ function" (see below). The eq function is next put as an argument in the integrate_const function. Finally, I can compile it in R with Rcpp and I get an R function called my_fun, which only depends on the vector vs. Everything works well.

// [[Rcpp::export]]
void my_fun22(Rcpp::NumericVector &x, const double t){
  Function f("mod_cpp");
  x=f(_["t"]=t,_["x"]=x);
}     

Rcpp::NumericVector nvec(130);
void eq(const state_type &x, state_type &dxdt, const double t){
  boost_array_to_nvec2(x, nvec);
  my_fun22(nvec,t);
  nvec_to_boost_array2(nvec, dxdt);
}

Rcpp::NumericMatrix my_fun(const Rcpp::NumericVector vs) {
  state_type x = nvec_to_boost_array(vs); // initial conditions
  integrate_const(make_dense_output( 1E-4 , 1E-4 , stepper_type () ) ,
                eq , x , 0.0 , 120.0 , 1.0 , write_cout);
  return data;
}

The problem is that all the parameters of my model contained in the mod_cpp function are fixed. What I want to do now is to create another function that does the same job as my_fun but that depends on some parameters of the model. More precisely, create a function called my_fun2(vs, theta) that depends on vs AND theta. I already tried to do so, but was struggling with redifining the eq function within the my_fun function as it is not allowed to define a function within another function in R. Did I miss something?

Answer 1

For achieving what you want I would use c++ 11 lambda functions and instead of creating eq as a function I would define it as a second order function: a function that returns a different function which argument can be the theta you mentioned. For instance:

...
Rcpp::NumericMatrix my_fun(const Rcpp::NumericVector vs, const float theta) {
  state_type x = nvec_to_boost_array(vs); // initial conditions
  integrate_const(make_dense_output( 1E-4 , 1E-4 , stepper_type () ) ,
            eq(theta) , x , 0.0 , 120.0 , 1.0 , write_cout);
...
}

and then the modified eq would be something like:

std::function<void(const state_type, state_type, const double)> eq(const float theta) {
   return [&theta](const state_type &x, state_type &dxdt, const double t) {
     // do something with theta here
     // for instance, modify nvec
     ...
     // then your old function body
     boost_array_to_nvec2(x, nvec);
     my_fun22(nvec,t);
     nvec_to_boost_array2(nvec, dxdt);
   }
}

This way you can pass parameters in execution time. I haven't tested it, but the idea should be fine I guess.

To use c++ 11 features you need to make sure you have installed a c++ compiler in your system path that can compile c++ 11 code and to edit the Makevars file in your R project and add the following at the top:

PKG_CXXFLAGS = -g -O3 -std=c++11

More about this here, for example.