### Present Remotely

Send the link below via email or IM

• Invited audience members will follow you as you navigate and present
• People invited to a presentation do not need a Prezi account
• This link expires 10 minutes after you close the presentation

Do you really want to delete this prezi?

Neither you, nor the coeditors you shared it with will be able to recover it again.

# R and C++ Integration

No description
by

## Suz Tolwinski-Ward

on 1 September 2015

Report abuse

#### Transcript of R and C++ Integration

by Suz Tolwinski-Ward
Souped-up R code with C++
Problem:
Despite being awesome in so many ways (open-source, extensible, user community, interactivity), sometimes R is too slow (it's interpreted, not compiled).
Solution:
Rewrite bottlenecks in C++ callable from R with the Rcpp package.
real modeling problem: tropical cyclone landfalls
Determine algorithmically whether the simulated TC tracks (red, green) intersect the coastline (blue) from their piecewise linear representations.
Solution: check whether each pair of segments from the coastline and the TC track intersect.
(Beware the NESTED LOOPS!!!!)
Want to be able to do this zillions of times for Monte Carlo sampling of zillions of tracks.
Moral of the demo: my C++ implementation, called from R via ,
is several hundreds of times faster than my plain-old R implementation!
Rcpp
orient3pts <- function(pt1,pt2,pt3){

orient <- (pt2[2]-pt1[2])*(pt3[1]-pt2[1]) -
(pt2[1]-pt1[1])*(pt3[2]-pt2[2]);

# orient == 0 ==> colinear triple of points
# orient > 0 ==> clockwise orientation
# orient < 0 ==> counter-clockwise orientation

return(sign(orient))

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

// [[Rcpp::export]]
int sgn_cpp(double x) {
int y = 0;
if (0 < x){y = 1;}
if (0 > x){y = -1;}
return y;
}

// [[Rcpp::export]]
int orient3pts_cpp(NumericVector pt1,NumericVector pt2,
NumericVector pt3){

double orient = (pt2[1]-pt1[1])*(pt3[0]-pt2[0]) -
(pt2[0]-pt1[0])*(pt3[1]-pt2[1]);

return(sgn_cpp(orient));

}
use the Rcpp namespace...
...and export these functions for
use within R sessions.
(Note: if you create a new C++ file from the Rstudio drop-down menu, it opens a template file for you with the "magical incantations" already filled in.)
#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
List rowrep_vec (int n, NumericVector x){

int m = x.size();

NumericMatrix y(n,m);

for(int i = 0; i < n; i++){
for(int j = 0; j < m; j++){
y(i,j) = x(j);
}
}

return List::create(
_["input_vec"] = x,
_["y"] = y);

}

}
For more "syntactic sugar," see: http://dirk.eddelbuettel.com/code/rcpp/Rcpp-sugar.pdf
Useful "Quick Reference" for your cube wall: https://cran.r-project.org/web/packages/Rcpp/vignettes/Rcpp-quickref.pdf
#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
NumericVector rnorm_trunc(int n, double mu,
double sig,NumericVector lb,NumericVector ub){

NumericVector u_lb = pnorm(lb,mu,sig);
NumericVector u_ub = pnorm(ub,mu,sig);
NumericVector u = runif(n,u_lb(0),u_ub(0));

return(qnorm(u,mu,sig));

}

using namespace Rcpp;
using namespace arma;
using namespace std;

// Compute the likelihood of an AR(1) model with 2 linear predictors
// [[Rcpp::export()]]
arma::vec lnL_ar1model_2pred(NumericVector y, NumericVector x1,
NumericVector x2,NumericVector beta,
double phi, double sig2){

int N = y.size();

arma::mat ImBphi(N,N);
for(int n = 1; n < N; n++){
ImBphi(n,n) = 1;
ImBphi(n,n-1) = -phi;
}
ImBphi(0,0) = 1;

arma::mat X(N,3);
for(int i = 0; i < N; i++){
X(i,0) = 1;
X(i,1) = x1(i);
X(i,2) = x2(i);
}

arma::mat Sigma = sig2 * arma::inv(arma::trans(ImBphi) * ImBphi);

arma::vec lnL = dmvn_arma(arma::trans(as<arma::vec>(y)),
arma::trans(X * as<arma::vec>(beta)),Sigma,true);
return(lnL);

}

// from http://gallery.rcpp.org/articles/dmvnorm_arma/
// [[Rcpp::export]]
arma::vec dmvn_arma(arma::mat x, arma::rowvec mean,
arma::mat sigma,
bool logd = false) {
int n = x.n_rows;
int xdim = x.n_cols;
arma::vec out(n);
arma::mat rooti = arma::trans(arma::inv(trimatu(arma::chol(sigma))));
double rootisum = arma::sum(log(rooti.diag()));
double constants = -(static_cast<double>(xdim)/2.0) * log2pi;

for (int i=0; i < n; i++) {
arma::vec z = rooti * arma::trans( x.row(i) - mean) ;
out(i) = constants - 0.5 * arma::sum(z%z) + rootisum;
}

if (logd == false) {
out = exp(out);
}
return(out);
}
Solution "Recipe":
Step 0. install Rcpp and import the library
> install.packages('Rcpp')
> library(Rcpp)
Step 1. Profile your R code to find the bottlenecks
Step 2. Port the offending chunks of code to C++ and save file (as, say )
myfunc.cpp
Step 3. Source the C++ code in your R session
> sourceCpp('mypath/myfunc.cpp')
Step 4. Call your C++ function in your R session and burn rubber!
> Rprof()
> for(i in 1:100){foo.test <- myfunc_r(args)}
> Rprof(NULL)
> summaryRprof()
(Examples of this forthcoming in the rest of the prezi!)
deliciously
fast
R code
An offending chunk of R code in Ex 1:
C++ port of offending code chunk:
Example 1
Example 2
Example 3
Example 4
Rcpp delivers awesome
speedups!
Rcpp has convenient matrix,
vector, list data types for I/O...
...and d/p/q/r random number gen. a la R
Matrix algebra is possible
#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]

<Run demo in Rcpp_tutorial.r here!>
variable declarations,
lines ending in semi-colons
built-in vector type!
indexing starts at 0
Much more complete intro slides to Rcpp and RcppArmadillo: http://learn.stat.ubc.ca/~andy.leung/files/seminars-talks/2012/03/RcppDemo.pdf