---
title: "Rcpp"
author: "JJB + course"
output:
   html_document:
     toc: true
     toc_float:
       collapsed: false
---

# Rcpp

## Example: Saying hello

```{Rcpp hello_world_ex}
#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
void hello_world_cpp() {
  Rcout << "Hello R/C++ World!\n";  
}
```

## Example: Embedding Rcpp into RMarkdown.

To use _Rcpp_ code within _RMarkdown, we change `{R}` to `{Rcpp}` within the code chunk block.

## Example: Calling the C++ Function within *R*

```{r}
# Call C++ Code like a normal R function
hello_world_cpp()    
```

# Practical Usage

## Example: Embedding `is_odd_cpp()` in _R_

```{Rcpp}
#include<Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
bool is_odd_cpp(int n = 10) {// function definition
  bool v = (n % 2  == 1);    // modulus check
  return v;                  // return result
}
```

## Example: Calling `is_odd_cpp()` in _R_

```{r}
is_odd_cpp()

is_odd_cpp(10)

is_odd_cpp(3)
```

## Example: Vectorizing Is Odd in C++

`int` can only hold one value. What if we wanted to process a vector?

```{Rcpp}
#include<Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
LogicalVector is_odd_vec_cpp(IntegerVector x) {
  int n = x.size();            // number of elements
  LogicalVector v(n);          // logical vector
  for(int i = 0; i < n; i++) { // process each value
    v[i] = (x[i] % 2  == 1);   // modulus check
  }
  return v;
}
```

Differences:

```{r, eval = FALSE}
is_odd_cpp(1:5)
# Error in is_odd_cpp(1:5) : Expecting a single value: [extent=5].
```

```{r}
is_odd_vec_cpp(1:5)
```

## Example: Loops

Consider the mean:

```{r}
mean_r = function(x) {
  sum_r = 0
  for (i in seq_along(x)) {
    sum_r = sum_r + x[i]
  }
  
  sum_r / length(x)
}
```

**C++** Implementation

```{Rcpp mean}
#include <Rcpp.h>
using namespace Rcpp;           // Import Statements

// [[Rcpp::export]]
double mean_cpp(NumericVector x) { // Declaration
  
  int n = x.size();              // Find the vector length
  double sum_x = 0;              // Set up storage

  for(int i = 0; i < n; ++i) {   // For Loop in C++
    // Shorthand for sum_x = sum_x + x[i];
    sum_x += x[i];   
  }
  return sum_x / n;             // Return division
}
```

Verifying equality:

```{r}
# Done in *R*

set.seed(112)                # Set seed

x = rnorm(10)                # Generate data

all.equal(mean_cpp(x), mean_r(x)) # Test Functions
```


### Exercise: Creating a scalar version of `is_divisible_by()`

Modify the `is_odd_cpp()` to `is_divisible_by()`, which checks to see if a number is divisible by another.

```{Rcpp}
#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
bool is_divisible_by(int num, int divisor) {
  
  bool result = (num % divisor == 0);
  
  return result;
}
```

```{r}
is_divisible_by
```


```{r}
is_divisible_by(4, 2)

is_divisible_by(5, 2)

is_divisible_by(5, 3)
```



## Example: Squaring Numbers

The squaring operation switches from `base^exp` to `pow(base, exp);`

```{Rcpp}
#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
NumericVector square_numbers(NumericVector x) {
  
  int n = x.size();
  NumericVector res(n);
  
  for(int i = 0; i < n; i++) {
    res[i] = pow(x[i], 2.0); // x^2
  }
  
  return res;
}
```

### Exercise: Compute the Sum of Squares

Modify the squaring operation given previously such that is computes
the sum of squares for $x$. 

\[\sum\limits_{i = 1}^n {{{\left( {{x_i} - \bar x} \right)}^2}} \]

**Recall:** C++ indices start at 0 **NOT** 1. 


```{Rcpp}

#include <Rcpp.h>
using namespace Rcpp;           // Import Statements

// [[Rcpp::export]]
double mean_cpp(NumericVector x) { // Declaration
  
  int n = x.size();              // Find the vector length
  double sum_x = 0;              // Set up storage

  for(int i = 0; i < n; ++i) {   // For Loop in C++
    // Shorthand for sum_x = sum_x + x[i];
    sum_x += x[i];   
  }
  return sum_x / n;             // Return division
}

// [[Rcpp::export]]
double sum_of_squares(NumericVector x){
  
  int n = x.size();
  
  double summed_value = 0;
  double x_mean = mean_cpp(x);
  
  for(int i = 0; i< n; i++) {
    summed_value = summed_value + pow((x[i] - x_mean), 2.0);
  }
  
  
  return summed_value;
}
```


# Behind the Scenes: Proxy Model

## Example: 

```{Rcpp}
// Written in C++

// Import Statements
#include <Rcpp.h>
using namespace Rcpp; 

// Assume numeric type
// [[Rcpp::export]]
void ref_update_num(NumericVector x) { 
    x = x + 1; 
} 

// Assume integer type
// [[Rcpp::export]]
void ref_update_int(IntegerVector x) { 
    x = x + 1; 
}
```


```{r}
# In R
x = c(-8L, 2L)

# Note: Neither function will output since there is no return

ref_update_num(x)  # Call numeric update
x                  # Same values!!! 

ref_update_int(x)  # Call numeric update
x                  # Different values!!!
```


## Example: Type Sensitivity

```{r}
add_r = function(a, b) {
  return(a + b)
}

add_r(0L, 2L)   # Remember L means integer!
add_r(2.5, 1.1) # Double/numeric
```

```{Rcpp, eval = TRUE}
#include <Rcpp.h>
using namespace Rcpp;           // Import Statements

// [[Rcpp::export]]
double add_cpp(double a, double b) { // Declaration
  double out = a + b; // Add `a` to `b`
  return out;         // Return output
}
```

```{r}
add_cpp(0L, 2L)   # Integers into double
add_cpp(2.5, 1.1) # Double into double
```

### Example: Calling a C++ Function with different types

```{Rcpp}
#include <Rcpp.h>
using namespace Rcpp;           // Import Statements

// [[Rcpp::export]]
double add_cpp(double a, double b) { // Declaration
  double out = a + b;                // Add `a` to `b`
  return out;                        // Return output
}

// [[Rcpp::export]]
int add_cpp_int(int a, int b) { // Declaration
  return add_cpp(a, b);         // Call previous function
}
```

```{r}
add_cpp_int(2.5, 1.1)  # Call in *R*
```