Loading presentation...

Present Remotely

Send the link below via email or IM

Copy

Present to your audience

Start remote presentation

  • 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
  • A maximum of 30 users can follow your presentation
  • Learn more about this feature in our knowledge base article

Do you really want to delete this prezi?

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

DeleteCancel

Make your likes visible on Facebook?

Connect your Facebook account to Prezi and let your likes appear on your timeline.
You can change this under Settings & Account at any time.

No, thanks

R and C++ Integration

No description
by

Susan Tolwinski-Ward

on 1 September 2015

Comments (0)

Please log in to add your comment.

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));

}
include the Rcpp header...
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
Getting started with porting your R to C++: http://adv-r.had.co.nz/Rcpp.html
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));

}

//[[Rcpp::depends(RcppArmadillo)]]

#include <RcppArmadillo.h>
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
Bookmark these links!
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
with RcppArmadillo
#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]

your C++ code here
<Run demo in Rcpp_tutorial.r here!>
variable declarations,
lines ending in semi-colons
built-in vector type!
indexing starts at 0
Great slides on RcppArmadillo: http://q-aps.princeton.edu/sites/default/files/q-aps/files/slides_day4_am.pdf
Much more complete intro slides to Rcpp and RcppArmadillo: http://learn.stat.ubc.ca/~andy.leung/files/seminars-talks/2012/03/RcppDemo.pdf
Armadillo C++ library documentation: http://arma.sourceforge.net/docs.html
a few more incantations
matrix inverse
matrix transpose
matrix multiplication
convert to an arma type
Cholesky
decomposition
interpret matrix as
triangular
element-wise
multiplication
Full transcript