Interfacing C++ code from R Sren Hjsgaard Department of - PowerPoint PPT Presentation

Interfacing C++ code from R Søren Højsgaard Department of Mathematical Sciences Aalborg University, Denmark November 20, 2012 Printed: November 20, 2012 File: interfaceCpp-slides.tex

2 Contents 1 Calling C++ from R 3 2 Some important libraries and packages 4 3 Example: The exponential function 6 3.1 Using Rcpp together with inline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 4 Example: Calculating Fibonacci recursively 12 5 Example: Matrix multiplication 15 5.1 Using Rcpp together with inline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 5.2 Using RcppArmadillo together with inline . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 5.3 Benchmarking – III . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 6 Compiling without using inline 21 6.1 Example: Matrix multiplication using Rcpp . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 6.2 Example: Matrix multiplication using RcppArmadillo . . . . . . . . . . . . . . . . . . . . . . 27 6.3 Example: Inverting a symmetric positive definite matrix using RcppArmadillo . . . . . . . . . . 28 7 Further examples on using RcppArmadillo with inline 30 7.1 Example: Extracting submatrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 8 Calling R from C++ 33 9 Building packages using Rcpp libraries 36 10 EXERCISES 37

3 1 Calling C++ from R C++ code can be called from R . • Easily done using the Rcpp package. (There are other ways not to be discussed here). • The Rcpp package combines nicely with the inline package. • There are several existing libraries available with a C++ interface that we may use – instead of“reinventing the wheel”(as we did when creating our own matrix multiplication functions).

4 2 Some important libraries and packages • Armadillo is an open-source C++ linear algebra library (matrix maths) aiming towards a good balance between speed and ease of use. http://arma.sourceforge.net/ R interface via RcppArmadillo • Eigen is a C++ template library for linear algebra: matrices, vectors, numerical solvers, and related algorithms. http://eigen.tuxfamily.org/ R interface via RcppEigen • The GNU Scientific Library (GSL) is a numerical library for C and C++ programmers. http://www.gnu.org/software/gsl/ R interface via RcppGSL

5 Notice: • To use these packages we must program in C++ rather than in C . • Dirk Eddelbuettel provides many many examples on his website: http://dirk.eddelbuettel.com/code/rcpp.html

6 3 Example: The exponential function Recall our implementation of the exponential function in R : > expfunR <- function(x, eps=.Machine$double.eps){ + if (x<0) { + 1/expfunR(-x) + } else { + n <- 0; term <- 1; ans <- 1 + while(abs(term)> eps) { + n = n + 1 + term = term * x / n + ans <- ans + term + } + ans + } + }

7 A pure C implementation (which ignores the numerical difficulties when x < 0 ) is: 1 #include <math.h> 2 double C_expfun2 (double x){ 3 double ans =1.0 , term =1.0 , eps =1e -16; 4 int n=0; 5 while (fabs(term)>eps ){ 6 n++; 7 term = term * x / n; 8 ans = ans + term; 9 } 10 return(ans ); 11 }

8 3.1 Using Rcpp together with inline > src <- ' + double x = as<double>(x_); + double ans=1.0, term=1.0, eps=1e-16; + int n=0; + while (fabs(term)>eps){ + n++; + term = term * x / n; + ans = ans + term; + } return(wrap(ans)); ' + Compiling, linking and loading is done with one function call > library(inline) > expfunC <- cxxfunction(signature(x_="numeric"), plugin="Rcpp", body=src) • R data object ( SEXP s) converted to C++ objects with as() • C++ data object converted R data objects ( SEXP s) with wrap()

9 An alternative implementation is to use the original function (accounting for nummerical difficulties when x < 0 has been resolved): > incltxt <- ' double C_expfun2 (double x){ + double ans=1.0, term=1.0, eps=1e-16; + int n=0; + while (fabs(term)>eps){ + n++; + term = term * x / n; + ans = ans + term; + } + return(ans); + } ' > expfunC2 <- cxxfunction(signature(x_="numeric"), includes=incltxt, plugin="Rcpp", body= ' + + double x = as<double>(x_); + if (x>0){ + return wrap(C_expfun2(x)); + } else { + return wrap(1/C_expfun2(-x)); + } + ' )

10 > xx <- 1 > c(expfunR(xx), expfunC(xx), expfunC2(xx)) [1] 2.718282 2.718282 2.718282 > xx <- 20 > c(expfunR(xx), expfunC(xx), expfunC2(xx)) [1] 485165195 485165195 485165195 > xx <- -20 > c(expfunR(xx), expfunC(xx), expfunC2(xx)) [1] 2.061154e-09 5.621884e-09 2.061154e-09 > (expfunC(-20)-exp(-20))/exp(-20) [1] 1.727543 > (expfunC2(-20)-exp(-20))/exp(-20) [1] 2.006596e-16 >

11 > library(rbenchmark) > cols <- c("test", "replications", "elapsed", "relative") > N <- 20000 > benchmark(expfunR(10), expfunC(10), expfunC2(10), + columns=cols, order="relative", replications=N) test replications elapsed relative 2 expfunC(10) 20000 0.08 1.0 3 expfunC2(10) 20000 0.08 1.0 1 expfunR(10) 20000 2.44 30.5

12 4 Example: Calculating Fibonacci recursively The standard definition of the Fibonacci sequence is f n = f n − 1 + f n − 2 where f 0 = 0 , f 1 = 1 . A simple recursive implementation in R is: > fibR <- function(n){ + if (n==0) + return(0) + else + if (n==1) + return(1) + else + return(fibR(n-1)+fibR(n-2)) + }

13 Easy to write fast C++ version: > incltxt <- ' + int fibonacci(int n){ + if (n==0) + return(0); + else + if (n==1) + return(1); + else + return(fibonacci(n-1)+fibonacci(n-2)); + } ' > fibC <- cxxfunction(signature(x_="int"), + plugin="Rcpp", + include=incltxt, + body= ' + int x = as<int>(x_); + return wrap( fibonacci(x) ); + ' )

14 > library(rbenchmark) > cols <- c("test", "replications", "elapsed", "relative") > M <- 30 > benchmark(fibR(M), fibC(M), + columns=cols, order="relative", replications=1) test replications elapsed relative 1 fibR(M) 1 7.86 NA 2 fibC(M) 1 0.00 NA

15 5 Example: Matrix multiplication Recall our first version of the matrix multiplication function: 1 /* File: matprod1.c: Calculate the product of matrices X and Y */ 2 void matprod1(double *X, int *dimX , double *Y, int *dimY , double *ans ){ 3 double sum; 4 int ii , jj , kk; 5 int nrX=dimX [0], ncX=dimX [1], nrY=dimY [0], ncY=dimY [1]; 6 7 for (ii =0; ii <nrX; ii ++){ 8 for (jj =0; jj <ncY; jj ++){ 9 sum = 0; 10 for (kk =0; kk <ncX; kk ++){ 11 sum = sum + X[ii+nrX*kk]*Y[kk+nrY*jj]; 12 } 13 ans[ii+nrX*jj] = sum; 14 } 15 } 16 }

16 An interface using SEXP s is: 1 /* File: matprod2.c: Calculates the product of matrices X and Y */ 2 #include <Rdefines.h> 3 #include "matprod1.h" 4 5 SEXP matprod2(SEXP X, SEXP Y) { 6 int nprot =0; 7 PROTECT(X = AS_NUMERIC (X)); nprot ++; /* Digest SEXPs from R */ 8 PROTECT(Y = AS_NUMERIC (Y)); nprot ++; 9 double *xptr; xptr = REAL(X); 10 double *yptr; yptr = REAL(Y); 11 int *dimX; dimX = INTEGER(GET_DIM(X)); 12 int *dimY; dimY = INTEGER(GET_DIM(Y)); 13 SEXP ans; /* Create SEXP to hold result */ 14 PROTECT(ans = allocMatrix (REALSXP , dimX [0], dimY [1])); nprot ++; 15 double *ansptr; ansptr = REAL(ans ); 16 matprod1(xptr , dimX , yptr , dimY , ansptr ); /* Calculate product */ 17 UNPROTECT(nprot ); /* Wrap up; */ 18 return(ans ); /* Return the result to R */ 19 } where 1 void matprod1(double *X, int *dimX , double *Y, int *dimY , double *ans );

17 5.1 Using Rcpp together with inline With Rcpp matrices can be indexed the usual way: > src <- ' + NumericMatrix X(X_); + NumericMatrix Y(Y_); + NumericMatrix ans (X.nrow(), Y.ncol()); + int ii, jj, kk; + for (ii=0; ii<X.nrow(); ii++){ + for (jj=0; jj<Y.ncol(); jj++){ + double sum = 0; + for (kk=0; kk<X.ncol(); kk++){ + sum += X(ii, kk) * Y(kk, jj); + } + ans(ii,jj) = sum; + } + } + return(ans); ' > library(inline) > mprod5_inline_Rcpp <- cxxfunction(signature(X_="numeric", Y_="numeric"), + body = src, plugin="Rcpp")

18 5.2 Using RcppArmadillo together with inline The Armadillo library is an excellent C++ package for linear algebra and RcppArmadillo makes this easy > src <- ' + arma::mat X = Rcpp::as<arma::mat>(X_); + arma::mat Y = Rcpp::as<arma::mat>(Y_); + arma::mat ans = X * Y; + return(wrap(ans)); + ' > mprod6_inline_RcppArma <- cxxfunction(signature(X_="numeric", Y_="numeric"), + body = src, plugin="RcppArmadillo")

19 5.3 Benchmarking – III > library(rbenchmark) > cols <- c("test", "replications", "elapsed", "relative") > N <- 10 > A <- matrix(rnorm(1000000), nr=1000); B<-A[,1:5] > benchmark(A %*% B, mprod5_inline_Rcpp(A,B), mprod6_inline_RcppArma(A,B), + columns=cols, order="relative", replications=N) test replications elapsed relative 1 A %*% B 10 0.05 1.0 3 mprod6_inline_RcppArma(A, B) 10 0.06 1.2 2 mprod5_inline_Rcpp(A, B) 10 1.26 25.2