[Rcpp]: Code Snippets

README.md

Rcpp Code Snippets

Learning Rcpp to speed things up!

simple_comparison.R

Example: Fibonacci Recursive expr min lq median uq max neval 1 fibR(n_fib) 3017.132 3200.2715 3354.8345 3562.343 5865.408 500 2 fibC(n_fib) 5.839 6.6255 8.7405 15.455 136.292 500 Speed Up = 383.8264x

Example: Fibonacci For Loop expr min lq median uq max neval 1 fibRiter(n_fib) 13.131 13.844 14.9400 15.430 47.256 500 2 fibCiter(n_fib) 1.349 1.666 1.9175 2.931 29.696 500 Speed Up = 7.791395x

Example: VAR Simulation expr min lq median uq max neval 1 rSim(a, u) 33539.167 35417.579 36168.264 37027.9425 77904.367 500 2 rcppSim(a, u) 164.803 173.787 187.624 204.4735 352.779 500 Speed Up = 192.7699x

hello_world.cpp

#include <Rcpp.h>
using namespace Rcpp;

// Below are simple examples of exporting a C++ function to R. 
// Source this function into an R session using the Rcpp::sourceCpp 

// Every function needs to start with "// [[Rcpp::export]]"

// [[Rcpp::export]]
int timesTwo(int x) {
   return x * 2;
}

// [[Rcpp::export]]
double sqrt(double x) {
   return sqrt(x);
}


// [[Rcpp::export]]
double meanC(NumericVector x) {
  int n = x.size();
  double total = 0;
  for(int i = 0; i < n; ++i) {
    total += x[i] / n;
    }
  return total;
}

simple_comparison.R

## =============================================================================
## Simple Examples: Pure R vs C++ Speed Comparison
## =============================================================================

rm(list=ls())
setwd("/media/SUPPORT/Repo/sandbox/cpp")
library(Rcpp)
library(microbenchmark)
set.seed(1234)

## =============================================================================
## Parameters
## =============================================================================

n_fib <- 15
n_times <- 500


## =============================================================================
## Fibonacci recursive
## =============================================================================

## Pure R ----------------------
fibR <- function(n) {
  if (n == 0) return(0)
  if (n == 1) return(1)
  return (fibR(n - 1) + fibR(n - 2))
}

## CPP --------------------------
cppFunction('
  int fibC(const int x) {
    if (x == 0) return (0);
    if (x == 1) return (1);
    return fibC(x - 1) + fibC(x - 2);
  }')

## Benchmark ------------------------
res <- summary(microbenchmark(fibR(n_fib), fibC(n_fib), times = n_times))
cat("\nExample: Fibonacci Recursive\n")
print(res)
cat("Speed Up = ", res[1,]$median / res[2,]$median, "x\n", sep = "")


## =============================================================================
## Fibonacci for loop
## =============================================================================

## linear / iterative solution ---------------
fibRiter <- function(n) {
  first <- 0
  second <- 1
  third <- 0
  for (i in seq_len(n)) {
    third <- first + second
    first <- second
    second <- third
  }
  return(first)
}

## CPP -----------------
cppFunction('
  int fibCiter(const int n) {
    double first = 0;
    double second = 1;
    double third = 0;
    for (int i = 0; i < n; i++) {
      third = first + second;
      first = second;
      second = first;
    }
    return(first);
  }')

## Benchmark ---------------------
res <- summary(microbenchmark(fibRiter(n_fib), fibCiter(n_fib), times = n_times))
cat("\nExample: Fibonacci For Loop\n")
print(res)
cat("Speed Up = ", res[1,]$median / res[2,]$median, "x\n", sep = "")


## =============================================================================
## VAR Simulation
## =============================================================================

## parameter and error terms used throughout ---------------------------
a <- matrix(c(0.5,0.1,0.1,0.5),nrow=2)
u <- matrix(rnorm(10000),ncol=2)

## Let’s start with the R version -------------------
rSim <- function(coeff, errors) {
  simdata <- matrix(0, nrow(errors), ncol(errors))
  for (row in 2:nrow(errors)) {
    simdata[row,] = coeff %*% simdata[(row-1),] + errors[row,]
  }
  return(simdata)
}
rData <- rSim(a, u) # generated by R

## Cpp -------------------------------------------------------------
## Now load ’inline’ to compile C++ code on the fly
suppressMessages(require(inline))
code <- '
  arma::mat coeff = Rcpp::as<arma::mat>(a);
  arma::mat errors = Rcpp::as<arma::mat>(u);
  int m = errors.n_rows;
  int n = errors.n_cols;
  arma::mat simdata(m,n);
  simdata.row(0) = arma::zeros<arma::mat>(1,n);
  for (int row=1; row<m; row++) {
    simdata.row(row) = simdata.row(row-1) * trans(coeff) + errors.row(row);
  }
  return Rcpp::wrap(simdata);'

## create the compiled function
rcppSim <- cxxfunction(signature(a = "numeric", u = "numeric"),
                       code, plugin="RcppArmadillo")
rcppData <- rcppSim(a,u) # generated by C++ code
stopifnot(all.equal(rData, rcppData)) # checking results

## Compare ------------------------------
res <- summary(microbenchmark(rSim(a,u), rcppSim(a,u), times = n_times))
cat("\nExample: VAR Simulation\n")
print(res)
cat("Speed Up = ", res[1,]$median / res[2,]$median, "x\n", sep = "")