With the Rcpp package it is possible to write, compile and call C++ code on the fly during an R session using the cppFunction() and sourceCpp() functions. A plugin has been written that allows the C++ components of FLasher to be used during an R session, including access to all of the FLCpp classes (the C++ implementation of the FLR classes) and automatic differentiation (AD) functionality through access to the CppAD library.
Here we demonstrate how the Rcpp functions cppFunction() and sourceCpp() can be used
cppFunction() is used for writing functions in C++ that you want to call from R. You write your C++ function using standard C++ types for the arguments and returned object and the automatic Rcpp as<> and wrap takes care of the conversion. The C++ function is passed as a string to cppFunction() during the R session:
cppFunction("
int my_add(int x, int y, int z) {
int sum = x + y + z;
return sum;
}")The C++ function can then be called as if it was an R function:
my_add(1L, 2L, 10L)## [1] 13It is possible to use C++11 functionality, for example, using range-based loops and auto types: To do this we need to use C++11 plugin. This function takes a vector of integers and increments each value in the vector.
cppFunction("
std::vector<int> rbl_demo(std::vector<int> v){
for (auto& i : v){
i++;
}
return v;
}",
plugins = "cpp11")We can call it as if it was a normal R function.
rbl_demo(c(1L, 2L, 3L))## [1] 2 3 4sourceCpp is for writing longer C++ scripts and can contain multiple functions and classes, not all of which need to be exposed to R. Exposing the desired functions to R is done using the Rcpp::attributes (see the vignette in the Rcpp package for details). The C++ code can either be included as a text string or written in a separate file. Writing the code in a separate file makes it easier to manage and also your text editor will highlight the syntax correctly. You need to include the include to get all the advantages of Rcpp. Ideally, the following source code should be in a separate script. However, for the purposes of this vignette we write the C++ code as a text string, save it as a temporary file and then source the file. Be careful that the #include line does not get interpreted as a comment by R! This is why it is not on a separate line.
source_code <- " #include <Rcpp.h>
// This function is not exposed to R
double foo(double x){
return 2.0 * x;
}
// This function is exposed to R and calls the unexposed one
// [[Rcpp::export]]
double call_foo(double x){
double y = foo(x);
return y;
}
"
cat(source_code, file = paste(tempdir(), "test-1.cpp", sep = "/"))
sourceCpp(file = paste(tempdir(), "test-1.cpp", sep = "/"))call_foo(3.5)## [1] 7C++11 code can be included using the C++11 plugin:
source_code <- " #include <Rcpp.h>
// [[Rcpp::plugins(cpp11)]]
// [[Rcpp::export]]
std::vector<double> rbl_demo2(std::vector<double> v){
for (auto& i : v){
i = i * 2.0;
}
return v;
}
"
cat(source_code, file = paste(tempdir(), "test-2.cpp", sep = "/"))
sourceCpp(file = paste(tempdir(), "test-2.cpp", sep = "/"))rbl_demo2(c(1.3, 2.6, 3.9))## [1] 2.6 5.2 7.8Using the FLasher plugin means that you have access to the methods and classes in the C++ code of the FLasher package. For example, you can pass in and manipulate FLQuant objects. In this example, we write a C++ function that takes two FLQuants adds them together and returns the resulting FLQuant.
To use it with cppFunction() you must specify it as a depends argument:
cppFunction("
FLQuant calc_catches(FLQuant landings, FLQuant discards){
FLQuant catches = landings + discards;
return catches;
}",
depends = "FLasher")data(ple4)
landings <- landings.n(ple4)[, ac(2000:2003)]
discards <- discards.n(ple4)[, ac(2000:2003)]The C++ function can be called as normal:
calc_catches(landings, discards)## An object of class "FLQuant"
## , , unit = unique, season = all, area = unique
##
## year
## age 2000 2001 2002 2003
## 1 128525.79 97495.00 273411.66 91796.79
## 2 138479.50 187737.70 172176.50 531815.50
## 3 147232.10 142364.40 199413.20 136174.20
## 4 242206.10 92615.80 77923.40 74693.40
## 5 32822.43 109946.10 41589.14 34002.05
## 6 12638.84 12252.56 47176.75 20056.13
## 7 3078.29 4697.26 5248.95 19344.91
## 8 1433.52 1455.38 2162.48 1819.76
## 9 1209.85 716.29 593.44 673.88
## 10 2395.19 1899.89 1503.05 1324.69
##
## units: 1000To use the FLasher plugin with sourceCpp() you must add a depends at the top of the script and include the FLasher header file. Again, be careful that the #include line does not interpreted as a comment by R. For this reason we place it on the same line as another line but include the line separator . This is not necessary if creating a stand alone C++ file from scratching instead of trying to create a text string to write to a file.
source_code <- "
// [[Rcpp::depends(FLasher)]] \n #include <FLasher.h>
// [[Rcpp::export]]
FLQuant calc_catches2(FLQuant landings, FLQuant discards){
FLQuant catches = landings + discards;
return catches;
}
"
cat(source_code, file = paste(tempdir(), "test-3.cpp", sep = "/"))
sourceCpp(file = paste(tempdir(), "test-3.cpp", sep = "/"))calc_catches2(landings, discards)## An object of class "FLQuant"
## , , unit = unique, season = all, area = unique
##
## year
## age 2000 2001 2002 2003
## 1 128525.79 97495.00 273411.66 91796.79
## 2 138479.50 187737.70 172176.50 531815.50
## 3 147232.10 142364.40 199413.20 136174.20
## 4 242206.10 92615.80 77923.40 74693.40
## 5 32822.43 109946.10 41589.14 34002.05
## 6 12638.84 12252.56 47176.75 20056.13
## 7 3078.29 4697.26 5248.95 19344.91
## 8 1433.52 1455.38 2162.48 1819.76
## 9 1209.85 716.29 593.44 673.88
## 10 2395.19 1899.89 1503.05 1324.69
##
## units: 1000As well as providing access to the FLCppad classes and methods, the plugin allows the AD library CppAD that FLasher uses to be accessed. Unfortunately, at the moment, the interface is a bit clunky.
Here we write C++ code that returns the value and the gradient of the banana function (see the R help page for optim for more information on the banana function). We can pass the exposed gradient function to R’s optim functions. There is also an exposed function that returns the Hessian.
The function func() can be rewritten to be any function that you want that derivatives for. The rest of the code remains the same (it would be good to have this other code in the package but it is not possible at the moment).
source_code <- "
// [[Rcpp::depends(FLasher)]] \n #include <FLasher.h>
// This is the function we want to solve - the banana function
// It is templated because we need versions of it that deal with
// types double (for normal evaluation) and adouble (for AD evaluation)
template <typename T>
std::vector<T> func(std::vector<T> params){
std::vector<T> res(1, 0.0);
res[0] = 100 * pow((params[1] - params[0] * params[0]), 2.0) + pow((1 - params[0]), 2.0);
return res;
}
// Evaluates the function
// [[Rcpp::export]]
std::vector<double> eval_function(std::vector<double> params){
return func(params);
}
// Uses CppAD magic to get the gradient of the function
// [[Rcpp::export]]
std::vector<double> eval_gradient(std::vector<double> params){
std::vector<adouble> x(params.begin(), params.end());
CppAD::Independent(x);
std::vector<adouble> res = func(x);
CppAD::ADFun<double> fun(x, res);
return fun.Jacobian(params);
}
// Uses CppAD magic to get the Hessian
// [[Rcpp::export]]
std::vector<double> eval_hessian(std::vector<double> params, unsigned int var = 0){
std::vector<adouble> x(params.begin(), params.end());
CppAD::Independent(x);
std::vector<adouble> res = func(x);
CppAD::ADFun<double> fun(x, res);
return fun.Hessian(params, var);
}
"
cat(source_code, file = paste(tempdir(), "test-4.cpp", sep = "/"))
sourceCpp(file = paste(tempdir(), "test-4.cpp", sep = "/"))We can test this by solving the function in R with optim() using an approximate gradient, the exact gradient function and the AD gradient.
# Rosenbrock Banana function
fr <- function(x) {
100 * (x[2] - x[1] * x[1])^2 + (1 - x[1])^2
}
# The exact gradient of the banana function Gradient of 'fr'
grr <- function(x) {
x1 <- x[1]
x2 <- x[2]
c(-400 * x[1] * (x[2] - x[1] * x[1]) - 2 * (1 - x[1]), 200 *
(x[2] - x[1] * x[1]))
}We then solve the function using optim using the three methods for calculating the gradient:
Using the approximate gradient in R.
res1 <- optim(c(-1.2, 1), fr, method = "BFGS")
res1[c("par", "value", "counts")]## $par
## [1] 0.9998 0.9996
##
## $value
## [1] 3.827e-08
##
## $counts
## function gradient
## 118 38Using the exact gradient in R.
res2 <- optim(c(-1.2, 1), fr, grr, method = "BFGS")
res2[c("par", "value", "counts")]## $par
## [1] 1 1
##
## $value
## [1] 9.595e-18
##
## $counts
## function gradient
## 110 43Using the the AD gradient we get from using the CppAD library
res3 <- optim(c(-1.2, 1), eval_function, eval_gradient, method = "BFGS")
res3[c("par", "value", "counts")]## $par
## [1] 1 1
##
## $value
## [1] 9.595e-18
##
## $counts
## function gradient
## 110 43The version with the AD gradient is exactly the same as the version with the exact gradient function.
We can also get the Hessian:
# Estimated by R
optimHess(res1$par, fr)## [,1] [,2]
## [1,] 801.7 -399.9
## [2,] -399.9 200.0# Estimated by R using the the gradient function
optimHess(res2$par, fr, grr)## [,1] [,2]
## [1,] 802 -400
## [2,] -400 200# Calculated using the AD function
eval_hessian(res3$par)## [1] 802 -400 -400 200The above C++ code can be used to provide the gradients and Hessians for any functions. All the user needs to do is write their own func() function (with the same arguments).