Using_the_FLasher_Plugin.Rmd

---
title: "Using the FLasher Rcpp plugin for Automatic Differentiation and other nerdy stuff"
author: Finlay Scott, Iago Mosqueira - European Commission Joint Research Center
date: "`r format(Sys.time(), '%d %B, %Y')`"
output:
  html_document:
    fig_caption: yes
    number_sections: yes
vignette: >
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteIndexEntry{FLasher}
tags: [FLR FLasher forecast fisheries]
abstract:
license: European Union Public Licence (EUPL) V.1.1
---

```{r, pkgs, echo=FALSE, message=FALSE, warning=FALSE}
library(Rcpp)
library(FLasher)
library(ggplotFL)
library(knitr)
opts_chunk$set(dev='png', cache=FALSE, fig.width=4, fig.height=4, tidy=TRUE, dpi=72)
options(width=60)
```

# Introduction

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.

# Using *cppFunction()* and *sourceCpp()*

Here we demonstrate how the **Rcpp** functions *cppFunction()* and *sourceCpp()* can be used

## *cppFunction()*

*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:

```{r}
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:

```{r}
my_add(1L, 2L, 10L)
```

It 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.

```{r, warning=FALSE}
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.

```{r}
rbl_demo(c(1L,2L,3L))
```

## *sourceCpp()*

*sourceCpp* 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.

```{r}
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="/"))
```

```{r}
call_foo(3.5)
```

C++11 code can be included using the C++11 plugin:

```{r, warning=FALSE}
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="/"))
```

```{r}
rbl_demo2(c(1.3, 2.6, 3.9))
```
# Using the **FLasher** plugin

## With *cppFunction()*

Using 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 *FLQuant*s adds them together and returns the resulting *FLQuant*.

To use it with *cppFunction()* you must specify it as a *depends* argument:

```{r}
cppFunction('
FLQuant calc_catches(FLQuant landings, FLQuant discards){
    FLQuant catches = landings + discards;
    return catches;
    }',
    depends="FLasher"
)
```

```{r}
data(ple4)
landings <- landings.n(ple4)[,ac(2000:2003)]
discards <- discards.n(ple4)[,ac(2000:2003)]
```

The C++ function can be called as normal:

```{r}
calc_catches(landings, discards)
```

## With *sourceCpp()*

To 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 _\n_.
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.

```{r}
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="/"))
```

```{r, "demo"}
calc_catches2(landings, discards)
```

# Using automatic differentiation

As 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).

```{r, adexample}
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.

```{r}
# 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
grr <- function(x) { ## Gradient of 'fr'
    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.
```{r}
res1 <- optim(c(-1.2,1), fr, method = "BFGS")
res1[c("par", "value", "counts")]
```

Using the exact gradient in R.
```{r}
res2 <- optim(c(-1.2,1), fr, grr, method = "BFGS")
res2[c("par", "value", "counts")]
```

Using the the AD gradient we get from using the CppAD library 
```{r}
res3 <- optim(c(-1.2,1), eval_function, eval_gradient, method = "BFGS")
res3[c("par", "value", "counts")]
```

The version with the AD gradient is exactly the same as the version with the exact gradient function.

We can also get the Hessian:

```{r}
# Estimated by R
optimHess(res1$par, fr)
# Estimated by R using the the gradient function
optimHess(res2$par, fr, grr)
# Calculated using the AD function
eval_hessian(res3$par)
```

The 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).