In this paragraph, a basic example demonstrates how to write the R code, translate it and call it from C++. Particular emphasis is placed on the C++ code. First of all, the R function is defined which accepts one argument called a, adds two to a and stores it into b. The variable b is returned at the end of the function. The R function called f is translated to an external pointer to the C++ function.
<- function(a) {
f <- a + 2
b return(b)
} library(ast2ast)
<- translate(f) f_cpp
The C++ function depends on RcppArmadillo and ast2ast therefore the required macros and headers were included. Moreover, ETR requires std=c++17 therefore the corresponding plugin is added. The function getXPtr is defined by the function RcppXPtrUtils::cppXPtr. In the last 5 lines, the translated code is depicted. The function f returns a sexp and gets one argument of type sexp called a. The body of the function looks almost identical to the R function. Except that the variable b is defined in the first line of the body with the type sexp. The function i2d converts an integer to a double variable. This is necessary since C++ would identify the 2 as an integer which is not what the user wants in this case.
// [[Rcpp::depends(ast2ast, RcppArmadillo)]]
#include <RcppArmadillo.h>
#include <Rcpp.h>
// [[Rcpp::plugins(cpp17)]]
using namespace Rcpp;
#include <etr.hpp>
// [[Rcpp::export]]
();
SEXP getXPtr
(sexp a) {
sexp f;
sexp b= a + i2d(2);
b return(b);
}
Afterwards, the translated R function has to be used in C++ code. This would be your package code for example. First, the macros were defined for RcppArmadillo and ast2ast. Subsequently, the necessary header files were included. As already mentioned ast2ast requires std=c++17 thus the required plugin is included. To use the function, it is necessary to dereference the pointer. The result of the dereferenced pointer has to be stored in a function pointer. Later the function pointer can be used to call the translated function. Therefore, a function pointer called fp is defined. It is critical that the signature of the function pointer matches the one of the translated function. Perhaps it would be a good idea to check the R function before it is translated. After defining the function pointer, a function is defined which is called later by the user (called call_package). This function accepts the external pointer. Within the function body, a variable f is defined of type fp and inp is assigned to it. Next, a sexp object called a is defined which stores a vector of length 3 containing 1, 2 and 3. The function coca is equivalent to the c function R. Afterwards a is printed. Followed by the call of the function f and storing the result in a. The variable a is printed again to show that the values are changed according to the code defined in the R function.
// [[Rcpp::depends(RcppArmadillo, ast2ast)]]
#include "etr.hpp"
// [[Rcpp::plugins("cpp17")]]
typedef sexp (*fp)(sexp a);
// [[Rcpp::export]]
void call_package(Rcpp::XPtr<fp> inp) {
= *inp;
fp f = coca(1, 2, 3);
sexp a (a);
print= f(a);
a ("a is now:");
print(a);
print}
The user can call now the package code and pass the R function to it. Thus, the user only has to install the compiler or Rtools depending on the operating system. But it is not necessary to write the function in Rcpp.
call_package(f_cpp)
## 1
## 2
## 3
## a is now:
## 3
## 4
## 5
In the last section, the usage of ast2ast was described. However, only sexp variables were defined. Which are most likely not used in your package. Therefore interfaces to common libraries are defined. First of all, ast2ast can communicate with Rcpp which alleviates working with the library substantially. The code below shows that it is possible to pass a sexp object to a variable of type NumericVector or NumericMatrix and vice versa. Here, the data is always copied.
// [[Rcpp::depends(ast2ast, RcppArmadillo)]]
#include <RcppArmadillo.h>
#include <Rcpp.h>
// [[Rcpp::plugins(cpp17)]]
using namespace Rcpp;
#include <etr.hpp>
// [[Rcpp::export]]
void fct() {
// NumericVector to sexp
{1, 2};
NumericVector aa_; // sexp a_ = a; Error!
sexp a_ = a;
(a_);
print
// sexp to NumericVector
b_ = coca(3, 4);
sexp = b_;
NumericVector b ::Rcout << b << std::endl;
Rcpp
// NumericMatrix to sexp
(3, 3);
NumericMatrix cc_; // sexp c_ = c; Error!
sexp c_ = c;
(c_);
print
// sexp to NumericMatrix
d_ = matrix(colon(1, 16), 4, 4);
sexp = d_;
NumericMatrix d ::Rcout << d << std::endl;
Rcpp}
<- fct() trash
## 1
## 2
## 3 4
## 0 0 0
## 0 0 0
## 0 0 0
## 1.00000 5.00000 9.00000 13.0000
## 2.00000 6.00000 10.0000 14.0000
## 3.00000 7.00000 11.0000 15.0000
## 4.00000 8.00000 12.0000 16.0000
Besides Rcpp types, sexp objects can transfer data to RcppArmadillo objects and it is also possible to copy the data from RcppArmadillo types to sexp objects using the operator =. The code below shows that it is possible to pass a sexp object to a variable of type vec or mat and vice versa. Here the data is always copied.
// [[Rcpp::depends(ast2ast, RcppArmadillo)]]
#include <RcppArmadillo.h>
#include <Rcpp.h>
// [[Rcpp::plugins(cpp17)]]
using namespace arma;
#include <etr.hpp>
// [[Rcpp::export]]
void fct() {
// vec to sexp
::vec a(4, fill::value(30.0));
armaa_; // sexp a_ = a; Error!
sexp a_ = a;
(a_);
print
// sexp to vec
b_ = coca(3, 4);
sexp = b_;
vec b .print();
b
// mat to sexp
(3, 3, fill::value(31.0));
mat cc_; // sexp c_ = c; Error!
sexp c_ = c;
(c_);
print
// sexp to mat
d_ = matrix(colon(1, 16), 4, 4);
sexp = d_;
mat d .print();
d}
<- fct() trash
## 30
## 30
## 30
## 30
## 3.0000
## 4.0000
## 31 31 31
## 31 31 31
## 31 31 31
## 1.0000 5.0000 9.0000 13.0000
## 2.0000 6.0000 10.0000 14.0000
## 3.0000 7.0000 11.0000 15.0000
## 4.0000 8.0000 12.0000 16.0000
You can pass the information of data stored on heap to a sexp object. The constructor for type vector accepts 3 arguments:
The constructor for the type matrix accepts 4 arguments:
If cob is 0 then the data is copied. Else if cob is 1 then the pointer itself is copied. Meaning that the ownership is transferred to the sexp object and the user should not call delete [] on the pointer. Be aware that only one sexp variable can take ownership of one vector otherwise the memory is double freed. Else if cob is 2 the ownership of the pointer is only borrowed. Meaning that the sexp object cannot be resized. The user is responsible for freeing the memory! The code below shows how the pointer interface works in general. Showing how sexp objects can be created by passing the information of pointers (double*) which hold data on the heap. Currently, only constructors are written which can use the pointer interface.
// [[Rcpp::depends(ast2ast, RcppArmadillo)]]
#include <RcppArmadillo.h>
#include <Rcpp.h>
// [[Rcpp::plugins(cpp17)]]
using namespace arma;
#include <etr.hpp>
// [[Rcpp::export]]
void fct() {
int size = 3;
// copy
double* ptr1;
= new double[size];
ptr1 int cob = 0;
(size, ptr1, cob);
sexp adelete [] ptr1;
= vector(3.14, 5);
a (a);
print
();
print
// take ownership
double* ptr2;
= new double[size];
ptr2 = 1;
cob (size, ptr2, cob);
sexp b= vector(5, 3);
b (b);
print
();
print
// borrow ownership
double* ptr3;
= new double[size];
ptr3 = 2;
cob (size, ptr3, cob);
sexp c//error calls resize
//c = vector(5, size + 1);
= vector(4, size);
c (c);
print
();
print(size, ptr3, cob);
sexp d= d + 10;
d (d);
print();
print
delete[] ptr3;
}
<- fct() trash
## 3.14
## 3.14
## 3.14
## 3.14
## 3.14
##
## 5
## 5
## 5
##
## 4
## 4
## 4
##
## 14
## 14
## 14
The pointer interface is particularly useful if the user function has to change the data of a vector or matrix of type NumericVector, vec, NumericMatrix or mat. Assuming that the user passes a function that accepts its arguments by reference it is easy to modify any variable which has a type that can return a pointer to its data. In the code below it is shown how sexp objects are constructed using the pointer interface. Thereby changing the content of variables which has an Rcpp type, a RcppArmadillo type, or is of type std::vector
// [[Rcpp::depends(ast2ast, RcppArmadillo)]]
#include <RcppArmadillo.h>
#include <Rcpp.h>
// [[Rcpp::plugins(cpp17)]]
using namespace Rcpp;
using namespace arma;
#include <etr.hpp>
typedef void (*fp)(sexp& a);
// [[Rcpp::export]]
void call_package(Rcpp::XPtr<fp> inp) {
= *inp;
fp f
// NumericVector
{1, 2, 3};
NumericVector a_rcpp(a_rcpp.size(), a_rcpp.begin(), 2);
sexp a(a);
f::Rcout << a_rcpp << std::endl;
Rcpp
//arma::vec
(2, fill::value(30));
vec a_arma(2, a_arma.memptr(), 2);
sexp b(b);
f.print();
a_arma
// NumericMatrix
(2, 2);
NumericMatrix c_rcpp(2, 2, c_rcpp.begin(), 2);
sexp c(c);
f::Rcout << c_rcpp << std::endl;
Rcpp
//arma::mat
(3, 2, fill::value(30));
mat d_arma(3, 2, d_arma.memptr(), 2);
sexp d(d);
f.print();
d_arma}
<- function(a) {
f <- a + 2
a
}
library(ast2ast)
<- translate(f, reference = TRUE)
fa2a <- call_package(fa2a) trash
## 3 4 5
## 32.0000
## 32.0000
## 2.00000 2.00000
## 2.00000 2.00000
##
## 32.0000 32.0000
## 32.0000 32.0000
## 32.0000 32.0000
In this section two examples are shown to illustrate how ast2ast could be used in R packages. First, it is shown how ODE-systems can be solved very efficiently. In order to do this the R package r2sundials is used. The second examples point out how ast2ast can be used together with the R package paropt to optimize parameters of ODE-systems. Both examples only construct wrapper functions that are used by the package code. Probably it would be more efficient when the package code itself would take care of the transfer of data between Rcpp/RcppArmadillo and ast2ast.
In this example it is shown how ODE-systems can be solved by using r2sundials and ast2ast. The code below shows how r2sundials is used normally. Either using an R function or an external pointer to a C++ function.
library(Rcpp)
library(ast2ast)
library(r2sundials)
## Loading required package: rmumps
library(RcppXPtrUtils)
library(microbenchmark)
# R version
<- seq(0, 5, length.out=101)
ti <- list(a = 2)
p <- c(nu = 2, a = 1)
p <- 0
y0 <- function(t, y, p, psens) {
frhs -p["nu"]*(y-p["a"])
}
<- r2cvodes(y0, ti,
res_exp param = p)
frhs, attributes(res_exp) <- NULL
# External pointer
<- cppXPtr(code = '
ptr_exp int rhs_exp(double t, const vec &y,
vec &ydot,
RObject ¶m,
NumericVector &psens) {
double a = 1;
double nu = 2;
ydot[0] = -nu*(y[0] - a);
return(CV_SUCCESS);
}
',
depends=c("RcppArmadillo",
"r2sundials","rmumps"),
includes="using namespace arma;\n
#include <r2sundials.h>",
cacheDir="lib", verbose=FALSE)
<- c(a = 1)
pv <- r2cvodes(y0, ti,
res_exp2 param = pv)
ptr_exp, attributes(res_exp2) <- NULL
In the code below is the wrapper function defined which is later called by r2sundials. This function is called rhs_exp_wrapper and has the correct function signature. Furthermore, a global function pointer named Fct is defined of type void (*user_fct) (sexp& y_, sexp& ydot_). Within rhs_exp_wrapper* the data of the vectors y and ydot are used to construct two sexp objects which are passed to Fct. Thus, the vector ydot is modified by the function passed from the user. Furthermore, another function called solve_ode is defined. Which calls the code from r2sundials, solves the ODE-system, and returns the output to R. In R the user defines the R function ode. Next, the function is translated and passed to solve_ode. Comparing the results shows that all three approaches (R, C++, ast2ast) generate the same result. Afterwards a benchmark is conducted showing that R is substantially slower than C++ and that the translated function is almost as fast as C++. Mentionable, it is possible to increase the speed of the ast2ast version by using the at function and the *_db* extension.
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::depends(rmumps)]]
// [[Rcpp::depends(r2sundials)]]
// [[Rcpp::depends(ast2ast)]]
// [[Rcpp::plugins("cpp17")]]
#include "etr.hpp"
#include "RcppArmadillo.h"
#include "r2sundials.h"
using namespace arma;
using namespace Rcpp;
typedef int (*fp)(double t, const vec &y,
&ydot, RObject ¶m,
vec &psens);
NumericVector
typedef void (*user_fct)(sexp& y_,
& ydot_);
sexp;
user_fct Fct
int rhs_exp_wrapper(double t, const vec &y,
&ydot, RObject ¶m,
vec &psens) {
NumericVector (param);
NumericVector pconst int size = y.size();
ydot_(size, ydot.memptr(), 2);
sexp
double* ptr = const_cast<double*>(
.memptr());
yy_(size, ptr, 2);
sexp (y_, ydot_);
Fctreturn(CV_SUCCESS);
}
// [[Rcpp::export]]
(XPtr<user_fct> inp,
NumericVector solve_ode,
NumericVector time) {
NumericVector y= *inp;
Fct <fp> ptr = XPtr<fp>(new fp(
XPtr&rhs_exp_wrapper));
=
Environment pkg ::namespace_env("r2sundials");
Environment= pkg["r2cvodes"];
Function solve = solve(y, time,
NumericVector output , time);
ptr
return output;
}
# ast2ast version
<- seq(0, 5, length.out=101)
ti <- 0
y0
library(ast2ast)
<- function(y, ydot) {
ode <- 2
nu <- 1
a 1] <- -nu*(y[1] - a)
ydot[
}<- translate(ode,
pointer_to_ode reference = TRUE)
<- solve_ode(pointer_to_ode,
res_exp3
ti, y0)attributes(res_exp3) <- NULL
stopifnot(identical(res_exp,
res_exp2,
res_exp3))<- microbenchmark(
out r2cvodes(y0, ti,
param = p),
frhs, r2cvodes(y0, ti,
param = pv),
ptr_exp, solve_ode(pointer_to_ode,
ti, y0))
boxplot(out, names=c("R", "C++", "ast2ast"))
In the second example, it is shown how ast2ast is used together with paropt. Again the code below consists of an Rcpp part and a part written in R. Within Rcpp the C++ code is defined which describes the Lotka-Volterra model (ode_system). For details check the documentation of paropt. Subsequently, a function is defined that returns an external pointer to the ODE-system (test_optimization). As already described in the last example a wrapper function is created which calls the actual ODE-system (ode_system_wrapper). Furthermore, a function is defined which calls the code from the package paropt called optimize_paropt and conducts the optimization. Within R the ODE-system is defined, translated and used for optimization. Both versions the pure C++ function and the translated R code yield the same result. The C++ function is almost as fast as the translated R code. Interesting is the fact that paropt calls the ODE-solver in parallel. This example demonstrates that the translated function can be called in parallel. Which is not possible with R functions. Admittedly, the printing is conducted by Rcpp. Therefore it is unknown how printing behaves in parallel.
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::depends(ast2ast)]]
// [[Rcpp::plugins("cpp17")]]
// [[Rcpp::depends(paropt)]]
#include "etr.hpp"
#include "RcppArmadillo.h"
using namespace Rcpp;
// paropt version
typedef int (*OS)(double &t,
std::vector<double> ¶ms,
std::vector<double> &states);
int ode_system(double &t,
std::vector<double> ¶ms,
std::vector<double> & states) {
double a = params[0];
double b = params[1];
double c = params[2];
double d = params[3];
double n1 = states[0];
double n2 = states[1];
[0] = n1*c*n2 - n1*d;
states[1] = n2*a - n2*b*n1;
states
return 0;
}
// [[Rcpp::export]]
<OS> test_optimization() {
XPtr<OS> xpfun = XPtr<OS>(new OS(
XPtr&ode_system));
return xpfun;
}
// ast2ast version
typedef int (*paropt_fct)(double &t,
std::vector<double> ¶ms,
std::vector<double> & states);
typedef void (*user_fct2)(sexp& p,
& y);
sexp
;
user_fct2 Fct_paropt
int ode_system_wrapper(
double &t,
std::vector<double> ¶ms,
std::vector<double> & states) {
(params.size(), params.data(), 2);
sexp p(states.size(), states.data(), 2);
sexp s(p, s);
Fct_paroptreturn 0;
}
// [[Rcpp::export]]
(XPtr<user_fct2> inp,
List optimize_paropt,
NumericVector time,
DataFrame lb,
DataFrame ub) {
DataFrame states= *inp;
Fct_paropt <paropt_fct> ptr = XPtr<paropt_fct>(
XPtrnew paropt_fct(&ode_system_wrapper));
=
Environment pkg ::namespace_env("paropt");
Environment=
Function optim ["optimizer_pointer"];
pkg
{1e-8, 1e-8};
NumericVector abs_tols
= optim(time, ptr, 1e-6,
List output , lb, ub,
abs_tols, 40,
states1000, 0.0001,
"bdf");
return output;
}
#states
<- system.file("examples",
path package = "paropt")
<- read.table(paste(
states
path,"/states_LV.txt",
sep = ""),
header = T)
# parameter
<- data.frame(time = 0,
lb a = 0.8,
b = 0.3,
c = 0.09,
d = 0.09)
<- data.frame(time = 0,
ub a = 1.3,
b = 0.7,
c = 0.4,
d = 0.7)
suppressMessages(library(paropt))
set.seed(1)
<- Sys.time()
start_time <- optimizer_pointer(
df_cpp integration_times = states$time,
ode_sys = test_optimization(),
relative_tolerance = 1e-6,
absolute_tolerances = c(1e-8, 1e-8),
lower = lb, upper = ub, states = states,
npop = 40, ngen = 1000, error = 0.0001,
solvertype = "bdf")
<- Sys.time()
end_time <- end_time - start_time
cpp_time
# ast2ast with at and _db
<- function(params, states) {
ode = at(params, 1)
a_db = at(params, 2)
b_db = at(params, 3)
c_db = at(params, 4)
d_db = at(states, 1)
n1_db = at(states, 2)
n2_db at(states, 1) = n1_db*c_db*n2_db -
*d_db;
n1_dbat(states, 2) = n2_db*a_db -
*b_db*n1_db;
n2_db
}
<- ast2ast::translate(
pointer_to_ode
ode,reference = TRUE)
set.seed(1)
<- Sys.time()
start_time <- optimize_paropt(pointer_to_ode,
df_ast2ast $time,
states
lb, ub, states)<- Sys.time()
end_time <- end_time - start_time
a2a_time
stopifnot(identical(df_cpp[[8]],
8]]) )
df_ast2ast[[
<- data.frame(cpp = cpp_time,
times ast2ast = a2a_time)
::kbl(times) kableExtra
cpp | ast2ast |
---|---|
7.576208 secs | 10.31115 secs |