Rcpp Gallery

Extending R with C++ and Fortran

Sun, 05 Feb 2023 00:00:00 +0000

A recent social-media question by James Curran inquired about the best, or recommended ways, to extend R with Fortran code. Part of the question was whether the .Fortran() interface was still recommended or not as there is ‘conflicting advice’ out there. Dirk then followed up and pointed to the (stunning!) performance gains reported by glmnet which switched from .Fortran() to a C++ interface using Rcpp and the (now much preferred) .Call() interface. One key reason behind the performance gains is that .Fortran() requires copies of all arguments, just like the (also effectively deprecated) .C() interface. Whereas .Call() works with SEXP objects which are pointers: this can be dramatically faster and more efficient as object sizes increase.

A few years earlier, and for a related question, JBrandon Duck-Mayr had written a very comprehensive answer on StackOverflow. It is backed by an example package mixedlang which implements the recommendation.

It starts from a Fortran90 function multiplying two ‘real’ aka double valued inputs:

REAL*8 FUNCTION MULTIPLY (X, Y) 
REAL*8 X, Y
MULTIPLY = X * Y
RETURN
END

This can be connected quite easily to C++ code using the common extern "C"declaration (specifying that a C calling convention is used from the C++ code). It still shows the Rcpp::depends() used when sourceCpp()-ing a function, it is not needed in a package like mixedlang.

#include "RcppArmadillo.h"

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

// First we'll declare the MULTIPLY Fortran function
// as multiply_ in an extern "C" linkage specification
// making sure to have the arguments passed as pointers.
extern "C" {
    double multiply_(double *x, double *y);
}

// Now our C++ function
// [[Rcpp::export]]
Rcpp::NumericVector test_function(Rcpp::NumericVector x) {
    // Get the size of the vector
    int n = x.size();
    // Create a new vector for our result
    Rcpp::NumericVector result(n);
    for ( int i = 0; i < n; ++i ) {
        // And for each element of the vector,
        // store as doubles the element and the index
        double starting_value = x[i], multiplier = (double)i;
        // Now we can call the Fortran function,
        // being sure to pass the address of the variables
        result[i] = multiply_(&starting_value, &multiplier);
    }
    return result;
}

Once both functions are compiled and loaded (as e.g. in package mixedlang) the wrapper function can be called from R as usual:

mixedlang::test_function(0:9)
# [1]  0  1  4  9 16 25 36 49 64 81

We hope the (recently updated) package at GitHub serves as starting point for others wanting to combine R and Fortran via Rcpp.

Benchmarking Rcpp code with RcppClock

Sat, 06 Nov 2021 00:00:00 +0000

Rcpp is all about the need for speed. However, Rcpp code is not magically fast, but is the result of careful coding and profiling. RcppClock makes it easy to profile C++ code from R:

On the C++ side, RcppClock wraps std::chrono::high_resolution_clock (C++11) to time chunks of code, and exports the results to an R object in the global environment.
On the R side, RcppClock provides convenient methods for summarizing results and plotting with ggplot2.

There are other Rcpp benchmarking tools, such as the Rcpp timer class (see RcppGallery vignettte here) and rcppgeiger (see repo here). RcppClock moves beyond both of these tools by supporting overlapping time points and direct access to timing results in the R environment.

By benchmarking Rcpp code with RcppClock, we avoid measuring latency associated with converting R types to C++ types (and back again), and we can profile at line-by-line resolution rather than entire functions.

Note: RcppClock is not compatible with OpenMP parallelization, and does not accurately measure processes taking less than ~1 microsecond (somewhat implementation-dependent).

A simple example

Install RcppClock from CRAN:

install.packages("RcppClock")

In your .cpp file, link the RcppClock header with //[[Rcpp::depends(RcppClock)]] and, if this is an R package, link to it in your DESCRIPTION file. RcppClock loads Rcpp.h so we don’t need to include it again.

//[[Rcpp::depends(RcppClock)]]
#include <RcppClock.h>
#include <thread>

//[[Rcpp::export]]
void sleepy(){
  Rcpp::Clock clock;
  
  clock.tick("both_naps");
  
  clock.tick("short_nap");
  std::this_thread::sleep_for(std::chrono::milliseconds(10));  
  clock.tock("short_nap");
  
  clock.tick("long_nap");
  std::this_thread::sleep_for(std::chrono::milliseconds(100));  
  clock.tock("long_nap");

  clock.tock("both_naps");
  
  // send the times to the R global environment variable, named "naptimes"
  clock.stop("naptimes");
}

.tick(std::string) starts a new timer. Provide a name to record what is being timed.

.tock(std::string) stops a timer. It is important to use the same name as declared in .tick().

.stop(std::string) calculates the duration between all .tick() and .tock() timing results, and creates an object in the R environment with the name provided.

On the R end, we can now summarize the results using the “naptimes” variable that was created in the above Rcpp function:

sleepy()
# global variable "naptimes" is now created in the environment
library(RcppClock)
naptimes

Unit: milliseconds 
    ticker   mean sd    min    max neval
 both_naps 110.19 NA 110.19 110.19     1
  long_nap 100.10 NA 100.10 100.10     1
 short_nap  10.09 NA  10.09  10.09     1

Multiple replicates

For extremely fast calculations it can be useful to repeat the calculation some number of times in order to improve the accuracy of the measurement.

If a .tick() with the same name is called multiple times, RcppClock automatically groups the results.

The following code reproduces the ?fibonacci function example included in the RcppClock package:

int fib(int n) {
  return ((n <= 1) ? n : fib(n - 1) + fib(n - 2));
}

//[[Rcpp::export]]
void fibonacci(std::vector<int> n, int reps) {
  Rcpp::Clock clock;
  
  while(reps-- > 0){
    for(auto number : n){
      clock.tick("fib" + std::to_string(number));
      fib(number);
      clock.tock("fib" + std::to_string(number));
    }
  }
  clock.stop("clock");
}

On the R end, we will get an object named “clock”:

fibonacci(n = 25:35, reps = 10)
# global variable "clock" is created in the R global environment
clock

Unit: milliseconds 
 ticker    mean      sd     min     max neval
  fib25  0.1734 0.04975  0.1354  0.2845    10
  fib26  0.2651 0.06472  0.2263  0.4273    10
  fib27  0.3932 0.07265  0.3493  0.5869    10
  fib28  0.6463 0.12682  0.5626  0.9919    10
  fib29  1.0641 0.22482  0.8996  1.5922    10
  fib30  1.7432 0.53559  1.3993  2.7752    10
  fib31  2.7441 0.67205  2.4041  4.6192    10
  fib32  4.0757 0.73495  3.6354  6.1233    10
  fib33  6.8620 0.49135  6.2664  7.7575    10
  fib34 10.0244 0.22531  9.7690 10.4643    10
  fib35 17.9276 1.05583 17.1352 20.5875    10

The result of summary() is an S3 object of class RcppClock (essentially a data.frame) and has a plot method using ggplot2:

plot(clock)

Unfortunately, some times were nearly impossible to measure. We can fix this by replicating the calculation many times within each timing bracket:

//[[Rcpp::export]]
void fibonacci_rep(std::vector<int> n, int outer_reps, int inner_reps) {
  Rcpp::Clock clock;
  
  while(outer_reps-- > 0){
    for(auto number : n){
      clock.tick("fib" + std::to_string(number));
      for(int rep = 0; rep < inner_reps; ++rep)
        fib(number);
      clock.tock("fib" + std::to_string(number));
    }
  }
  clock.stop("clock");
}

Now on the R end we can account for this:

fibonacci_rep(n = 10:20, outer_reps = 10, inner_reps = 10000)
clock$timer <- clock$timer / 10000
clock

Unit: nanoseconds 
 ticker     mean      sd      min      max neval
  fib10    77.66   4.469    71.67    87.31    10
  fib11   116.92   8.336   109.19   133.93    10
  fib12   206.95   7.735   194.62   218.52    10
  fib13   306.46  18.058   290.51   349.29    10
  fib14   637.75  69.371   582.90   815.85    10
  fib15   909.88  45.979   864.89  1020.59    10
  fib16  1782.64  45.661  1735.67  1866.92    10
  fib17  2642.99  84.156  2562.64  2832.20    10
  fib18  5170.67 126.798  5063.44  5456.16    10
  fib19  7485.84 107.298  7319.87  7685.55    10
  fib20 14326.20 214.079 14063.28 14684.91    10

And plot the results:

plot(clock)

In conclusion, RcppClock is a non-intrusive method for benchmarking Rcpp code in R.

Simulation Smoother using RcppArmadillo

Sun, 05 Sep 2021 00:00:00 +0000

Simulation smoothing involves drawing latent state variables in discrete-time state–space models from their conditional distribution given parameters and data as defined by McCausland, Miller, & Pelletier (2011). This tool greatly facilitates Bayesian estimation in the case of a popular family of conditionally linear Gaussian state-space models. This post presents Rcpp implementations of a sampling algorithm for a multivariate normal distribution. RcppArmadillo is used to efficiently calculate the Cholesky decomposition and implement the back-substitution to solve linear equations, which can greatly improves computation times when using routines for banded or tridiagonal matrices.

Precision sampler

Consider the problem of sampling random numbers from a multivariate normal distribution specified by Chan & Jeliazkov (2009) as: \begin{equation}\label{eq1} \mathcal{N}_T\left(\Omega^{-1}\alpha, \Omega^{-1}\right) \end{equation} where:

the covariance matrix, denoted by $\Omega^{-1}$, is a $T\times T$ symmetric positive-definite matrix given by the inverse of a precision matrix $\Omega$ that, depending on the autoregressive lag order of the state equation, is band diagonal or tridiagonal and easy to compute, and
the mean of the distribution, denoted by $\Omega^{-1}\alpha$, is constructed as the product of the inverse of the precision matrix times an easy-to-compute $T$-vector $\alpha$.

Let $L$ be a lower-triangular matrix obtained by the Cholesky decomposition of the precision matrix $\Omega=LL’$. Consequently, the mean of the normal distribution above can be presented as $(LL’)^{-1}\alpha = L^{-1\prime}L^{-1}\alpha$. Let a $T$-column vector $\epsilon$ contain independent random draws from the standard normal distribution. A random draw from the target distribution can be computed then by: \begin{equation}\label{eq2} L^{-1\prime}L^{-1}\alpha + L^{-1\prime}\epsilon = L^{-1\prime}\left(L^{-1}\alpha + \epsilon\right) \end{equation} The formulation of the problem on the right-hand side of the equation above follows the proposition of McCausland, Miller, & Pelletier (2011) and leads to the sampling algorithm following the steps:

Compute the Cholesky decomposition of matrix, $\Omega = LL’$. Use dedicated computational routines for band or tridiagonal matrices if possible.
Sample $T$ independent draws from the standard normal distribution and collect them in the column vector $\epsilon$.
Solve the linear equation $La = \alpha$ for $a$ using forward-substitution.
Solve the linear equation $L’X = a + \epsilon$ for $X$ using back-substitution.
Return $X$ as a draw from the target distribution.

The computational gains from applying the algorithm presented above are twofold. First, taking into account that matrix $\Omega$ is band or tridiagonal makes computing the Cholesky decomposition fast if dedicated routines are applied and leads to a lower-triangular matrix $L$ with as many non-zero subdiagonals as in the precision matrix. Secondly, solving the linear equations using back-substitution bypasses the computation of the inverse of $L$ and takes advantage of the special type of matrix $L$. This strategy greatly reduces the computational time required to obtain the random draws and, therefore, was recommended by Frühwirth-Schnatter & Kastner (2014) and Bitto & Frühwirth-Schnatter (2019) following the developments by Chan & Jeliazkov (2009) and McCausland, Miller, & Pelletier (2011).

Three implementations in RcppArmadillo

All the presented solutions are designed to sample n draws from the target distribution and return them in columns of the output matrix draw. The results are required to be reproducible. Therefore, the compatibility with R’s set.seed is obtained by generating normal random draws using Rcpp::rnorm rather than arma::randn. All the considered functions take three arguments: the number of requested draws n, the precision matrix, and the covector location. Then they read the dimension of the normal vector T, fill matrix epsilon with normal random draws, create the location_matrix of appropriate dimensions whose columns are filled with location using each_col, and proceed to the implementation of steps 3 and 4 from the sampling algorithm presented above.

The first function rmvnorm_arma_inv is based on the Cholesky decomposition computed by chol subsequently declared as an upper-triangular matrix using trimatu, and uses function inv to compute the inverse. This matrix is then used to compute the draws from the target distribution using matrix multiplication and summation.

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]

using namespace Rcpp;
using namespace arma;

// [[Rcpp::export]]
mat rmvnorm_arma_inv (int n, mat precision, vec location){
  
  int T = precision.n_rows;
  
  mat epsilon(T, n);
  for (int i=0; i<n; i++){
    epsilon.col(i) = as<vec>(rnorm(T));
  }
  mat location_matrix(T, n, fill::zeros);
  location_matrix.each_col() += location;
  
  mat precision_chol_inv = trans(inv(trimatu(chol(precision))));
  mat draw    = trans(precision_chol_inv) * (
    precision_chol_inv * location_matrix + epsilon
    );

  return draw;
}

The second function rmvnorm_arma_solve uses chol and trimatu the same way as the first one. However, this implementation uses solve searching for the solutions to the linear equations using back-substitution. Additionally, function solve checks if the matrix is band and applies dedicated numerical procedures as explained in the API Documentation for Armadillo 10.6. This additional check facilitates the computational gains. The development of the two functions considered thus far was somehow inspired by the RcppGallery article by Hardt, Ahmadou, & Christoffersen (2013-2020).

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]

using namespace Rcpp;
using namespace arma;

// [[Rcpp::export]]
mat rmvnorm_arma_solve(int n, mat precision, vec location){
  
  int T = precision.n_rows;
  
  mat epsilon(T, n);
  for (int i=0; i<n; i++){
    epsilon.col(i) = as<vec>(rnorm(T));
  }
  mat location_matrix(T, n);
  location_matrix.each_col() = location;
  
  mat precision_chol  = trimatu(chol(precision));
  mat draw            = solve(precision_chol, 
                              solve(trans(precision_chol), 
                                    location_matrix) + epsilon);
  
  return draw;
}

The third function rmvnorm_arma_stochvol is based on three utility functions from package stochvol by Kastner (2016) and Hosszejni & Kastner. These functions implement dedicated numerical algorithms described in McCausland, Miller, & Pelletier (2011) taking advantage of the fact that, in many applications of the state-space models, all the elements on the first sub-diagonal of the precision matrix are given by the same number. In this case, the algorithms for the Cholesky decomposition, forward- and back-substitution simplify granting further computational time gains. The three functions taken from the stochvol package namespace are cholesky_tridiagonal, forward_algorithm, and backward_algorithm. The only changes implemented in the functions presented in the listing below consider the definitions of object types that were simplified. The function rmvnorm_arma_stochvol exhibits a slightly adjusted structure relative to the first two presented functions for the simulation smoother accommodating the fact that the functions forward_algorithm, and backward_algorithm accept only a vector as the third argument and not a matrix.

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]

using namespace Rcpp;
using namespace arma;

/*
 The three functions below, cholesky_tridiagonal, forward_algorithm, 
 and backward_algorithm were copied from the source code of package 
 stochvol by Georg Kastner that was shared under licence GPL (>= 2). 
 They were modified by the author of the current note on 30 August 2021
 and are shared under the GPL (>= 2) licence as well.
 */
List cholesky_tridiagonal(const vec& omega_diag, const double& omega_offdiag) {
  const int T = omega_diag.n_elem - 1;
  vec chol_diag(T+1);
  vec chol_offdiag(T+1);
  chol_diag[0] = std::sqrt(omega_diag[0]);
  for (int j = 1; j < T+1; j++) {
    chol_offdiag[j-1] = omega_offdiag/chol_diag[j-1];
    chol_diag[j] = std::sqrt(omega_diag[j]-chol_offdiag[j-1]*chol_offdiag[j-1]);
  }
  return List::create(_["chol_diag"]=chol_diag, _["chol_offdiag"]=chol_offdiag);
}

vec forward_algorithm(const vec& chol_diag, const vec& chol_offdiag, const vec& covector) {
  const int T = chol_diag.n_elem - 1;
  vec htmp(T+1);
  htmp[0] = covector[0]/chol_diag[0];
  for (int j = 1; j < T+1; j++) {
    htmp[j] = (covector[j] - chol_offdiag[j-1]*htmp[j-1])/chol_diag[j];
  }
  return htmp;
}

vec backward_algorithm(const vec& chol_diag, const vec& chol_offdiag, const vec& htmp) {
  const int T = chol_diag.size() - 1;
  vec h(T+1);
  h[T] = htmp[T] / chol_diag[T];
  for (int j = T-1; j >= 0; j--) {
    h[j] = (htmp[j] - chol_offdiag[j] * h[j+1]) / chol_diag[j];
  }
  return h;
}

// [[Rcpp::export]]
mat rmvnorm_arma_stochvol(int n, mat precision, vec location){
  
  int     T                 = precision.n_rows;
  vec     precision_diag    = precision.diag();
  double  precision_offdiag = precision(1,0);
  
  List  precision_chol  = cholesky_tridiagonal(precision_diag, precision_offdiag);
  vec   aa              = forward_algorithm(precision_chol["chol_diag"],
                                            precision_chol["chol_offdiag"],
                                                          location);
  mat draw(T, n);
  vec epsilon;
  for (int i=0; i<n; i++){
    epsilon     = rnorm(T);
    draw.col(i) = backward_algorithm(precision_chol["chol_diag"], 
                                     precision_chol["chol_offdiag"],
                                                   aa + epsilon);
  }

  return draw;
}

Three implementations in R

In addition to the cpp functions, three R functions are presented for the sake of comparing the computational time of their execution. These functions exhibit a similar structure as the cpp functions with respect to the definition of the arguments, outputs, and the order of implementing the sampling algorithm steps.

The first function rmvnorm_solve uses the functionality of R’s base package and is very similar in its structure to function rmvnorm_arma_inv. The Cholesky decomposition is performed using function chol and the inversion of a matrix using function solve. These functions seem to not check the type of the input matrix.

rmvnorm_solve = function(n, precision, location){
  T           = dim(precision)[1]
  precision.chol.inv    = solve(t(chol(precision)))
  
  epsilon     = matrix(rnorm(n * T), ncol=n)
  draw        = t(precision.chol.inv) %*% (matrix(rep(precision.chol.inv %*% location, n), ncol=n) + epsilon)
  
  return(draw)
}

The remaining two functions use LAPACK library R wrappers included in package mgcv by Wood (2021). This package provides dedicated functions for the Cholesky decomposition for band and tridiagonal matrices. Additionally, these implementations use R’s base package wrappers for the level-3 BLAS routine dtrsm for solving linear equations by forward- and back-substitution given by functions forwardsolve and backsolve respectively. They recognize a special type of the input matrix to be inverted to solve the linear equations.

The function rmvnorm_bandchol is written for a band precision matrix and uses function bandchol, being the implementation of LAPACK’s dpbtrf routine, to compute the Cholesky decomposition.

library(mgcv)

rmvnorm_bandchol = function(n, precision, location){
  T           = dim(precision)[1]
  precision.L = t(bandchol(precision))
  
  epsilon     = matrix(rnorm(n * T), ncol=n)
  a           = forwardsolve(precision.L, location)
  draw        = backsolve(t(precision.L), matrix(rep(a, n), ncol=n) + epsilon)
  
  return(draw)
}

The third function rmvnorm_trichol is written for a tridiagonal precision matrix and uses function trichol from package mgcv to compute the Cholesky decomposition using LAPACK’s dpttrf routine.

rmvnorm_trichol = function(n, precision, location){
  T           = dim(precision)[1]
  lead.diag   = diag(precision)
  sub.diag    = sdiag(precision, -1)
  
  precision.chol    = trichol(ld = lead.diag, sd=sub.diag)
  precision.L       = diag(precision.chol$ld)
  sdiag(precision.L,-1) = precision.chol$sd
  
  epsilon     = matrix(rnorm(n * T), ncol=n)
  a           = forwardsolve(precision.L, location)
  draw        = backsolve(t(precision.L), matrix(rep(a, n), ncol=n) + epsilon)
  
  return(draw)
}

The compiled versions of the three R functions were constructed using package compiler. However, they improved the computational time of their execution only by a small margin relative to the original functions and, thus, are skipped in the reporting below.

Another package SparseM by Koenker & Ng (2003) provides useful functions to sample random draws from the target distribution. Its application requires defining the precision matrix as sparse using function as.matrix.csr, computation of the Cholesky decomposition using function chol, and solving the linear equations using function backsolve. The computational gains from using such an implementation are significant relative to the basic R’s implementation rmvnorm_solve. However, as the package vignette points out, function chol returns a solution subjected to a random permutation of rows and columns, which makes the results uncomparable with the reproducible implementations considered here.

Performance comparison

The functions are compared using routine microbenchmark from package microbenchmark in an exercise consisting of sampling n=10 draws from the target distribution with the location matrix filled with standard normal random draws, and the precision matrix imitating the structure of a precision matrix for a state-space equation given by an autoregressive process of order 1. It implies that this matrix is tridiagonal with all the elements on the first sub-diagonal being equal to the same number. The results are reported for the precision and location of size T=250 which is motivated by the number of observations available for quarterly macroeconomic aggregates for many developed countries.

set.seed(12345)
n           = 10
T           = 250
s           = rgamma(1, shape=10, scale=10)
precision   = rgamma(1, shape=10, scale=10) * diag(T) + 2 * s * diag(T)
sdiag(precision, 1)  = -s
sdiag(precision, -1) = -s
location    = as.matrix(rnorm(T))

Function microbenchmark includes the set.seed specification and returns the tabular outcome only if the matrices sampled using alternative functions are equal to one another.

library(microbenchmark)

microbenchmark(cpp.inv      = rmvnorm_arma_inv(n, precision, location),
               cpp.solve    = rmvnorm_arma_solve(n, precision, location),
               cpp.stochvol = rmvnorm_arma_stochvol(n, precision, location),
               r.solve      = rmvnorm_solve(n, precision, location),
               r.bandchol   = rmvnorm_bandchol(n, precision, location),
               r.trichol    = rmvnorm_trichol(n, precision, location),
               check  = "equal",
               setup  = set.seed(123)
               )

Unit: microseconds
         expr      min       lq      mean   median        uq       max neval cld
      cpp.inv  790.184  882.677  3582.257  938.960  1678.659  49951.59   100  a 
    cpp.solve  504.267  914.869  1016.789  965.231  1020.201   1711.95   100  a 
 cpp.stochvol  172.530  322.851   401.402  346.279   364.765   4501.92   100  a 
      r.solve 2352.417 2583.490 17877.209 3918.437 17634.294 279008.16   100   b
   r.bandchol 1834.725 4234.104  6738.529 4432.209  5653.852  93943.39   100  a 
    r.trichol  484.612  813.264  4031.655  891.851  1605.639  97691.21   100  a

Unsurprisingly, function rmvnorm_arma_stochvol that is custom-made for this example wins the performance comparison. It is twice as fast as another cpp implementation rmvnorm_arma_solve, three times faster than the fastest R implementation rmvnorm_trichol, and over thirteen times faster than R’s rmvnorm_bandchol.

Two other unreported runs of this exercise were performed for T=750 and T=2500, the sample sizes corresponding to the number of observations available for monthly macroeconomic aggregates and ten years of daily financial data respectively. The ordering of the functions in terms of their average execution time remained the same while the multipliers increased. For instance, when T=2500 the function rmvnorm_arma_stochvol was over six times faster than both rmvnorm_arma_solve and rmvnorm_trichol and over seventeen times faster than rmvnorm_arma_inv and rmvnorm_bandchol.

Conclusion

One can expect substantial gains in the computation times by writing the functions in Rcpp for the simulation smoother. The more fine-tuned the simulation algorithm to the application at hand the more improvements are provided. The computational gains are significant in the context of using the simulation smoother in the Gibbs sampler, an iterative Bayesian estimation procedure in which a single or multiple draws from the multivariate normal distribution are sampled at each iteration, where the number of iterations ranges between several thousand up to even a million.

Tomasz thanks Alex Ballentine, Matthieu Droumaguet, and Gregor Kastner for their useful feedback.

Constructing a Sparse Matrix Class in Rcpp

Fri, 16 Apr 2021 00:00:00 +0000

Introduction

It is no secret that sparse matrix operations are faster in C++ than in R. RcppArmadillo and RcppEigen do a great job copying sparse matrices from R to C++ and back again. But note the word “copy”. Every time RcppArmadillo converts an R sparse matrix to an arma::SpMat<T> object, it has to creates a deep copy due to the difference in representation between dense matrices (usually occupying contiguous chunks of memory) and sparse matrices (which do not precisely because of sparsity). Similarly, every time RcppEigen converts an R sparse matrix to an Eigen::SparseMatrix<T> object, it also creates a deep copy.

The price of instantiating sparse matrix object is measurable in both time and memory. But one of the advantages of types such as Rcpp::NumericVector is that they simply re-use the underlying memory of the R object by wrapping around the underlying SEXP representation and its contiguous dense memory layout so no copy is created! Can we create a sparse matrix class using Rcpp::NumericVector and Rcpp::IntegerVector that uses them similarly as references rather than actual deep-copy of each element?

Structure of a `dgCMatrix`

We will focus on the dgCMatrix type, the most common form of compressed-sparse-column (CSC) matrix in the Matrix package. The dgCMatrix class is an S4 object with several slots:

library(Matrix)
set.seed(123)
str(rsparsematrix(5, 5, 0.5), vec.len = 12)

Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  ..@ i       : int [1:12] 2 4 0 3 4 0 3 4 2 3 0 2
  ..@ p       : int [1:6] 0 2 5 8 10 12
  ..@ Dim     : int [1:2] 5 5
  ..@ Dimnames:List of 2
  .. ..$ : NULL
  .. ..$ : NULL
  ..@ x       : num [1:12] 0.5 -1 0.83 -0.056 2.5 0.24 1.7 1.3 0.55 -1.7 -0.78 1.3
  ..@ factors : list()

Here we have:

Slot i is an integer vector giving the row indices of the non-zero values of the matrix
Slot p is an integer vector giving the index of the first non-zero value of each column in i
Slot x gives the non-zero elements of the matrix corresponding to rows in i and columns delineated by p

Realize that i and p are integer vectors, and x is a numeric vector (stored as a double), corresponding to Rcpp::IntegerVector and Rcpp::NumericVector. That means that we can construct a sparse matrix class in C++ using Rcpp vectors!

A dgCMatrix reference class for Rcpp

A dgCMatrix is simply an S4 object containing integer and double vectors, and Rcpp already has implemented exporters and wrappers for S4 objects, integer vectors, and numeric vectors. That makes class construction easy:

#include <Rcpp.h>

namespace Rcpp {
    class dgCMatrix {
    public:
        IntegerVector i, p, Dim;
        NumericVector x;
        List Dimnames;

        // constructor
        dgCMatrix(S4 mat) {
            i = mat.slot("i");
            p = mat.slot("p");
            x = mat.slot("x");
            Dim = mat.slot("Dim");
            Dimnames = mat.slot("Dimnames");
        };

        // column iterator
        class col_iterator {
          public:
            int index;
            col_iterator(dgCMatrix& g, int ind) : parent(g) { index = ind; }
            bool operator!=(col_iterator x) { return index != x.index; };
            col_iterator& operator++(int) { ++index; return (*this); };
            int row() { return parent.i[index]; };
            int col() { return column; };
            double& value() { return parent.x[index]; };
          private:
            dgCMatrix& parent;
            int column;
        };
        col_iterator begin_col(int j) { return col_iterator(*this, p[j]); };
        col_iterator end_col(int j) { return col_iterator(*this, p[j + 1]); };
        
    };

    template <> dgCMatrix as(SEXP mat) { return dgCMatrix(mat); }

    template <> SEXP wrap(const dgCMatrix& sm) {
        S4 s(std::string("dgCMatrix"));
        s.slot("i") = sm.i;
        s.slot("p") = sm.p;
        s.slot("x") = sm.x;
        s.slot("Dim") = sm.Dim;
        s.slot("Dimnames") = sm.Dimnames;
        return s;
    }
}

In the above code, first we create a C++ class for a dgCMatrix in the Rcpp namespace with public members corresponding to the slots in an R Matrix::dgCMatrix. Second, we add a constructor for the class that receives an S4 R object (which should be a valid dgCMatrix object). Third, we add a simple forward-only sparse column iterator with read/write access for convenient traversal of non-zero elements in the matrix. Finally, we use Rcpp::as and Rcpp::wrap for seamless conversion between R and C++ and back again.

Using the Rcpp::dgCMatrix class

We can now simply pull a dgCMatrix into any Rcpp function thanks to our handy class methods and the magic of Rcpp::as.

//[[Rcpp::export]]
Rcpp::dgCMatrix R_to_Cpp_to_R(Rcpp::dgCMatrix& mat){
    return mat;
}

And call the following from R:

library(Matrix)
spmat <- abs(rsparsematrix(100, 100, 0.1))
spmat2 <- R_to_Cpp_to_R(spmat)
all.equal(spmat, spmat2)

[1] TRUE

Awesome!

Passing by reference versus value

Rcpp structures such as NumericVector are wrappers around the existing R objects. This means that if we modify our sparse matrix in C++, it will modify the original R object. For instance:

//[[Rcpp::export]]
Rcpp::dgCMatrix Rcpp_square(Rcpp::dgCMatrix& mat){
    for (Rcpp::NumericVector::iterator i = mat.x.begin(); i != mat.x.end(); ++i)
        (*i) *= (*i);
    return mat;
}

used in

cat("sum before squaring: ", sum(spmat))

sum before squaring:  801.059

spmat2 <- Rcpp_square(spmat)
cat("sum of original object after squaring: ", sum(spmat))

sum of original object after squaring:  1007.97

cat("sum of assigned object after squaring: ", sum(spmat2))

sum of assigned object after squaring:  1007.97

Notice how the original object AND the returned object are identical, yet they don’t point to the same memory address (because of copy on write):

tracemem(spmat)

[1] "<0x55f9609a3ae0>"

tracemem(spmat2)

[1] "<0x55f960f449f0>"

You might further inspect memory addresses within the structure using .Internal(inspect()) and indeed, you will see the memory addresses are not the same. What this all means is that we can simply call the function in R and modify the object implicitly without an assignment operator.

set.seed(123)
spmat <- abs(rsparsematrix(100, 100, 0.1))
sum_before <- sum(spmat)
Rcpp_square(spmat)
sum_after <- sum(spmat)

sum_before =  793.861

sum_after  =  970.174

This principle of course applies to other Rcpp classes such as NumericVector as well. But especially when working with very large sparse matrices, it is useful to avoid deep copies and pass by reference only. If you do need to operate on a copy of your matrix in C++, there is no reason to be using an Rcpp container when you can be using RcppArmadillo or RcppEigen. These libraries are extremely capable and well-documented—and generally give you access to specific sparse Matrix algorithms.

Sparse iterator class

One of the best ways to read and/or write non-zero values in a sparse matrix is with an iterator. A basic column forward iterator with read/write access was presented in the declaration of our sparse matrix class. We can use this iterator in a very Armadillo-esque manner to compute things like column sums:

//[[Rcpp::export]]
Rcpp::NumericVector Rcpp_colSums(Rcpp::dgCMatrix& mat){
    int n_col = mat.p.size() - 1;
    Rcpp::NumericVector sums(n_col);
    for (int col = 0; col < n_col; col++)
       for (Rcpp::dgCMatrix::col_iterator it = mat.begin_col(col); it != mat.end_col(col); it++)
           sums[col] += it.value();
    return sums;
}

On the R end:

sums <- Rcpp_colSums(spmat)
all.equal(Matrix::colSums(spmat), sums)

[1] TRUE

Great—but is it faster???

library(microbenchmark)
microbenchmark(Rcpp_colSums(spmat), Matrix::colSums(spmat))

Unit: microseconds
                   expr    min      lq     mean  median      uq    max neval cld
    Rcpp_colSums(spmat)  2.510  2.7885  3.32225  3.1355  3.2685 21.842   100  a 
 Matrix::colSums(spmat) 11.068 11.5380 13.26691 11.7775 12.0525 79.324   100   b

Extending the class

There are many ways we can extend Rcpp::dgCMatrix. For example, we can use const_iterator, row_iterator and iterators that traverse all values in addition to col_iterator. We can also add support for subview copies, dense copies of columns and rows, basic element-wise operations, cross-product, etc.

We have implemented and documented these methods, and hopefully there will be a Rcpp-extending CRAN package in the not-too-distant future that allows you to seamlessly interface with the dgCMatrix class. For now, see the Github repo for RcppSparse for the in-progress package.

Using iconv via the R exported header

Sat, 06 Mar 2021 00:00:00 +0000

Introduction

Character encodings can be tricky and frustrating to deal with. Several newer languages such as Go or Julia default to native UTF-8 on all platforms, greatly facilitating and easing use of characters with languages other than English on all common platforms. With R we are not quite as lucky. UTF-8 is mostly working as desired on “operating systems with an x” but sadly, there are still a lot of Windows users out there for whom native UTF-8 is not quite in reach. A very detailed discussion of the issues involved was provided last summer on the R Developers Blog in this post.

More recently, another blog post covered character conversion. The useful idea presented in the post is to rely on the (public) header R_ext/Riconv.h which then transparently passes on to the iconv library R itself uses. (Strictly speaking this is an optional feature, see capabilities("iconv") to check your build of R.)

In order to test this, we wrapped up a little (GitHub-only) package RcppIconvExample so that we could toss the file at the accessiable Windows builders for tests (given that we do not have a physical Windows machine around). Together with an input file encoded in windows-1252 (taken from the uchardet CRAN package wrapping Mozilla’s uchardet library) we can then read and convert text in these ‘foreign’ encoding:

win1252file <- system.file("rawdata", "windows-1252.txt", package="RcppIconvExample")
win1252txt <- RcppIconvExample::read_file(win1252file, "windows-1252")
cat(win1252txt)

L’œuf de volaille est un produit agricole servant d'ingrédient entrant dans la
composition de nombreux plats, dans de nombreuses cultures gastronomiques du
monde.

Our implementation of read_file() follows. It refactors the two functions in the blog post into a single function with an optional encoding argument:

// cf https://fishandwhistle.net/post/2021/using-rs-cross-platform-iconv-wrapper-from-cpp11
std::string read_file(std::string filename, std::string encoding = "") {
    const int len = 2048;
    char buffer[len/2];

    std::ifstream file;
    file.open(filename, std::ifstream::in | std::ifstream::binary);

    file.read(buffer, len/2);
    size_t n_read = file.gcount();
    file.close();

    if (encoding == "") {       // no encoding given so return 'as is'
        return std::string(buffer, n_read);
    }

    std::string str_source(buffer, n_read);

    void* iconv_handle = Riconv_open("UTF-8", encoding.c_str());
    if (iconv_handle == ((void*) -1)) {
        Rcpp::stop("Can't convert from '%s' to 'UTF-8'", encoding.c_str());
    }

    const char* in_buffer = str_source.c_str();
    char utf8_buffer[len];
    char* utf8_buffer_mut = utf8_buffer;
    size_t in_bytes_left = n_read;
    size_t out_bytes_left = len;

    size_t result = Riconv(iconv_handle, &in_buffer, &in_bytes_left, &utf8_buffer_mut, &out_bytes_left);
    Riconv_close(iconv_handle);

    if (result == ((size_t) -1) || (in_bytes_left != 0)) {
        Rcpp::stop("Failed to convert file contents to UTF-8");
    }

    return std::string(utf8_buffer, len - out_bytes_left);
}

The entire function body is plain C++ code in a basic C++1998 standard, calls the C API of R to access iconv if a conversion is selected, and relies on Rcpp for the convenience of automating the interface and translating strings to SEXP objects and back.

The function has also been implemented in a bareboned C++ version on top of the C API of R (using the helper definitions in the tidyCpp package) in package dang as function readAndConvert().

Postscriptum: An astute reader pointed out that while iconv appears to be platform-independent, it is not. Results on Windows can differ from results on Linux. The only consistent cross-platform encoding conversions are offered by library stringi (as used by package stringr) but that one is also a known pain to install (as it is not tiny).

Using Scheme from R via Guile and Rcpp

Thu, 17 Dec 2020 00:00:00 +0000

Introduction

GNU Guile is a programming language “designed to help programmers create flexible applications that can be extended by users or other programmers with plug-ins, modules, or scripts”. It is an extension language platform, and contains an efficient compiler and virtual machine allowing both out of the box use for programs in Scheme, as well as integration with C and C++ programs.

Matías Guzmán Naranjo posed an almost-complete question on StackOverflow about how to combine Guile with R via Rcpp. It turns out that the code by Matías came from this gist showing both a Guile Scheme script, and a C use case example.

The Guile Scheme Script `script.scm`

This file is posted in this gist.

(define simple-func
  (lambda ()
    (display "Script called, now I can change this") (newline)))

(define quick-test
  (lambda ()
    (display "Adding another function, can modify without recompilation") 
    (newline)
    (adding-without-recompilation)))

(define adding-without-recompilation
  (lambda ()
    (display "Called this, without recompiling the C code") (newline)))

The C wrapper `testguile.c`

This file is also posted in this gist. We used (as shown below) simpler link instructions which (on Ubuntu) can be ontained via guile-config link and are just -lguile-3.0 -lgc on our Ubuntu 20.04 machine.

/*
 * Compile using: 
 * gcc testguile.c \
 * -I/usr/local/Cellar/guile/1.8.8/include \
 * -D_THREAD_SAFE -L/usr/local/Cellar/guile/1.8.8/lib -lguile -lltdl -lgmp -lm -lltdl -o testguile
 */

#include <stdio.h>
#include <libguile.h>

int main(int argc, char** argv) {
    SCM func, func2;
    scm_init_guile();

    scm_c_primitive_load("script.scm");

    func = scm_variable_ref(scm_c_lookup("simple-func"));
    func2 = scm_variable_ref(scm_c_lookup("quick-test"));

    scm_call_0(func);
    scm_call_0(func2);

    return 0;
}

The Rcpp Caller

This is an edited and updated version of Matías’ initial file.

#include <Rcpp.h>
#include <stdio.h>
#include <libguile.h>

using namespace Rcpp;

// [[Rcpp::export]]
int test_guile(std::string file) {
    SCM func, func2;
    scm_init_guile();

    scm_c_primitive_load(file.c_str());

    func = scm_variable_ref(scm_c_lookup("simple-func"));
    func2 = scm_variable_ref(scm_c_lookup("quick-test"));

    scm_call_0(func);
    scm_call_0(func2);

    return 0;
}

The System Integration

On our Ubuntu 20.04 box, Guile 3.0 is well supported and is easy installable. Header files are in a public location, as are the library files. The configuration helper guile-config (also guile-config-3.0 to permit continued use of 2.2 and 2.0) provides the required compiler and linker flags which we simply hardwired into src/Makevars. A proper package would need to discover these from a script, typically configure, or Makefile invocations. We leave this for another time; if anybody reading this is interested in taking this further please get in touch.

PKG_CXXFLAGS = -I"/usr/include/guile/3.0"
PKG_LIBS = -lguile-3.0 -lgc

Package

The rest of package was created via one call Rcpp.package.skeleton("RcppGuile") after which we

simply copied in the C++ file above into src/guile.cpp;
added the src/Makevars snippet above; and
copied in the example Scheme scripts as inst/guile/script.scm

and that’s it! Then any of the usual steps with R will build and install a working package.

The complete package is available in this directory.

Guile Demo

The package can be used as follows, retrieving the Scheme script from within the package and clearly running it:

> library(RcppGuile)
> test_guile(system.file("guile", "script.scm", package="RcppGuile"))
Script called, now I can change this
Adding another function, can modify without recompilation
Called this, without recompiling the C code
[1] 0
> 

GNU Guile is made for building extensions. This short note showed how easy it is to connect R with Guile via a simple Rcpp package.

Handling R6 objects in C++

Fri, 31 Jul 2020 00:00:00 +0000

Introduction

When we are using R6 objects and want to introduce some C++ code in our project, we may also want to interact with these objects using Rcpp. With this objective in mind, the key to interacting with R6 objects is that they are essentially environments, and so they can be treated as such.

In this example, we will define a ‘Person’ class with some private attributes and public methods. It will allow us to show how to create R6 objects in C++ and how to interact with them.

Creating a dummy R6 class

First, we define the class ‘Person’, with a name, an id, an item and some methods:

library(R6)

#' R6 class that defines a person
Person <<- R6::R6Class("Person",
  public = list(
    #' @description 
    #' Constructor of the 'Person' class
    #' @param name a string that defines this person's name
    #' @param id an integer that defines this person's id
    #' @return A new 'Person' object
    initialize = function(name, id){
      private$name <- name
      private$id <- id
    },
    
    get_name = function(){return(private$name)},
    
    get_id = function(){return(private$id)},
    
    #' @description 
    #' Gives an item to the person
    #' @param item a string that defines the item
    give_item = function(item){private$item <- item},
    
    #' @description 
    #' A public function that calls a private one
    good_foo = function(){
      return(paste0("Wrapped inside: {", private$bad_foo(), "}"))
    }
  ),
  private = list(
    #' @field name the name of the person
    name = NULL,
    #' @field id the id of the person
    id = NULL,
    #' @field item some item that the person has
    item = NULL,
    
    #' @description 
    #' A private function that should not be called from outside the object
    bad_foo = function(){return("This is a private function")}
  )
)

With this simple class, we will be able to create R6 objects, initialize them and call some of their methods.

Creating R6 objects in C++

To create R6 objects, we must first get the ‘new’ function that initializes them. In this case, given that we have loaded the Person class into the global environment, we will get it from there:

#include<Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
List initialize_list(StringVector &names, IntegerVector &ids){
  List res(names.size());
  Environment g_env = Environment::global_env();
  Environment p_env = g_env["Person"];
  Function new_person = p_env["new"];
  
  for(unsigned int i = 0; i < names.size(); i++){
    Environment new_p;
    String name = names[i];
    int id = ids[i];
    new_p = new_person(name, id);
    res[i] = new_p;
  }
  
  return res;
}

names <- c("Jake", "Anne", "Alex")
ids <- 1:3

res <- initialize_list(names, ids)

print(res)

[[1]]
<Person>
  Public:
    clone: function (deep = FALSE) 
    get_id: function () 
    get_name: function () 
    give_item: function (item) 
    good_foo: function () 
    initialize: function (name, id) 
  Private:
    bad_foo: function () 
    id: 1
    item: NULL
    name: Jake

[[2]]
<Person>
  Public:
    clone: function (deep = FALSE) 
    get_id: function () 
    get_name: function () 
    give_item: function (item) 
    good_foo: function () 
    initialize: function (name, id) 
  Private:
    bad_foo: function () 
    id: 2
    item: NULL
    name: Anne

[[3]]
<Person>
  Public:
    clone: function (deep = FALSE) 
    get_id: function () 
    get_name: function () 
    give_item: function (item) 
    good_foo: function () 
    initialize: function (name, id) 
  Private:
    bad_foo: function () 
    id: 3
    item: NULL
    name: Alex

The previous example initializes a list with Persons. It is important to notice that the only function relevant to us from the global environment is the ‘new’ function created by the R6 class. All the functions defined inside Person will be contained in each instance of the class.

If the R6 object is defined inside some package, included our own package if we are developing one, we must get it from there. To do this, the previous function has to be modified as follows:

#include<Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
List initialize_list(unsigned int size){
  List res(size);
  Environment package_env("package:some_package");
  Environment class_env = package_env["some_class"];
  Function new_instance = class_env["new"];
  
  for(unsigned int i = 0; i < size; i++){
    Environment new_i;
    new_i = new_instance();
    res[i] = new_i;
  }
  
  return res;
}

Calling the methods of an instance

Once we have instantiated an R6 object as an environment, we can call its methods after getting them from the instance:

#include<Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
Environment create_person_item(){
  Environment g_env = Environment::global_env();
  Environment p_env = g_env["Person"];
  Function new_person = p_env["new"];
  Environment new_p;
  String name = "Jake";
  int id = 1;
  String item = "shovel";
  
  new_p = new_person(name, id);
  Function give_i = new_p["give_item"];
  give_i(item);

  return new_p;
}

res <- create_person_item()

print(res)

<Person>
  Public:
    clone: function (deep = FALSE) 
    get_id: function () 
    get_name: function () 
    give_item: function (item) 
    good_foo: function () 
    initialize: function (name, id) 
  Private:
    bad_foo: function () 
    id: 1
    item: shovel
    name: Jake

Note that if we are creating multiple instances, we have to get the desired method from each one of them.

Private attributes and methods cannot be accessed in this manner, trying to do so results in an error:

#include<Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
String access_method(std::string foo_name){
  Environment g_env = Environment::global_env();
  Environment p_env = g_env["Person"];
  Function new_person = p_env["new"];
  Environment new_p;
  String name = "Jake";
  int id = 1;
  
  new_p = new_person(name, id);
  Function foo = new_p[foo_name];
  String res = foo();
  
  return res;
}

tryCatch({res <- access_method("bad_foo")},
         error=function(cond){print(paste0("Exception raised: ", cond))})

[1] "Exception raised: Error in access_method(\"bad_foo\"): Cannot convert object to a function: [type=NULL; target=CLOSXP, SPECIALSXP, or BUILTINSXP].\n"

res <- access_method("good_foo")

print(res)

[1] "Wrapped inside: {This is a private function}"

The error is telling us that there is no function called “bad_foo” in the environment of the object “new_p”. In that environment, only public attributes and methods are present on the surface. Theoretically, we could still try to find the private attributes and methods inside the object, but it goes against the paradigm and it would be much easier to just make the method we want to call public.

Passing down R6 objects

Passing R6 objects as arguments to C++ functions is also straight forward by treating them as environments:

#include<Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
void give_shovel(Environment &person){
  String item = "shovel";

  Function give_i = person["give_item"];
  give_i(item);
}

p <- Person$new("Jake", 1)

give_shovel(p)

print(p)

<Person>
  Public:
    clone: function (deep = FALSE) 
    get_id: function () 
    get_name: function () 
    give_item: function (item) 
    good_foo: function () 
    initialize: function (name, id) 
  Private:
    bad_foo: function () 
    id: 1
    item: shovel
    name: Jake

Note that there is no need to return the object, as the environment is passed by reference and changes inside it will be reflected outside the function.

Mixing Rcpp modules and Rcpp attributes

Mon, 08 Jun 2020 00:00:00 +0000

Introduction

With Rcpp attributes Rcpp modules (described in the Rcpp vignettes) it is easy to expose C++ classes and functions to R. This note describes how to use classes exported by modules in conjunction with functions exported using Rcpp attributes through the use of RCPP_EXPOSED* macros.

In the following snippets, a simple example is given of a simple class and a function that takes that class as an argument. The C++ function is exported using Rcpp attributes as follows:

#include <Rcpp.h>

// [[Rcpp::export]]
void shout(std::string message) {
    Rcpp::Rcout << message << "!" << std::endl;
}

Now, calling this function from R is as easy as one can hope:

shout("Hello World")

Hello World!

C++ classes can be exported using Rcpp modules. The simple class Echo below has a get() method which returns the original constructor parameter.

#include <Rcpp.h>
#include <string>

using namespace Rcpp;

class Echo {
private:
      std::string message;
public:
      Echo(std::string message) : message(message) {}

      std::string get() {
        return message;
      }
};

This class can now be exposed to R by specifing the constructors and the methods that should be callable from R with

RCPP_MODULE(echo_module) {
      class_<Echo>("Echo")
      .constructor<std::string>()
      .method("get", &Echo::get)
      ;
};

Unfortunately, combining these two snippets as above creates a problem. The Rcpp attributes machinery that exports shout() will not be automagically aware of the Echo class. This will cause an error when the package is loaded by R as the required functionality that transforms the class between a SEXP and a regular C++ object can’t be loaded. The solution is simple: instruct the compiler to do so explicitly using the RCPP_EXPOSED* family of macros. In the current case it suffices to add

RCPP_EXPOSED_AS(Echo)

Now, constructing and using the class from R is again straightforward

echo <- new(Echo, "Hello World")
echo$get()

[1] "Hello World"

shout(echo$get())

Hello World!

The `RCPP_EXPOSED*` macros

Rcpp defines a number RCPP_EXPOSED* macros in inst/include/Rcpp/macros/module.h, the most important ones are

RCPP_EXPOSED_AS which allows passing objects from R to C++. As seen above, this is needed when exported functions want to take a C++ object as argument. Other uses include methods and constructors of other Rcpp modules classes that take a C++ object as argument;
RCPP_EXPOSED_WRAP which allows the other way around; This is needed when a exported function or method wants to return a C++ object;
RCPP_EXPOSED_CLASS which allows both.

Recreating rmultinom and rpois with Rcpp

Sat, 06 Jun 2020 00:00:00 +0000

Sometimes one needs to mimic the exact behavior of R’s Distributions within C++ code. The incredible Rcpp team has provided access to these distributions through Rmath.h (in the R:: namespace), as well as through the Rcpp:: namespace where there can be two forms: scalar as in R, and vectorized via Rcpp sugar. The behavior of these functions may not always exactly match what the user expects from the standard R behavior, particularly if attempting to use the functions in Rmath.h. In particular, the functions in Rmath.h are not vectorized. In what follows, I will use Rcpp to mimic the behavior of both the rmultinom and rpois functions available in base R so that this functionality and behavior is provided in native C++.

The multinomial distribution generalizes the binomial distribution to k discrete outcomes instead of 2; consequently, it is parameterized in terms of k probabilities that must sum to 1. The base R function rmultinom used for generating multinomial data takes three arguments: n the number of simulated data sets to produce, size, the number of multinomial outcomes to sample for each data set, and prob a numeric vector of probabilities. The function returns a k $\times$ n integer matrix.

The following C++ code uses the R::rmultinom function available in Rmath.h to generate size multinomial outcomes. The R::rmultinom function relies on referencing a pointer to an IntegerVector to store the results. We create a helper function, rmultinom_1, that draws size multinomial outcomes from the multinomial distribution based on the probabilities in prob. We then do this n independent times in the function rmultinom_rcpp. To match the base R functionality, rmultinom_rcpp returns a k $\times$ n IntegerMatrix.

#include <Rcpp.h>
using namespace Rcpp;

IntegerVector rmultinom_1(unsigned int &size, NumericVector &probs, unsigned int &N) {
    IntegerVector outcome(N);
    rmultinom(size, probs.begin(), N, outcome.begin());
    return outcome;
}

// [[Rcpp::export]]
IntegerMatrix rmultinom_rcpp(unsigned int &n, unsigned int &size, NumericVector &probs) {
    unsigned int N = probs.length();
    IntegerMatrix sim(N, n);
    for (unsigned int i = 0; i < n; i++) {
        sim(_,i) = rmultinom_1(size, probs, N);
    }
    return sim;
}

We now check if the rmultinom and rmultinom_rcpp functions produce the same results. We generate a vector of 200 probabilities that sum to 1. We will sample 500 multinomial outcomes and do this independently 20 times.

prob <- runif(200)
prob <- prob/sum(prob) # standardize the probabilities
size <- 500
n <- 20

set.seed(10)
sim_r <- rmultinom(n, size, prob)
set.seed(10)
sim_rcpp <- rmultinom_rcpp(n, size, prob)
all.equal(sim_r, sim_rcpp)

[1] TRUE

A benchmark of the functions suggests that the rmultinom_rcpp function is very slightly slower than the rmultinom function, but that is not really a concern for our purposes.

microbenchmark::microbenchmark(
    rmultinom(1000, size, prob),
    rmultinom_rcpp(1000, size, prob)
)

Unit: milliseconds
                             expr     min      lq    mean  median      uq     max neval cld
      rmultinom(1000, size, prob) 10.9042 11.1841 11.7729 11.6485 12.1532 14.1841   100  a 
 rmultinom_rcpp(1000, size, prob) 11.1452 11.3780 12.0209 11.8841 12.2702 14.9434   100   b

The poisson distribution is a non-negative discrete distribution characterized by having identical mean and variance. The base R function rpois used for generating Poisson data takes two arguments: n the number of simulated values to produce, and lambda, a positive numeric vector. The rpois function cycles (and recycles) through the values in lambda for each successive value simulated. The function produces an integer vector of length n. We provide similar functionality using the R::rpois function available in Rmath.h. Note that we cycle through the values of lambda so that if the end of the lambda vector is reached before we have generated n values, then we restart at the beginning of the lambda vector.

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
IntegerVector rpois_rcpp(unsigned int &n, NumericVector &lambda) {
    unsigned int lambda_i = 0;
    IntegerVector sim(n);
    for (unsigned int i = 0; i < n; i++) {
        sim[i] = R::rpois(lambda[lambda_i]);
        // update lambda_i to match next realized value with correct mean
        lambda_i++;
        // restart lambda_i at 0 if end of lambda reached
        if (lambda_i == lambda.length()) {
            lambda_i = 0;
        }
    }
    return sim;
}

We now evaluate whether the rpois and rpois functions produce the same results. We generate a positive vector with 200 values for lambda and draw length(lambda) + 5 independent Poisson values.

lambda <- runif(200, 0.5, 3)
set.seed(10)
pois_sim_r <- rpois(length(lambda) + 5, lambda)
set.seed(10)
pois_sim_rcpp <- rpois_rcpp(length(lambda) + 5, lambda)
all.equal(pois_sim_r, pois_sim_rcpp)

[1] TRUE

A benchmark of the two functions suggests the rpois_rcpp function may be slightly faster, but once again, that is not our primary concern here.

microbenchmark::microbenchmark(
    rpois(length(lambda) + 5, lambda),
    rpois_rcpp(length(lambda) + 5, lambda)
)

Unit: microseconds
                                   expr   min     lq    mean median     uq    max neval cld
      rpois(length(lambda) + 5, lambda) 7.455 7.7825 8.02154  7.909 8.2425 11.145   100   b
 rpois_rcpp(length(lambda) + 5, lambda) 6.737 6.9860 7.31607  7.182 7.4515 16.328   100  a

Nullable Optional Arguments in Rcpp functions

Mon, 18 May 2020 00:00:00 +0000

Introduction

Often we need to have optional arguments in R of Rcpp functions with default values. Sometimes, the default value for the optional parameters is set to be NULL. Rcpp provides the Nullable <> to set default value as to be R_NilValue (equivalent of NULL in Rcpp). There have been several StackOverflow posts on using the Nullable behavior. As seen from quite a few posts, the key step in using Rcpp::Nullable<> is to cast it to the underlying type first (i.e., instantiation) after checking that the input is not NULL.

Nullability of Vector, Matrix or Logical Vector

// Checking setting Vector, Matrix and LogicalVector to NULL by default and 
// using the input if not set to NULL
#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
void nullable1(Nullable<NumericVector> NV_ = R_NilValue, 
               Nullable<NumericMatrix> NM_ = R_NilValue, 
               Nullable<LogicalVector> LG_ = R_NilValue){
  
    if (NV_.isNotNull()) {
        NumericVector NV(NV_);        // casting to underlying type NumericVector
        Rcout << "Numeric Vector is set to not NULL." << std::endl;
        Rcout << NV << std::endl;
    } else if (NM_.isNotNull()){
        NumericMatrix NM(NM_);       // casting to underlying type NumericMatrix
        Rcout << "Numeric Matrix is set to not NULL." << std::endl;
        Rcout << NM << std::endl;
    } else if (LG_.isNotNull()){
        LogicalVector LG(LG_);       // casting to underlying type Boolean
        Rcout << "Logical Vector is set to not NULL." << std::endl;
        Rcout << LG << std::endl;
    } else {
        warning("All arguments are set to NULL.\n");
    }
}

Running a few examples with setting NULL for a vector, matrix or a boolean value gives us the expected results.

nullable1(c(1,2), NULL, NULL)

Numeric Vector is set to not NULL.
1 2

m <- matrix(-0.5, 3, 3)
nullable1(NULL, m, NULL)

Numeric Matrix is set to not NULL.
-0.500000 -0.500000 -0.500000
-0.500000 -0.500000 -0.500000
-0.500000 -0.500000 -0.500000

nullable1(NULL, NULL, FALSE)

Logical Vector is set to not NULL.
0

nullable1(NULL, NULL, NULL)

Warning in nullable1(NULL, NULL, NULL): All arguments are set to NULL.

nullable1(c(), NULL, NULL)

Warning in nullable1(c(), NULL, NULL): All arguments are set to NULL.

We get the same result when the input to the NumericVector argument is not NULL but an empty vector, i.e., c(), which is also expected since is.null(c()) is TRUE in R.

A stricter test whether the input is usable can be (aptly named) isUsable().

// Testing another check, isUsable for a Nullable Vector
#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
void nullable2(Nullable<NumericVector> NV_ = R_NilValue) {
  
    if (NV_.isUsable()) {
        NumericVector NV(NV_);        // casting to underlying type NumericVector
        Rcout << "Input is usable." << std::endl;
        Rcout << NV << std::endl;
    } else {
        Rcout << "Input is either NULL or not usable." << std::endl;
    }
}

Nullability of DataFrame and List

Rcpp::Nullable<> works for SEXP based Rcpp types, so Rcpp::DataFrame and Rcpp::List can also be set to Nullable and instantiated if not NULL.

// Checking setting List and DataFrame to NULL by default and 
// using the input if not set to NULL
#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
void nullable3(Nullable<List> ls_ = R_NilValue, Nullable<DataFrame> df_ = R_NilValue){

    if (ls_.isNotNull()){
        Rcpp::List ls(ls_);                          // casting to underlying type List
        Rcout << "List is not NULL." << std::endl;
        Rcout << "List length of " << ls.length() << " elements." << std::endl;
    } else if(df_.isNotNull()) {
        Rcpp::DataFrame df(df_);                    // casting to underlying type DataFrame
        Rcout << "DataFrame is not NULL." << std::endl;
        Rcout << "DataFrame of " << df.nrows() << " rows and " << df.length() << " columns." << std::endl;
    } else {
        warning("Both inputs are NULL.\n");
    }
}

Testing with Rcpp::List and Rcpp::DataFrame gives expected results, i.e.,

mylist <- list(A = 1:10, B = letters[1:10])
nullable3(mylist, NULL)

List is not NULL.
List length of 2 elements.

df  <- data.frame(A = 1:20, B = letters[1:20])
nullable3(NULL, df)

DataFrame is not NULL.
DataFrame of 20 rows and 2 columns.

Nullability of `RcppGSL::Matrix`

In addition to Rcpp types, RcppGSL::Matrix can also be set with Nullable type (e.g, in the mvlabund package): e.g.,

// Checking setting RcppGSL Matrix to NULL by default and 
// using the input if not set to NULL
// [[Rcpp::depends(RcppGSL)]]
#include <RcppGSL.h>

using namespace Rcpp;

// [[Rcpp::export]]
void nullable4(Rcpp::Nullable<RcppGSL::Matrix> M_ = R_NilValue) {
  
    if (M_.isNotNull()){
        RcppGSL::Matrix M(M_);      // casting to underlying type RcppGSL::Matrix
        Rcout << "Input is not NULL." << std::endl;
        Rcout << "Input GSL matrix has " << M.nrow() << " and " << M.ncol() << " columns.\n";
    } else {
        warning("Input GSL Matrix is NULL.\n");
    }
}

Finally, testing with RcppGSL::Matrix which can also be set to Nullable<>, i.e.,

nullable4(NULL)  # testing with NULL

Warning in nullable4(NULL): Input GSL Matrix is NULL.

m <- matrix(-0.5, 3, 3) # testing with a non-NULL matrix
nullable4(m)

Input is not NULL.
Input GSL matrix has 3 and 3 columns.

Summary

Rcpp provides a convenient construct to set datatypes to NULL using R_NilValue and application of the datatype if not set to NULL using the .isNotNull() check. This construct to applied to set datatypes to NULL as default values and possible simple simplification.

Creating a data.table from C++

Tue, 05 May 2020 00:00:00 +0000

Introduction

Rcpp provides the DataFrame class which enables us to pass data.frame object between C++ and R. DataFrame objects are key to R and used very widely. They also provide the basis from which two key packages extend them. One of these, the tibble package, operates in a similar fashion to data.frame and treats the data as immutable. Upon a change to the data, a new version is created internally—this is frequently referred to as “copy-on-write”. Immutable data structures have some desirable properties in terms of reasoning with and about data, but the copying comes at a price in terms of performance especially once data becomes sizeable.

The other package, data.table, offers a contrasting approach. It generally modifies the data in-place; this is often referred to as by reference. As (generally) no copies of the data are needed, the runtime is often reduced. This design difference is (along with a lot of attention to performance and optimisations in general) one of the reasons why data.table is so looking rather attractive in benchmarks.

As objects of type data.table also inherit from data.frame, we can use all data.frame functions on a data.table. Internally, however, the two are very different. The question now is, how can we create a data.table efficiently inside of Rcpp without the need for a deep copy of the data?

This question is important if the data generation process takes place in C++, e.g. when parsing, simulating, or otherwise generating data. One example for which question came up is the RITCH package which parses (large) binary ITCH files (financial messages) in Rcpp and needs a fast data conversion to data.table.

Very Short Answer

We first provide a very short answer, mostly for future reference. A more detailed answer is provided below.

Essentially two steps have to be taken:

set the S3 class of the object to "data.table", "data.frame" in Rcpp

Rcpp::List create_data_table() {
    Rcpp::List res;
  
    // Populate the list
    // ...
  
    res.attr("class") = Rcpp::CharacterVector::create("data.table", "data.frame");
    return res;
}

in R call the data.table::setalloccol(df) function on the “partial” data.table

foo <- function() {
    dt <- create_data_table()
    dt <- data.table::setalloccol(dt)
    return(dt)
}

A look inside the data.table::setalloccol() functions shows that it itself calls only a few data.table functions needed to do some housekeeping and optimisation on internal state. It would be useful to access these functions directly from another package at the C level, and we plan to discuss this with the data.table team.

Detailed Answer

A first thought would be to use the Rcpp::DataFrame or the Rcpp::List class in Rcpp to return the data to R and then call data.table::setDT() on the object to properly convert it to a data.table. To give an example, lets consider creating a dataset of some data in Rcpp and returing it to R as a Data.Frame

#include<Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
DataFrame simulate_data(int n = 10, int nobs = 10000) {
    // create a list with n slots
    List res(n);
    CharacterVector list_names = CharacterVector(n);
  
    // populate the list elements
    for (int i = 0; i < n; i++) {
        list_names[i] = "col_" + std::to_string(i);
        res[i] = rnorm(nobs);
    }
  
    // set the names for the list elements
    res.names() = list_names;
    return res;
}

Using setDT(), we can convert the list to a data.table:

create_dt_naive <- function(n = 10, nobs = 10000) {
    dt <- simulate_data(n, nobs)
    data.table::setDT(dt)
    return(dt)
}
set.seed(123)
dt <- create_dt_naive()
str(dt)

Classes 'data.table' and 'data.frame':	10000 obs. of  10 variables:
 $ col_0: num  -0.5605 -0.2302 1.5587 0.0705 0.1293 ...
 $ col_1: num  2.371 -0.167 0.927 -0.568 0.225 ...
 $ col_2: num  -0.836 -0.221 -2.104 -1.668 -1.098 ...
 $ col_3: num  -0.194 0.258 -0.538 -1.179 0.901 ...
 $ col_4: num  0.4825 0.7214 -0.5078 -0.0647 1.3021 ...
 $ col_5: num  0.26 0.918 -0.722 -0.808 -0.141 ...
 $ col_6: num  -0.3883 0.0274 -0.2761 -0.0867 2.1477 ...
 $ col_7: num  0.414 -0.641 0.281 -0.694 -0.367 ...
 $ col_8: num  0.783 -0.424 -0.844 0.876 1.125 ...
 $ col_9: num  0.5732 0.0183 -0.022 -0.4278 -0.4776 ...
 - attr(*, ".internal.selfref")=<externalptr>

Inspecting the data.table code reveals that internally, data.table uses the setalloccol() function to convert a data.frame to data.table after doing some checks. This can be leveraged in the following way:

Set the class attribute of the list to "data.table", "data.frame"

#include<Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
List simulate_data_partial_dt(int n = 10, int nobs = 10000) {
    // create a list with n slots
    List res(n);
    CharacterVector list_names = CharacterVector(n);
  
    // populate the list elements
    for (int i = 0; i < n; i++) {
        list_names[i] = "col_" + std::to_string(i);
        res[i] = rnorm(nobs);
    }
    res.names() = list_names;
    // Set the class attribute of the list
    res.attr("class") = CharacterVector::create("data.table", "data.frame");
    return res;
}

and then use setalloccol() to reallocate the data by reference

create_dt_correct <- function(n = 10, nobs = 10000) {
    dt <- simulate_data_partial_dt(n, nobs)
    # reallocate the data.table by reference
    dt <- data.table::setalloccol(dt)
    return(dt)
}
set.seed(123)
dt2 <- create_dt_correct()
str(dt2)

Classes 'data.table' and 'data.frame':	10000 obs. of  10 variables:
 $ col_0: num  -0.5605 -0.2302 1.5587 0.0705 0.1293 ...
 $ col_1: num  2.371 -0.167 0.927 -0.568 0.225 ...
 $ col_2: num  -0.836 -0.221 -2.104 -1.668 -1.098 ...
 $ col_3: num  -0.194 0.258 -0.538 -1.179 0.901 ...
 $ col_4: num  0.4825 0.7214 -0.5078 -0.0647 1.3021 ...
 $ col_5: num  0.26 0.918 -0.722 -0.808 -0.141 ...
 $ col_6: num  -0.3883 0.0274 -0.2761 -0.0867 2.1477 ...
 $ col_7: num  0.414 -0.641 0.281 -0.694 -0.367 ...
 $ col_8: num  0.783 -0.424 -0.844 0.876 1.125 ...
 $ col_9: num  0.5732 0.0183 -0.022 -0.4278 -0.4776 ...
 - attr(*, ".internal.selfref")=<externalptr>

all.equal(dt, dt2)

[1] TRUE

Comparing the two methods, we see that the second implementation is a lot faster on smaller datasets.

On larger datasets, the differences is less dominant, which hints that data.table is able to use the data in place and has no need to copy it.

# small dataset
rbenchmark::benchmark(
    create_dt_naive(15, 1e3),
    create_dt_correct(15, 1e3),
    replications = 10
)

                         test replications elapsed relative user.self sys.self
2 create_dt_correct(15, 1000)           10   0.007    1.000     0.008        0
1   create_dt_naive(15, 1000)           10   0.013    1.857     0.013        0
  user.child sys.child
2          0         0
1          0         0

# larger dataset
rbenchmark::benchmark(
    create_dt_naive(15, 1e6),
    create_dt_correct(15, 1e6),
    replications = 10
)

                          test replications elapsed relative user.self sys.self
2 create_dt_correct(15, 1e+06)           10   4.202    1.000     4.035    0.167
1   create_dt_naive(15, 1e+06)           10   4.265    1.015     4.147    0.118
  user.child sys.child
2          0         0
1          0         0

R 4.0.0 brings a new function list2DF() which we could consider as well.

Note: This post draws upon, and extends, an earlier writeup at his repo.

Documenting Rcpp functions and classes in R packages

Thu, 02 Jan 2020 00:00:00 +0000

Introduction

Roxygen2 is a convenient way to document functions in R packages. Based on the Doxygen model, it parses relevant information from the comments and generates the corresponding man/*.Rd files. The major strength of this model – besides not having to write *.Rd files by hand – is that documentation ends up living right next to the functionality it describes, thus enabling easy maintenance and synchronization between the two.

When using Rcpp in package development, the R functions in RcppExports.R are generated automatically from the source .cpp files. Rcpp also faithfully transcribes comment blocks from .cpp to .R. We can take advantage of this fact to add documentation directly to the .cpp files and have it appear in the package reference manual.

This vignette demonstrates four different ways to document Rcpp functions and classes in R packages. We will create an R package from scratch, add some very basic functionality and show how Roxygen2 blocks appear in the final documentation. Throughout the vignette, we will focus on base R, Rcpp and roxygen2 with no additional dependencies. There are a number of popular add-on packages helping with package development (e.g. devtools, remotes, usethis as of January 2020), and these additional helpers may be reviewed in a follow-up vignette. In case Rcpp or roxygen2 are not yet installed, do

install.packages(c("Rcpp","roxygen2"))

Create a package

We begin by initializing a package skeleton. The following functions create a new directory and establish a barebones structure for an R package with Rcpp support:

Rcpp::Rcpp.package.skeleton("testpkg", path="~/test/")

Using your favorite text editor, create a new file ~/test/testpkg/R/zzz.R and add the following lines to it:

Rcpp::loadModule("double_cpp", TRUE)

This exposes the functionality of an Rcpp module that we will use to encapsulate our C++ class. Replace "double_cpp" with your desired module name, but note that the name has to match what appears in the .cpp file below.

Add C++ functionality and Roxygen2 comment blocks

All code blocks presented in this section are part of the same contiguous .cpp file. For the purposes of this tutorial, assume that the file lives in ~/test/testpkg/src/mult.cpp. Let’s review the four different ways to introduce Roxygen2 comment blocks inside .cpp files.

Approach 1: Standalone functions

The first approach is the most straightforward and mimics traditional usage of Roxygen2 in pure R packages. The comments are placed immediately before a stand-alone function and begin with //', which Rcpp will translate into #' in the corresponding RcppExports.R file:

#include <Rcpp.h>

//' Multiplies two doubles
//'
//' @param v1 First value
//' @param v2 Second value
//' @return Product of v1 and v2
// [[Rcpp::export]]
double mult( double v1, double v2 ) {return v1*v2;}

Approach 2: Nested field structure

Classes impose an additional layer of hierarchy by grouping multiple functions and variables. We need to decide how this hierarchy should be represented in the R help pages, which themselves follow a flat format. If a class is relatively small, a simple solution is to use a nested field structure to describe individual member functions. The entire class will then appear as a single entry in the final package reference manual.

//' @name Double
//' @title Encapsulates a double
//' @description Type the name of the class to see its methods
//' @field new Constructor
//' @field mult Multiply by another Double object \itemize{
//' \item Parameter: other - The other Double object
//' \item Returns: product of the values
//' }
class Double {
public:
  Double( double v ) : value(v) {}
  double mult( const Double& other ) {return value * other.value;}
private:
  double value;
};

Approach 3: Stand-alone pages for individual class methods

Occasionally, a member function may be complex enough to require its own ? reference entry. A major advantage of Roxygen2 is that such entries can be created with relative ease by placing the corresponding comment block at the top level in a file.

//' @name Double$new
//' @title Constructs a new Double object
//' @param v A value to encapsulate

The only disadvantage is that the documentation does not live directly next to the function it describes. Unforunately, this is due to a limitation that only the top-level comment blocks are exported into R by Rcpp.

Approach 4: Rcpp module docstrings

The final place for function documentation is inside the docstring feature provided by the Rcpp modules themselves. This works well for relatively simple classes. Unfortunately, the documentation becomes overly verbose, if a class makes heavy use of templates.

RCPP_EXPOSED_CLASS(Double)
RCPP_MODULE(double_cpp) {
  using namespace Rcpp;

  class_<Double>("Double")
    .constructor<double>("Wraps a double")
    .method("mult", &Double::mult, "Multiply by another Double object")
    ;
}

Compile, install and test

After finishing ~/test/testpkg/R/zzz.R and ~/test/testpkg/src/mult.cpp, run the following commands in an R session somewhere inside ~/test/testpkg:

Rcpp::compileAttributes()           # this updates the Rcpp layer from C++ to R
roxygen2::roxygenize(roclets="rd")  # this updates the documentation based on roxygen comments

followed by the usual R CMD build and R CMD install (or equivalent helper functions via RStudio, the usethis or devtools package or alike).

The first command compiles the .cpp code and generates R/RcppExports.R. You will notice that the roxygen comment blocks are faithfully transcribed from .cpp to .R. The second command then generates the man/*.Rd files based on these roxygen blocks. Having prepared the sources, we then build and install package as usual, making it available through the standard library(testpkg) interface.

To test the package and view its documentation, start a fresh R session and examine the help pages.

library( help=testpkg )     # Lists our functions: Double, Double$new and mult
                            #     along with those generated by Rcpp.package.skeleton()
?testpkg::mult              # Approach 1: stand-alone functions
?testpkg::Double            # Approach 2: nested field structure
?testpkg::`Double$new`      # Approach 3: individual class methods
testpkg::Double             # Approach 4: docstrings

We can also make sure that the package functions as expected.

testpkg::mult(2, 3)    
# [1] 6

d1 = testpkg::Double$new(5)
d2 = testpkg::Double$new(3)
d1$mult(d2)
# [1] 15

Going Further

As mentioned above, additional packages (or tools like the RStudio IDE) offer to help with package creation, documentation, build and more. As well, the roxygen2 can do more than we showed here by auto-generating the NAMESPACE file, collating R files as needed and more. However, we feel it helps to understand what each relevant tool offers in and by itself, and also appreciate the relative simplicity of the tools describe in the Writing R Extensions manual that is part of base R. A follow-up vignette may describe the additional tools.

Using RcppParallel to aggregate to a vector

Tue, 05 Nov 2019 00:00:00 +0000

This article demonstrates using the RcppParallel package to aggregate to an output vector. It extends directly from previous demonstrations of single-valued aggregation, through providing necessary details to enable aggregation to a vector, or by extension, to any arbitrary form.

The General Problem

Many tasks require aggregation to a vector result, and many such tasks can be made more efficient by performing such aggregation in parallel. The general problem is that the vector in which results are to be aggregated has to be shared among the parallel threads. This is a parallelReduce task - we need to split the singular task into effectively independent, parallel tasks, perform our aggregation operation on each of those tasks, yielding as many instances of our aggregate result vector as there are parallel tasks, and then finally join all of those resultant vectors from the parallel tasks into our desired singular result vector. The general structure of the code demonstrated here extends from the previous Gallery article on parallel vector sums, through extending to summation to a vector result, along with the passing of additional variables to the parallel worker. The following code demonstrates aggregation to a vector result that holds the row sums of a matrix, noting at the output that is not intended to represent efficient code, rather it is written to explicitly emphasise the principles of using RcppParallel to aggregate over a vector result.

The parallelReduce Worker

The following code defines our parallel worker, in which the input is presumed for demonstration purposes to be a matrix stored as a single vector, and so has of total length nrow * ncol. The demonstration includes a few notable features:

The main input simply provides an integer index into the rows of the matrix, with the parallel job splitting the task among elements of that index. This explicit specification of an index vector is not necessary, but serves here to clarify what the worker is actually doing. An alternative would be for input to be the_matrix, and subsequently call the parallel worker only over [0 ... nrow] of that vector which has a total length of nrow * ncol.
We are passing two additional variables specifying nrow and ncol. Although one of these could be inferred at run time, we pass them simply to demonstrate how this is done. Note in particular the form in the second constructor, called for each Split job, which accepts as input the variables as defined by the main constructor, and so all variable definitions are of the form, nrow(oneJob.nrow). The initial constructor also has input variables explicitly defined with _in suffices, to clarify exactly how such variable passing works.
No initial values for the output are passed to the constructors. Rather, output must be resized to the desired size by each of those constructors, and so each repeats the line output.resize(nrow, 0.0), which also initialises the values. (This is more readily done using a std::vector than an Rcpp vector, with final conversion to an Rcpp vector result achieved through a simple Rcpp::wrap call.)

#include <Rcpp.h>
// [[Rcpp::depends(RcppParallel)]]
#include <RcppParallel.h>
using namespace Rcpp;
using namespace RcppParallel;

struct OneJob : public Worker
{
    RVector<int> input;

    const NumericVector the_matrix;
    const size_t nrow;
    const size_t ncol;

    std::vector<double> output;

    // Constructor 1: The main constructor
    OneJob (
            const IntegerVector input_in,
            const NumericVector the_matrix_in,
            const size_t nrow_in,
            const size_t ncol_in) :
        input(input_in), the_matrix(the_matrix_in),
        nrow(nrow_in), ncol(ncol_in), output()
    {
        output.resize(nrow, 0.0);
    }

    // Constructor 2: Called for each split job
    OneJob (
            const OneJob &oneJob,
            Split) :
        input(oneJob.input), the_matrix(oneJob.the_matrix),
        nrow(oneJob.nrow), ncol(oneJob.ncol), output()
    {
        output.resize(nrow, 0.0);
    }

    // Parallel function operator
    void operator() (std::size_t begin, std::size_t end)
    {
        for (size_t i = begin; i < end; i++)
        {
            // Very inefficient yet explicit way to calculate row sums:
            for (size_t j = 0; j < ncol; j++) {
                // static_cast becuase (i,j,nrow) are size_t, aka unsigned long,
                // but Rcpp vectors require `R_xlen_t`, aka long.
                output[i] += the_matrix[static_cast<R_xlen_t>(i + j * nrow)];
            }
        }
    } // end parallel function operator

    void join (const OneJob &rhs)
    {
        for (size_t i = 0; i < nrow; i++) {
            output[i] += rhs.output[i];
        }
    }
};

The worker can then be called via parallelReduce with the following code, in which static_casts are necessary because .size() applied to Rcpp objects returns an R_xlen_t or long value, but we need to pass unsigned long or size_t values to the worker to use as indices into standard C++ vectors. The output of oneJob is a std::vector<double>, which is converted to an Rcpp::NumericVector through a simple call to Rcpp::wrap.

// [[Rcpp::export]]
NumericVector vector_aggregator (IntegerVector index, NumericVector x)
{
    const size_t nrow = static_cast <size_t> (index.size ());
    const size_t ncol = static_cast <size_t> (x.size ()) / nrow;
    OneJob oneJob (index, x, nrow, ncol);
    parallelReduce (0, nrow, oneJob);
    return wrap (oneJob.output);
}

Demonstration

Finally, the following code demonstrates that this parallel worker correctly returns the row sums of the input matrix.

# allocate a vector
nrow <- 1e5
ncol <- 10
x <- runif (nrow * ncol) # input matrix
res <- vector_aggregator (seq(nrow), x)

# confirm that this equals rowsums of the matrix:
xmat <- matrix(x, ncol = ncol)
identical(res, rowSums(xmat))

[1] TRUE

You can learn more about using RcppParallel at https://rcppcore.github.com/RcppParallel.

Performance considerations with sparse matrices in Armadillo

Wed, 18 Apr 2018 00:00:00 +0000

Introduction

Besides outstanding support for dense matrices, the Armadillo library also provides a great way to manipulate sparse matrices in C++. However, the performance characteristics of dealing with sparse matrices may be surprising if one is only familiar with dense matrices. This is a collection of observations on getting best performance with sparse matrices in Armadillo.

All the timings in this article were generated using Armadillo version 8.500. This version adds a number of substantial optimisations for sparse matrix operations, in some cases speeding things up by as much as two orders of magnitude.

General considerations: sparsity, row vs column access

Perhaps the most important thing to note is that the efficiency of sparse algorithms can depend strongly on the level of sparsity in the data. If your matrices and vectors are very sparse (most elements equal to zero), you will often see better performance even if the nominal sizes of those matrices remain the same. This isn’t specific to C++ or Armadillo; it applies to sparse algorithms in general, including the code used in the Matrix package for R. By contrast, algorithms for working with dense matrices usually aren’t affected by sparsity.

Similarly, the pattern of accesssing elements, whether by rows or by columns, can have a significant impact on performance. This is due to caching, which modern CPUs use to speed up memory access: accessing elements that are already in the cache is much faster than retrieving them from main memory. If one iterates or loops over the elements of a matrix in Armadillo, one should try to iterate over columns first, then rows, to maximise the benefits of caching. This applies to both dense and sparse matrices. (Technically, at least for dense matrices, whether to iterate over rows or columns first depends on how the data is stored in memory. Both R and Armadillo store matrices in column-major order, meaning that elements in the same column are contiguous in memory. Sparse matrices are more complex but the advice to iterate by columns is basically the same; see below.)

Matrix multiplication

We start with a simple concrete example: multiplying two matrices together. In R, this can be done using the %*% operator which (via the Matrix package) is able to handle any combination of sparse and dense inputs. However, let us assume we want to do the multiplication in Armadillo: for example if the inputs are from other C++ functions, or if we want more precise control of the output.

In Armadillo, the * operator multiplies two matrices together, and this also works for any combination of sparse and dense inputs. However, the speed of the operation can vary tremendously, depending on which of those inputs is sparse. To see this, let us define a few simple functions:

// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>

// [[Rcpp::export]]
arma::sp_mat mult_sp_sp_to_sp(const arma::sp_mat& a, const arma::sp_mat& b) {
    // sparse x sparse -> sparse
    arma::sp_mat result(a * b);
    return result;
}

// [[Rcpp::export]]
arma::sp_mat mult_sp_den_to_sp(const arma::sp_mat& a, const arma::mat& b) {
    // sparse x dense -> sparse
    arma::sp_mat result(a * b);
    return result;
}

// [[Rcpp::export]]
arma::sp_mat mult_den_sp_to_sp(const arma::mat& a, const arma::sp_mat& b) {
    // dense x sparse -> sparse
    arma::sp_mat result(a * b);
    return result;
}

The outputs of these functions are all the same, but they take different types of inputs: either two sparse matrices, or a sparse and a dense matrix, or a dense and a sparse matrix (the order matters). Let us call them on some randomly generated data:

library(Matrix)
set.seed(98765)
n <- 5e3
# 5000 x 5000 matrices, 99% sparse
a <- rsparsematrix(n, n, 0.01, rand.x=function(n) rpois(n, 1) + 1)
b <- rsparsematrix(n, n, 0.01, rand.x=function(n) rpois(n, 1) + 1)

a_den <- as.matrix(a)
b_den <- as.matrix(b)

system.time(m0 <- a %*% b)

   user  system elapsed 
  0.230   0.091   0.322

system.time(m1 <- mult_sp_sp_to_sp(a, b))

   user  system elapsed 
  0.407   0.036   0.442

system.time(m2 <- mult_sp_den_to_sp(a, b_den))

   user  system elapsed 
  1.081   0.100   1.181

system.time(m3 <- mult_den_sp_to_sp(a_den, b))

   user  system elapsed 
  0.826   0.087   0.913

all(identical(m0, m1), identical(m0, m2), identical(m0, m3))

[1] TRUE

While all the times are within an order of magnitude of each other, multiplying a dense and a sparse matrix takes about twice as long as multiplying two sparse matrices together, and multiplying a sparse and dense matrix takes about three times as long. (The sparse-dense multiply is actually one area where Armadillo 8.500 makes massive gains over previous versions. This operation used to take much longer due to using an unoptimised multiplication algorithm.)

Let us see if we can help the performance of the mixed-type functions by creating a temporary sparse copy of the dense input. This forces Armadillo to use the sparse-sparse version of matrix multiply, which as seen above is much more efficient. For example, here is the result of tweaking the dense-sparse multiply. Creating the sparse copy does take some extra time and memory, but not enough to affect the result.

// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>

// [[Rcpp::export]]
arma::sp_mat mult_sp_den_to_sp2(const arma::sp_mat& a, const arma::mat& b) {
    // sparse x dense -> sparse
    // copy dense to sparse, then multiply
    arma::sp_mat temp(b);
    arma::sp_mat result(a * temp);
    return result;
}

system.time(m4 <- mult_sp_den_to_sp2(a, b_den))

   user  system elapsed 
  0.401   0.047   0.448

identical(m0, m4)

[1] TRUE

This shows that mixing sparse and dense inputs can have significant effects on the efficiency of your code. To avoid unexpected slowdowns, consider sticking to either sparse or dense classes for all your objects. If one decides to mix them, it is worth remembering to test and profile the code.

Row vs column access

Consider another simple computation: multiply the elements of a matrix by their row number, then sum them up. (The multiply by row number is to make it not completely trivial.) That is, we want to obtain:

\[\sum_{i=1}^n{\sum_{j=1}^n{ ix_{ij} }}\]

Armadillo lets us extract individual rows and columns from a matrix, using the .row() and .col() member functions. We can use row extraction to do this computation:

// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>

// [[Rcpp::export]]
arma::sp_mat row_slice(const arma::sp_mat& x, const int n) {
    return x.row(n - 1);
}

system.time({
    result <- sapply(1:nrow(a),
        function(i) i * sum(row_slice(a, i)))
    print(sum(result))
})

[1] 1248361320

   user  system elapsed 
  1.708   0.000   1.707

For a large matrix, this takes a not-insignificant amount of time, even on a fast machine. To speed it up, we will make use of the fact that Armadillo uses the compressed sparse column (CSC) format for storing sparse matrices. This means that the data for a matrix is stored as three vectors: the nonzero elements; the row indices of these elements (ordered by column); and a vector of offsets for the first row index in each column. Since the vector of row indices is ordered by column, and we have the starting offsets for each column, it turns out that extracting a column slice is actually very fast. We only need to find the offset for that column, and then pull out the elements and row indices up to the next column offset. On the other hand, extracting a row is much more work; we have to search through the indices to find those matching the desired row.

We can put this knowledge to use on our problem. Row slicing a matrix is the same as transposing it and then slicing by columns, so let us define a new function to return a column slice. (Transposing the matrix takes only a tiny fraction of a second.)

// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>

// [[Rcpp::export]]
arma::sp_mat col_slice(const arma::sp_mat& x, const int n) {
    return x.col(n - 1);
}

a_t <- t(a)
system.time({
    result <- sapply(1:nrow(a_t),
        function(i) i * sum(col_slice(a_t, i)))
    print(sum(result))
})

[1] 1248361320

   user  system elapsed 
  0.766   0.000   0.766

The time taken has come down by quite a substantial margin. This reflects the ease of obtaining column slices for sparse matrices.

The R–C++ interface

We can take the previous example further. Each time R calls a C++ function that takes a sparse matrix as input, it makes a copy of the data. Similarly, when the C++ function returns, its sparse outputs are copied into R objects. When the function itself is very simple—as it is here—all this copying and memory shuffling can be a significant proportion of the time taken.

Rather than calling sapply in R to iterate over rows, let us do the same entirely in C++:

// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>

// [[Rcpp::export]]
double sum_by_row(const arma::sp_mat& x) {
    double result = 0;
    for (size_t i = 0; i < x.n_rows; i++) {
        arma::sp_mat row(x.row(i));
        for (arma::sp_mat::iterator j = row.begin(); j != row.end(); ++j) {
            result += *j * (i + 1);
        }
    }
    return result;
}

system.time(print(sum_by_row(a)))

[1] 1248361320

   user  system elapsed 
  0.933   0.000   0.935

This is again a large improvement. But what if we do the same with column slicing?

// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>

// [[Rcpp::export]]
double sum_by_col(const arma::sp_mat& x) {
    double result = 0;
    for (size_t i = 0; i < x.n_cols; i++) {
        arma::sp_mat col(x.col(i));
        for (arma::sp_mat::iterator j = col.begin(); j != col.end(); ++j) {
            result += *j * (i + 1);
        }
    }
    return result;
}

system.time(print(sum_by_col(a_t)))

[1] 1248361320

   user  system elapsed 
  0.005   0.000   0.006

Now the time is less than a tenth of a second, which is faster than the original code by roughly three orders of magnitude. The moral to the story is: rather than constantly switching between C++ and R, you should try to stay in one environment for as long as possible. If your code involves a loop with a C++ call inside (including functional maps like lapply and friends), consider writing the loop entirely in C++ and combine the results into a single object to return to R.

(It should be noted that this interface tax is less onerous for built-in Rcpp classes such as NumericVector or NumericMatrix, which do not require making copies of the data. Sparse data types are different, and in particular Armadillo’s sparse classes do not provide constructors that can directly use auxiliary memory.)

Iterators vs elementwise access

Rather than taking explicit slices of the data, let us try using good old-fashioned loops over the matrix elements. This is easily coded up in Armadillo, and it turns out to be quite efficient, relatively speaking. It is not as fast as using column slicing, but much better than row slicing.

// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>

// [[Rcpp::export]]
double sum_by_element(const arma::sp_mat& x) {
    double result = 0;
    // loop over columns, then rows: see comments at the start of this article
    for (size_t j = 0; j < x.n_cols; j++) {
        for (size_t i = 0; i < x.n_rows; i++) {
            result += x(i, j) * (i + 1);
        }
    }
    return result;
}

system.time(print(sum_by_element(a)))

[1] 1248361320

   user  system elapsed 
  0.176   0.000   0.176

However, we can still do better. In Armadillo, the iterators for sparse matrix classes iterate only over the nonzero elements. Let us compare this to our naive double loop:

// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>

// [[Rcpp::export]]
double sum_by_iterator(const arma::sp_mat& x) {
    double result = 0;
    for (arma::sp_mat::const_iterator i = x.begin(); i != x.end(); ++i) {
        result += *i * (i.row() + 1);
    }
    return result;
}

system.time(print(sum_by_iterator(a)))

[1] 1248361320

   user  system elapsed 
  0.002   0.000   0.002

This is the best time achieved so far, to the extent that system.time might have difficulty capturing it. The timings are now so low that we should use the microbenchmark package to get accurate measurements:

library(microbenchmark)
microbenchmark(col=sum_by_col(a_t), 
               elem=sum_by_element(a), 
               iter=sum_by_iterator(a),
               times=20)

Unit: milliseconds
 expr       min        lq      mean    median        uq       max neval
  col   4.78921   4.88444   5.05229   4.99184   5.18450   5.50579    20
 elem 172.84830 177.20431 179.87007 179.06447 182.08075 188.11256    20
 iter   1.02268   1.05447   1.12611   1.12627   1.16482   1.30800    20

Thus, using iterators represents a greater than three-order-of-magnitude speedup over the original row-slicing code.

Suppressing Call Stack Info in Rcpp-Generated Errors and Warnings

Tue, 27 Mar 2018 00:00:00 +0000

Introduction

Rcpp has an elegant mechanism of exception handling whereby C++ exceptions are automatically translated to errors in R. For most projects, the Rcpp::stop wrapper (in conjunction with the BEGIN_RCPP and END_RCPP macros automatically inserted by RcppAttributes) is sufficient and easy to use, providing an Rcpp equivalent of base::stop.

By default, it captures the call stack and attaches it to the exception in R, giving informative error messages:

#include "Rcpp.h"
using namespace Rcpp; 

//[[Rcpp::export]]         
NumericVector add1(NumericVector x, NumericVector y){
    if(x.size() != y.size()){
        stop("x and y are not the same length!");
    }
    return x + y; 
}

add1(1:5, 1:3)

Error in add1(1:5, 1:3): x and y are not the same length!

This matches the default behavior of base::stop() which captures the call info.

For complex calling patterns (e.g., creating an argument list and calling the Rcpp function with do.call), the resulting error messages are less helpful:

#include "Rcpp.h"
using namespace Rcpp; 

// [[Rcpp::export]]
NumericVector internal_function_name(NumericVector x, NumericVector y){
    if(x.size() != y.size()){
        stop("x and y are not the same length!");
    }
    return x + y; 
}

add2 <- function(x, y){
    if(!is.numeric(x)){
        x <- as.numeric(x)
    }

    do.call(internal_function_name, list(x, y))
}

add2(1:5, 1:3)

Error in (function (x, y) : x and y are not the same length!

If the internal error were being generated in R code, we might choose to use the call.=FALSE argument to base::stop to suppress the unhelpful (function (x, y) part of the error message, but we don’t (immediately) have a corresponding option in Rcpp. In this gallery post, we show how to suppress the call-stack capture of Rcpp::stop to give cleaner error messages.

Error Messages

The key functionality was added to Rcpp by Jim Hester in Rcpp Pull Request #663. To generate an R-level exception without a call stack, we pass an optional false flag to Rcpp::exception. For example,

#include "Rcpp.h"
using namespace Rcpp; 

// [[Rcpp::export]]
NumericVector internal_function_name2(NumericVector x, NumericVector y){
    if(x.size() != y.size()){
        throw Rcpp::exception("x and y are not the same length!", false);
    }
    return x + y; 
}

add3 <- function(x, y){
    if(!is.numeric(x)){
        x <- as.numeric(x)
    }

    do.call(internal_function_name2, list(x, y))
}

add3(1:5, 1:3)

Error: x and y are not the same length!

This can’t capture the R level call stack, but it is at least cleaner than the error message from the previous example.

Note that here, as elsewhere in C++, we need to handle exceptions using a try/catch structure, but we do not add it explicitly because RcppAttributes automatically handles this for us.

Warnings

Similar to Rcpp::stop, Rcpp also provides a warning function to generate R level warnings. It has the same call-stack capture behavior as stop.

For the direct call case:

#include "Rcpp.h"
using namespace Rcpp; 

//[[Rcpp::export]]         
NumericVector add4(NumericVector x, NumericVector y){
    if(x.size() != y.size()){
        warning("x and y are not the same length!");
    }
    return x + y; 
}

add4(1:5, 1:3)

Warning in add4(1:5, 1:3): x and y are not the same length!

[1] 2 4 6 4 5

For the indirect call case:

#include "Rcpp.h"
using namespace Rcpp; 

// [[Rcpp::export]]
NumericVector internal_function_name3(NumericVector x, NumericVector y){
    if(x.size() != y.size()){
        warning("x and y are not the same length!");
    }
    return x + y; 
}

add5 <- function(x, y){
    if(!is.numeric(x)){
        x <- as.numeric(x)
    }

    do.call(internal_function_name3, list(x, y))
}

add5(1:5, 1:3)

Warning in (function (x, y) : x and y are not the same length!

[1] 2 4 6 4 5

If we want to suppress the call stack info in this warning, we have to drop down to the C-level R API. In particular, we use the Rf_warningcall function, which takes the call as the first argument. By passing a NULL, we suppress the call:

#include "Rcpp.h"
using namespace Rcpp; 

// [[Rcpp::export]]
NumericVector internal_function_name5(NumericVector x, NumericVector y){
    if(x.size() != y.size()){
        Rf_warningcall(R_NilValue, "x and y are not the same length!");
    }
    return x + y; 
}

add6 <- function(x, y){
    if(!is.numeric(x)){
        x <- as.numeric(x)
    }

    do.call(internal_function_name5, list(x, y))
}

add6(1:5, 1)

Warning: x and y are not the same length!

[1] 2 2 3 4 5

A C++11 Implementation

The above methods work, but they are not as clean as their Rcpp::stop and Rcpp::warning counterparts. We can take advantage of C++11 to provide similar functionality for our call-free versions.

Basing our implementation on the C++11 implementation of Rcpp::stop and Rcpp::warning we can define our own stopNoCall and warningNoCall

#include "Rcpp.h"
using namespace Rcpp; 

// [[Rcpp::plugins(cpp11)]]
template <typename... Args>
inline void warningNoCall(const char* fmt, Args&&... args ) {
    Rf_warningcall(R_NilValue, tfm::format(fmt, std::forward<Args>(args)... ).c_str());
}

template <typename... Args>
inline void NORET stopNoCall(const char* fmt, Args&&... args) {
    throw Rcpp::exception(tfm::format(fmt, std::forward<Args>(args)... ).c_str(), false);
}

// [[Rcpp::export]]
NumericVector internal_function_name6(NumericVector x, NumericVector y, bool warn){
    if(x.size() != y.size()){
        if(warn){
            warningNoCall("x and y are not the same length!");  
        } else {
            stopNoCall("x and y are not the same length!");
        }
        
    }
    return x + y; 
}

add7 <- function(x, y, warn=TRUE){
    if(!is.numeric(x)){
        x <- as.numeric(x)
    }

    do.call(internal_function_name6, list(x, y, warn))
}

add7(1:5, 1:3, warn=TRUE)

Warning: x and y are not the same length!

[1] 2 4 6 4 5

add7(1:5, 1:3, warn=FALSE)

Error: x and y are not the same length!

Note that we used C++11 variadic templates here – if we wanted to do something similar in C++98, we could use essentially the same pattern, but would need to implement each case individually.

Introducing RcppArrayFire

Wed, 07 Mar 2018 00:00:00 +0000

Introduction

The RcppArrayFire package provides an interface from R to and from the ArrayFire library, an open source library that can make use of GPUs and other hardware accelerators via CUDA or OpenCL.

The official R bindings expose ArrayFire data structures as S4 objects in R, which would require a large amount of code to support all the methods defined in ArrayFire’s C/C++ API. RcppArrayFire instead, which is derived from RcppFire by Kazuki Fukui, follows the lead of packages like RcppArmadillo or RcppEigen to provide seamless communication between R and ArrayFire at the C++ level.

Installation

Please note that RcppArrayFire is developed and tested on Linux systems. There is preliminary support for Mac OS X.

In order to use RcppArrayFire you will need the ArrayFire library and header files. While ArrayFire has been packaged for Debian, I currently prefer using upstream’s binary installer or building from source.

RcppArrayFire is not on CRAN, but you can install the current version via drat:

#install.packages("drat")
drat::addRepo("daqana")
install.packages("RcppArrayFire")

If you have installed ArrayFire in a non-standard directory, you have to use the configure argument --with-arrayfire, e.g.:

install.packages("RcppArrayFire", configure.args = "--with-arrayfire=/opt/arrayfire-3")

A first example

Let’s look at the classical example of calculating $\pi$ via simulation. The basic idea is to generate a large number of random points within the unit square. An approximation for $\pi$ can then be calculated from the ratio of points within the unit circle to the total number of points. A vectorized implementation in R might look like this:

piR <- function(N) {
    x <- runif(N)
    y <- runif(N)
    4 * sum(sqrt(x^2 + y^2) < 1.0) / N
}

set.seed(42)
system.time(cat("pi ~= ", piR(10^6), "\n"))

pi ~=  3.13999

   user  system elapsed 
  0.047   0.008   0.055

A simple way to use C++ code in R is to use the inline package or cppFunction() from Rcpp, which are both possible with RcppArrayFire. An implementation in C++ using ArrayFire might look like this:

src <- '
double piAF (const int N) {
    array x = randu(N, f32);
    array y = randu(N, f32);
    return 4.0 * sum<float>(sqrt(x*x + y*y) < 1.0) / N;
}'
Rcpp::cppFunction(code = src, depends = "RcppArrayFire", includes = "using namespace af;")

RcppArrayFire::arrayfire_set_seed(42)
cat("pi ~= ", piAF(10^6), "\n") # also used for warm-up

pi ~=  3.1451

system.time(piAF(10^6))

   user  system elapsed 
  0.001   0.000   0.001

Several things are worth noting:

The syntax is almost identical. Besides the need for using types and a different function name when generating random numbers, the argument f32 to randu as well as the float type catches the eye. These instruct ArrayFire to use single precision floats, since not all devices support double precision floating point numbers. If you want and can to use double precision, you have to specify f64 and double.
The results are not the same since ArrayFire uses a different random number generator.
The speed-up can be quite impressive. However, the first invocation of a function is often not as fast as expected due to the just-in-time compilation used by ArrayFire. This can be circumvented by using a warm-up call with (normally) fewer computations.

Arrays as parameters

Up to now we have only considered simple types like double or int as function parameters and return values. However, we can also use arrays. Consider the case of an European put option that was recently handled with R, Rcpp and RcppArmadillo. The Armadillo based function from this post reads:

#include <RcppArmadillo.h>

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

using arma::colvec;
using arma::log;
using arma::normcdf;
using std::sqrt;
using std::exp;


// [[Rcpp::export]]
colvec put_option_pricer_arma(colvec s, double k, double r, double y, double t, double sigma) {
  
  colvec d1 = (log(s / k) + (r - y + sigma * sigma / 2.0) * t) / (sigma * sqrt(t));
  colvec d2 = d1 - sigma * sqrt(t);
  
  // Notice the use of % to represent element wise multiplication
  colvec V = normcdf(-d2) * k * exp(-r * t) - s * exp(-y * t) % normcdf(-d1); 
  
  return V;
}

This function can be applied to a range of spot prices:

put_option_pricer_arma(s = 55:60, 60, .01, .02, 1, .05)

        [,1]
[1,] 5.52021
[2,] 4.58142
[3,] 3.68485
[4,] 2.85517
[5,] 2.11883
[6,] 1.49793

Porting this code to RcppArrayFire is straight forward:

#include <RcppArrayFire.h>

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

using af::array;
using af::log;
using af::erfc;
using std::sqrt;
using std::exp;

array normcdf(array x) {
  return erfc(-x / sqrt(2.0)) / 2.0;
}

// [[Rcpp::export]]
array put_option_pricer_af(RcppArrayFire::typed_array<f32> s, double k, double r,
                           double y, double t, double sigma) {
  
  array d1 = (log(s / k) + (r - y + sigma * sigma / 2.0) * t) / (sigma * sqrt(t));
  array d2 = d1 - sigma * sqrt(t);
  
  return normcdf(-d2) * k * exp(-r * t) - s * exp(-y * t) * normcdf(-d1); 
}

Compared with the implementations in R, Rcpp and RcppArmadillo the syntax is again almost the same. One exception is that ArrayFire does not contain a function for the cumulative normal distribution function. However, the closely related error function is available.

Since an object of type af::array can contain different data types, the templated wrapper class RcppArrayFire::typed_array<> is used to indicate the desired data type when converting from R to C++. Again single precision floats are used with ArrayFire, which leads to differences of the order $10^{-6}$ compared to the results from R, Rcpp and RcppArmadillo:

put_option_pricer_af(s = 55:60, 60, .01, .02, 1, .05)

[1] 5.52021 4.58142 3.68485 2.85517 2.11883 1.49793

Performance

The reason to use hardware accelerators is of course the quest for increased performance. How does ArrayFire fare in this respect? Using the same benchmark as in the R, Rcpp and RcppArmadillo comparison:

s <- matrix(seq(0, 100, by = .0001), ncol = 1)
rbenchmark::benchmark(Arma = put_option_pricer_arma(s, 60, .01, .02, 1, .05),
                      AF = put_option_pricer_af(s, 60, .01, .02, 1, .05),
                      order = "relative", 
                      replications = 100)[,1:4]

  test replications elapsed relative
2   AF          100   0.489    1.000
1 Arma          100   1.928    3.943

Here a Nvidia GeForce GT 1030 is used together with ArrayFire’s CUDA backend. With a build-in Intel HD Graphics 520 using the OpenCL backend the ArrayFire solution is about 6 times faster. Even without a high performance GPU the performance boost from using ArrayFire can be quite impressive. However, the results change dramatically, if fewer options are evaluated:

s <- matrix(seq(0, 100, by = 1), ncol = 1)
# use more replications to get run times of more than 10 ms
rbenchmark::benchmark(Arma = put_option_pricer_arma(s, 60, .01, .02, 1, .05),
                      AF = put_option_pricer_af(s, 60, .01, .02, 1, .05),
                      order = "relative", 
                      replications = 1000)[,1:4]

  test replications elapsed relative
1 Arma         1000   0.007    1.000
2   AF         1000   0.088   12.571

But is it realistic to process $10^6$ options at once? Probably not in the way used in the benchmark where only the spot price is allowed to vary. However, one can alter the function to process not only arrays of spot prices but also arrays of strikes, risk free rates etc.:

#include <RcppArrayFire.h>

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

using af::array;
using af::log;
using af::erfc;
using af::sqrt; // arrayfire function instead of standard function
using af::exp;  // arrayfire function instead of standard function

array normcdf(array x) {
  return erfc(-x / sqrt(2.0)) / 2.0;
}

// [[Rcpp::export]]
array put_option_pricer_af(RcppArrayFire::typed_array<f32> s,
                           RcppArrayFire::typed_array<f32> k, 
                           RcppArrayFire::typed_array<f32> r, 
                           RcppArrayFire::typed_array<f32> y, 
                           RcppArrayFire::typed_array<f32> t, 
                           RcppArrayFire::typed_array<f32> sigma) {
  
  array d1 = (log(s / k) + (r - y + sigma * sigma / 2.0) * t) / (sigma * sqrt(t));
  array d2 = d1 - sigma * sqrt(t);
  
  return normcdf(-d2) * k * exp(-r * t) - s * exp(-y * t) * normcdf(-d1); 
}

Note that ArrayFire does not recycle elements if arrays with non-matching dimensions are combined. In this particular case this means that all arrays must have the same length. One can ensure that by using a data frame for the values:

set.seed(42)
# 1000 * 21 * 3 * 3 * 3 * 3 = 1,701,000 different options
options <- expand.grid(
  s = rnorm(1000, mean = 60, sd = 20),
  k = 50:70,
  r = c(0.01, 0.005, 0.02),
  y = c(0.02, 0.01, 0.04),
  t = c(1, 0.5, 2),
  sigma = c(0.05, 0.025, 0.1)
)
head(within(options,
            p <- put_option_pricer_af(s, k, r, y, t, sigma)))

        s  k    r    y t sigma           p
1 87.4192 50 0.01 0.02 1  0.05 7.45937e-29
2 48.7060 50 0.01 0.02 1  0.05 2.09402e+00
3 67.2626 50 0.01 0.02 1  0.05 2.34672e-09
4 72.6573 50 0.01 0.02 1  0.05 6.84430e-14
5 68.0854 50 0.01 0.02 1  0.05 5.26394e-10
6 57.8775 50 0.01 0.02 1  0.05 2.57756e-03

Conclusion

The ArrayFire library provides a convenient way to use hardware accelerators without the need to write low-level OpenCL or CUDA code. The C++ syntax is actually quite similar to properly vectorized R code. The RcppArrayFire package makes this available to useRs. However, one still has to be careful: using hardware accelerators is not a “silver bullet” due to the inherent memory transfer overhead.

Using RcppArmadillo to price European Put Options

Wed, 28 Feb 2018 00:00:00 +0000

Introduction

In the quest for ever faster code, one generally begins exploring ways to integrate C++ with R using Rcpp. This post provides an example of multiple implementations of a European Put Option pricer. The implementations are done in pure R, pure Rcpp using some Rcpp sugar functions, and then in Rcpp using RcppArmadillo, which exposes the incredibly powerful linear algebra library, Armadillo.

Under the Black-Scholes model The value of a European put option has the closed form solution:

\[V = K e^{-rt} N(-d_2) - S e^{-yt} N(-d_1)\]

where

\[\begin{equation} \begin{aligned} V &= \text{Value of the option} \\ r &= \text{Risk free rate} \\ y &= \text{Dividend yield} \\ t &= \text{Time to expiry} \\ S &= \text{Current stock price} \\ K &= \text{Strike price} \\ N(.) &= \text{Normal CDF} \end{aligned} \end{equation}\]

and

\[\begin{equation} \begin{aligned} d_1 &= \frac{log(\frac{S}{K}) + (r - y + \frac{1}{2} \sigma^2)t}{\sigma \sqrt{t}} \\ d_2 &= d_1 - \sigma \sqrt{t}\\ \end{aligned} \end{equation}\]

Armed with the formulas, we can create the pricer using just R.

put_option_pricer <- function(s, k, r, y, t, sigma) {

    d1 <- (log(s / k) + (r - y + sigma^2 / 2) * t) / (sigma * sqrt(t))
    d2 <- d1 - sigma * sqrt(t)

    V <- pnorm(-d2) * k * exp(-r * t) - s * exp(-y * t) * pnorm(-d1)

    V
}

# Valuation with 1 stock price
put_option_pricer(s = 55, 60, .01, .02, 1, .05)

[1] 5.52021

# Valuation across multiple prices
put_option_pricer(s = 55:60, 60, .01, .02, 1, .05)

[1] 5.52021 4.58142 3.68485 2.85517 2.11883 1.49793

Let’s see what we can do with Rcpp. Besides explicitely stating the types of the variables, not much has to change. We can even use the sugar function, Rcpp::pnorm(), to keep the syntax as close to R as possible. Note how we are being explicit about the symbols we import from the Rcpp namespace: the basic vector type, and well the (vectorized) ‘Rcpp Sugar’ calls log() and pnorm() calls. Similarly, we use sqrt() and exp() for the calls on an atomic double variables from the C++ Standard Library. With these declarations the code itself is essentially identical to the R code (apart of course from requiring both static types and trailing semicolons).

#include <Rcpp.h>
                                        
using Rcpp::NumericVector;
using Rcpp::log;
using Rcpp::pnorm;
using std::sqrt;
using std::log;

// [[Rcpp::export]]
NumericVector put_option_pricer_rcpp(NumericVector s, double k, double r, double y, double t, double sigma) {

    NumericVector d1 = (log(s / k) + (r - y + sigma * sigma / 2.0) * t) / (sigma * sqrt(t));
    NumericVector d2 = d1 - sigma * sqrt(t);
    
    NumericVector V = pnorm(-d2) * k * exp(-r * t) - s * exp(-y * t) * pnorm(-d1);
    return V;
}

We can call this from R as well:

# Valuation with 1 stock price
put_option_pricer_rcpp(s = 55, 60, .01, .02, 1, .05)

[1] 5.52021

# Valuation across multiple prices
put_option_pricer_rcpp(s = 55:60, 60, .01, .02, 1, .05)

[1] 5.52021 4.58142 3.68485 2.85517 2.11883 1.49793

Finally, let’s look at RcppArmadillo. Armadillo has a number of object types, including mat, colvec, and rowvec. Here, we just use colvec to represent a column vector of prices. By default in Armadillo, * represents matrix multiplication, and % is used for element wise multiplication. We need to make this change to element wise multiplication in 1 place, but otherwise the changes are just switching out the types and the sugar functions for Armadillo specific functions.

Note that the arma::normcdf() function is in the upcoming release of RcppArmadillo, which is 0.8.400.0.0 at the time of writing and still in CRAN’s incoming. It also requires the C++11 plugin.

#include <RcppArmadillo.h>

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

using arma::colvec;
using arma::log;
using arma::normcdf;
using std::sqrt;
using std::log;


// [[Rcpp::export]]
colvec put_option_pricer_arma(colvec s, double k, double r, double y, double t, double sigma) {
  
    colvec d1 = (log(s / k) + (r - y + sigma * sigma / 2.0) * t) / (sigma * sqrt(t));
    colvec d2 = d1 - sigma * sqrt(t);
    
    // Notice the use of % to represent element wise multiplication
    colvec V = normcdf(-d2) * k * exp(-r * t) - s * exp(-y * t) % normcdf(-d1); 

    return V;
}

Use from R:

# Valuation with 1 stock price
put_option_pricer_arma(s = 55, 60, .01, .02, 1, .05)

        [,1]
[1,] 5.52021

# Valuation across multiple prices
put_option_pricer_arma(s = 55:60, 60, .01, .02, 1, .05)

        [,1]
[1,] 5.52021
[2,] 4.58142
[3,] 3.68485
[4,] 2.85517
[5,] 2.11883
[6,] 1.49793

Finally, we can run a speed test to see which comes out on top.

s <- matrix(seq(0, 100, by = .0001), ncol = 1)

rbenchmark::benchmark(R = put_option_pricer(s, 60, .01, .02, 1, .05),
                      Arma = put_option_pricer_arma(s, 60, .01, .02, 1, .05),
                      Rcpp = put_option_pricer_rcpp(s, 60, .01, .02, 1, .05), 
                      order = "relative", 
                      replications = 100)[,1:4]

  test replications elapsed relative
2 Arma          100   6.409    1.000
3 Rcpp          100   7.917    1.235
1    R          100   9.091    1.418

Interestingly, Armadillo comes out on top here on this (multi-core) machine (as Armadillo uses OpenMP where available in newer versions). But the difference is slender, and there is certainly variation in repeated runs. And the nicest thing about all of this is that it shows off the “embarassment of riches” that we have in the R and C++ ecosystem for multiple ways of solving the same problem.

Custom progress bars for RcppProgress

Thu, 28 Dec 2017 00:00:00 +0000

RcppProgress is a tool to help you monitor the execution time of your C++ code, by providing a way to interrupt the execution inside the C++ code, and also to display a progress bar indicative of the state of your computation. Additionally, it is compatible with multi-threaded code, for example using OpenMP. The initial (yet updated) article explains the basic setup.

Since version 0.4 it became more simple to create custom progress bars. In this new article we will show how to do this. Our final example displays a progress bar which provides an estimation of the remaining time (ETA) to finish a computation.

A minimal example

Imagine you added a progress bar with RcppProgress to your function long_computation() following the example from the first article mentioned above.

// [[Rcpp::depends(RcppProgress)]]
#include <progress.hpp>
#include <progress_bar.hpp>
// [[Rcpp::export]]
double long_computation(int nb, bool display_progress=true) {
    double sum = 0;
    Progress p(nb*nb, display_progress);
    for (int i = 0; i < nb; ++i) {
        if (Progress::check_abort() )
            return -1.0;
        for (int j = 0; j < nb; ++j) {
            p.increment(); 
            sum += R::dlnorm(i+j, 0.0, 1.0, 0);
        }
    }
    return sum + nb;
}

long_computation(10)

[1] 12.2038

What you get is a basic and useful console visualization that looks like this:

0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
******************************

That’s the default, platform independent display in RcppProgress defined in SimpleProgressBar.hpp. It’s OK for most purposes to give you an idea how much work is done and it also allows you to make a very intuitive estimation about how long it’s going to take to finish. But of course that’s not everything a progress bar could show you. A progress bar could give you information about the running progress or about performance parameters of your system. It could contain calculated estimates of passed and remaining time. After all it could just look much more fancy to impress your colleagues.

RcppProgress makes it now easy to create your own implementation of a progress bar class. Your own class has to be derived from the abstract class ProgressBar that defines some basic virtual methods:

class ProgressBar {
  public:
    virtual ~ProgressBar() = 0;
    virtual void display() = 0;
    virtual void update(float progress) = 0;
    virtual void end_display() = 0;
};

display() starts the display that will be updated by subsequent calls of update(). end_display finalizes it. Your progress bar implementation should not rely on the destructor to finalize the display.

A very minimal setup could look something like this:

#include <R_ext/Print.h>

// [[Rcpp::depends(RcppProgress)]]
#include <progress.hpp>
#include "progress_bar.hpp"

class MinimalProgressBar: public ProgressBar{
  public:
    MinimalProgressBar()  {
        _finalized = false;
    }
    
    ~MinimalProgressBar() {}
    
    void display() {
      REprintf("Progress: ");
    }
    
    void update(float progress) {
      if (_finalized) return;
      REprintf("+");
    }
    
    void end_display() {
      if (_finalized) return;
      REprintf("\n");
      _finalized = true;
    }
    
  private:
  
    bool _finalized;
          
};

// [[Rcpp::export]]
double long_computation_minimal(int nb, bool display_progress=true) {
    MinimalProgressBar pb;
    double sum = 0;
    Progress p(nb, display_progress, pb);
    for (int i = 0; i < nb; ++i) {
        if (Progress::check_abort() )
            return -1.0;
        for (int j = 0; j < nb; ++j) {
            sum += R::dlnorm(i+j, 0.0, 1.0, 0);
        }
        p.increment(); 
    }
    return sum + nb;
}

long_computation_minimal(10)

[1] 12.2038

The display() method in this example does nothing more than printing the word Progress. update() concatenates a + symbol every time Progress::increment() is called. The result looks like this:

Progress: ++++++++++

In comparison to the example of the default progress bar above, I moved the call to increment() out of the second level and into the first level loop to keep the amount of console output at bay. update() also checks if the display is _finalized. end_display triggers the finalization.

Remaining time estimation

Based on the minimal setup above, you can implement more sophisticated progress bars. Here’s an example of one that looks exactly like the default SimpleProgressBar, but adds an estimation of the remaining time for the process to finish. You can find a complete package setup with the code for this ETAProgressBar here. In this article we only highlight some crucial aspects of the implementation.

We use the Rinterface.h header to update the display dynamically. Unfortunately this header is only available for Unix-like systems. A less cool, old version of an ETA progress bar that also works on windows can be found here. The following preprocessor statements load Rinterface.h if the code is compiled on a non-windows computer.

#if !defined(WIN32) && !defined(__WIN32) && !defined(__WIN32__)
#include <Rinterface.h>
#endif

The class ETAProgressBar inherits from the abstract class ProgressBar. It has an integer variable _max_ticks that controls the amount of individual tick symbols necessary to reach the 100% mark of the progress bar. That depends on the display you want to craft. ETAProgressBar also has a boolean flag variable _timer_flag that acts as a switch to separate the initial starting turn where the time measurement starts and the following turns where the time is picked off. The measured time values are stored in two variables start and end of class time_t (from ctime).

class ETAProgressBar: public ProgressBar{

  // ...
  
  private: 
    int _max_ticks;
    bool _finalized;
    bool _timer_flag;
    time_t start,end;
    
  // ...  
    
}

The display() function initializes the progress bar visualization. The first two lines are hard coded ASCII art.

void display() {
  REprintf("0%%   10   20   30   40   50   60   70   80   90   100%%\n");
  REprintf("[----|----|----|----|----|----|----|----|----|----|\n");
  flush_console();
}

update() is the most important function for the progress bar mechanism. The if clause allows to separate the initial call of update() from the following ones to start the time counter. Afterwards the time passed is calculated and transformed to a human readable string by the custom function _time_to_string(). _current_ticks_display() is another custom function to transform the progress information to a string with the correct amount of * symbols and filling whitespaces. The progress string and the time string are concatenated to create the additional third line below the initial two lines drawn by display(). A string with sufficient whitespaces is also added to ensure that this dynamically updated line is overwritten completely from turn to turn. REprintf("\r"); triggers a carriage return to make this continuous overwriting possible.

void update(float progress) {
  
    // stop if already finalized
    if (_finalized) return;
  
    // start time measurement when update() is called the first time
    if (_timer_flag) {
        _timer_flag = false;
        // measure start time
        time(&start);
    } else {
    
        // measure current time
        time(&end);
    
        // calculate passed time and remaining time (in seconds)
        double pas_time = std::difftime(end, start);
        double rem_time = (pas_time / progress) * (1 - progress);
    
        // convert seconds to time string
        std::string time_string = _time_to_string(rem_time);
    
        // create progress bar string 
        std::string progress_bar_string = _current_ticks_display(progress);
    
        // ensure overwriting of old time info
        int empty_length = time_string.length();
        std::string empty_space = std::string(empty_length, ' ');
    
        // merge progress bar and time string
        std::stringstream strs;
        strs << "|" << progress_bar_string << "| " << time_string << empty_space;
        std::string temp_str = strs.str();
        char const* char_type = temp_str.c_str();
    
        // print: remove old display and replace with new
        REprintf("\r");
        REprintf("%s", char_type);
    
        // finalize display when ready
        if(progress == 1) {
           _finalize_display();
        }  
    }
}

_time_to_string() parses time information given in form of a floating point number of seconds to a human-readable string. The basic algorithm is based on an example from programmingnotes.org.

std::string _time_to_string(double seconds) {
  
    int time = (int) seconds;
  
    int hour = 0;
    int min = 0;
    int sec = 0;
  
    hour = time / 3600;
    time = time % 3600;
    min = time / 60;
    time = time % 60;
    sec = time;
  
    std::stringstream time_strs;
    if (hour != 0) time_strs << hour << "h ";
    if (min != 0) time_strs << min << "min ";
    if (sec != 0) time_strs << sec << "s ";
    std::string time_str = time_strs.str();
  
    return time_str;
}

_current_ticks_display() relies on _compute_nb_ticks() to first of all transform the progress information (floating point number between 0 and 1) to a natural number that expresses the correct fraction of _max_ticks. _construct_ticks_display_string() takes this value and parses a string with * symbols and whitespaces that can be plotted as a visual progress indication.

std::string _current_ticks_display(float progress) {
    int nb_ticks = _compute_nb_ticks(progress);
    std::string cur_display = _construct_ticks_display_string(nb_ticks);
    return cur_display;
}

int _compute_nb_ticks(float progress) {
    return int(progress * _max_ticks);
}

std::string _construct_ticks_display_string(int nb) {
    std::stringstream ticks_strs;
    for (int i = 0; i < (_max_ticks - 1); ++i) {
        if (i < nb) {
            ticks_strs << "*";
        } else {
            ticks_strs << " ";
        }
    }
    std::string tick_space_string = ticks_strs.str();
    return tick_space_string;
}

flush_console() is a wrapper around R_FlushConsole() which is called to flush any pending output to the system console. It’s necessary to do this when the display is started in display() and when it’s closed in _finalize_display().

void flush_console() {
#if !defined(WIN32) && !defined(__WIN32) && !defined(__WIN32__)
          R_FlushConsole();
#endif
}

The output of an ETAProgressBar looks like this:

0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
|*******                                          | 49s

Define a custom print method for exposed C++ classes

Tue, 19 Dec 2017 00:00:00 +0000

When writing an R package wrapping some C++ data structures, using Rcpp Modules is a convenient option. After exposing a class to R, it can be used to easily create new instances of that class.

As an example, let us look at the Uniform class of the Rcpp Modules vignette:

// cf Rcpp Modules vignette

#include <Rcpp.h>

using namespace Rcpp;

class Uniform {
  public:
    Uniform(double min_, double max_) : min(min_), max(max_) {}
    NumericVector draw(int n) const {
      RNGScope scope;
      return runif( n, min, max );
    }
    double min, max;
};

double uniformRange( Uniform* w) {
  return w->max - w->min;
}

RCPP_MODULE(unif_module) {
  class_<Uniform>( "Uniform" )
  
  .constructor<double,double>()

  .field( "min", &Uniform::min )
  .field( "max", &Uniform::max )

  .method( "draw", &Uniform::draw )
  .method( "range", &uniformRange )
  ;
}

After sourcing the file the Uniform class can be used:

library(methods)    ## needed for S4 

## Create new instance myuniform:
myuniform <- new(Uniform, 0, 10)

# Print the new uniform isntance:
myuniform

C++ object <0x55fd63177450> of class 'Uniform' <0x55fd63788a20>

What happens now, is that the uniform instance calls its default print method – which results in a fairly uninformative display. It would be nice if the printer of the instances could be customized to provide more information about the specific object.

Customize the Printer

It is possible to make use of the underlying S4 structure of the exposed C++ classes:

# Check exposed S4 structure:
isS4(myuniform)

[1] TRUE

# Get class name:
class(myuniform)

[1] "Rcpp_Uniform"
attr(,"package")
[1] ".GlobalEnv"

For the Uniform class this is Rcpp_Uniform. To obtain a custom printer the last step now is to set the method show and define the function which should be used as printer:

# Define the printer:
ignoreMe <- setMethod("show", "Rcpp_Uniform", function (object) {
  cat("\n Hi, I am an uniform object!\n")
  cat("\n I was initialized with a minimum value of", object$min)
  cat("\n and a maximum value of ", object$max, ".\n", sep = "")
  cat("\n Therefore my range is ", object$range(), ".", sep = "")
  cat("\n\n")
})

# Test the printer:
myuniform


 Hi, I am an uniform object!

 I was initialized with a minimum value of 0
 and a maximum value of 10.

 Therefore my range is 10.

This works very nicely. Now it is possible to provide some more informations about a new Uniform instance. One thing to note is that setMethod returns the method name as string, which we assign to an otherwise unused variable.

Use in Packages

To get that print method as a default printer after exposing the C++ class to R within a package, it is sufficient to create a R file (e.g. R/uniform_printer.R) and put the following code in there:

ignoreMe <- setMethod("show", "Rcpp_Uniform", function (object) {
  cat("\n Hi, I am an uniform object!\n")
  cat("\n I was initialized with a minimum value of", object$min)
  cat("\n and a maximum value of ", object$max, ".\n", sep = "")
  cat("\n Therefore my range is ", object$range(), ".", sep = "")
  cat("\n\n")
})

Creating integer64 and nanotime vectors in C++

Sat, 11 Nov 2017 00:00:00 +0000

Motivation: More Precise Timestamps

R has excellent facilities for dealing with both dates and datetime objects. For datetime objects, the POSIXt time type can be mapped to POSIXct and its representation of fractional seconds since the January 1, 1970 “epoch” as well as to the broken-out list representation in POSIXlt. Many add-on packages use these facilities.

POSIXct uses a double type to provide 53 bits of resolution. That is generally good enough for timestamps down to just above a microsecond, and has served the R community rather well.

But increasingly, time increments are measured in nanoseconds. Other languages uses a (signed) 64-bit integer to represent (integer) nanoseconds since the epoch. A bit over a year ago I realized that we have this in R too—by combining the integer64 type in the very nice bit64 package by Jens Oehlschlaegel with the CCTZ-based parser and formatter in my own RcppCCTZ package. And thus the nanotime package was created.

Since then, Leonardo Silvestri joined in and significantly enhanced nanotime by redoing it as a S4 class.

A simple example:

library(nanotime)
n <- nanotime(42)
n

[1] "1970-01-01T00:00:00.000000042+00:00"

Here we used a single element with value 42, and created a nanotime vector from it—which is taken to me 42 nanoseconds since the epoch, or basically almost at January 1, 1970. See the nanotime page and package for more.

Step 1: Large Integer Types

So more recently I had a need to efficiently generate such (many such) integer vectors from int64_t data. Both Leonardo and Dan helped with initial discussions and tests. One can either use a reinterpret_cast<>, or a straight memcpy as the key trick in bit64 is to (re-)use the underlying int64_t representation (which we do not have in R) via the 64-bit double representative. Just never access it as a double. So we have the space, we just need to ensure we copy the bits (i.e. actual binary content) rather than their value (when “mapped” to a type). This leads to the following function to create an integer64 vector for use in R at the C++ level:

#include <Rcpp.h>
                                        
Rcpp::NumericVector makeInt64(std::vector<int64_t> v) {
    size_t len = v.size();
    Rcpp::NumericVector n(len);         // storage vehicle we return them in

    // transfers values 'keeping bits' but changing type
    // using reinterpret_cast would get us a warning
    std::memcpy(&(n[0]), &(v[0]), len * sizeof(double));

    n.attr("class") = "integer64";
    return n;
}

This uses the standard trick of setting a class attribute to set an S3 class. Now the values in v will return to R (exactly how is treated below), and R will treat the vector as integer64 object (provided the bit64 package has been loaded).

As mentioned, reinterpret_cast<>() can be used too, but leads to a compiler warning (under g++-6). Per Matt’s excellent compiler explorer, both approaches lead to the same mov semantics, so we prefer the variant that does not yell at us.

Step 2: Nanotime

A nanotime vector is creating using an internal integer64 vector. So the previous function almost gets us there. But we need to set the S4 type correctly. So that needed some extra work—and the following function seems to do it right:

#include <Rcpp.h>

Rcpp::S4 makeNanotime(std::vector<int64_t> v) {
    size_t len = v.size();
    Rcpp::NumericVector n(len);         // storage vehicle we return them in

    // transfers values 'keeping bits' but changing type
    // using reinterpret_cast would get us a warning
    std::memcpy(&(n[0]), &(v[0]), len * sizeof(double));

    // do what needs to be done for the S4-ness: class, and .S3Class
    // this was based on careful reading of .Internal(inspect(nanotime(c(0,1))))
    Rcpp::CharacterVector cl = Rcpp::CharacterVector::create("nanotime");
    cl.attr("package") = "nanotime";
    n.attr(".S3Class") = "integer64";
    n.attr("class") = cl;
    SET_S4_OBJECT(n);

    return Rcpp::S4(n);
}

This creates a nanotime vector as a proper S4 object. As before, we set some class attributes (though in a nested fashion as S4 is that fancy) and also invoke one R macro.

Step 3: Returning them R via data.table

The astute reader will have noticed that neither one of the functions presented so far had an Rcpp::export tag. This is because of their function argument: int64_t is not representable natively by R, which is why we need a workaround.

Matt Dowle has been very helpful in providing excellent support for nanotime in data.table (even after we, ahem, borked it by switching from S3 to S4). This support was of course relatively straightforward because data.table already had support for the underlying integer64, and we had the additional formatters etc.

#include <Rcpp.h>

// Enable C++11 via this plugin (Rcpp 0.10.3 or later)
// [[Rcpp::plugins("cpp11")]]

Rcpp::NumericVector makeInt64(std::vector<int64_t> v) {
    size_t len = v.size();
    Rcpp::NumericVector n(len);         // storage vehicle we return them in

    // transfers values 'keeping bits' but changing type
    // using reinterpret_cast would get us a warning
    std::memcpy(&(n[0]), &(v[0]), len * sizeof(double));

    n.attr("class") = "integer64";
    return n;
}

Rcpp::S4 makeNanotime(std::vector<int64_t> v) {
    size_t len = v.size();
    Rcpp::NumericVector n(len);         // storage vehicle we return them in

    // transfers values 'keeping bits' but changing type
    // using reinterpret_cast would get us a warning
    std::memcpy(&(n[0]), &(v[0]), len * sizeof(double));

    // do what needs to be done for the S4-ness: class, and .S3Class
    // this was based on careful reading of .Internal(inspect(nanotime(c(0,1))))
    Rcpp::CharacterVector cl = Rcpp::CharacterVector::create("nanotime");
    cl.attr("package") = "nanotime";
    n.attr(".S3Class") = "integer64";
    n.attr("class") = cl;
    SET_S4_OBJECT(n);

    return Rcpp::S4(n);
}

// [[Rcpp::export]]
Rcpp::DataFrame getDT() {
    std::vector<int64_t> d = { 1L, 1000L, 1000000L, 1000000000L };
    std::vector<int64_t> ns = { 1510442294123456789L, 1510442295123456789L, 
        1510442296123456789L, 1510442297123456789L };
        
    Rcpp::DataFrame df = Rcpp::DataFrame::create(Rcpp::Named("int64s") = makeInt64(d),
                                                 Rcpp::Named("nanos") = makeNanotime(ns));
    df.attr("class") = Rcpp::CharacterVector::create("data.table", "data.frame");
    return(df);
}

Example

The following example shows the output from the preceding function:

suppressMessages(library("data.table"))
dt <- getDT()
print(dt)

       int64s                               nanos
1:          1 2017-11-11T23:18:14.123456789+00:00
2:       1000 2017-11-11T23:18:15.123456789+00:00
3:    1000000 2017-11-11T23:18:16.123456789+00:00
4: 1000000000 2017-11-11T23:18:17.123456789+00:00

dt[[1]]

integer64
[1] 1          1000       1000000    1000000000

dt[[2]]

[1] "2017-11-11T23:18:14.123456789+00:00"
[2] "2017-11-11T23:18:15.123456789+00:00"
[3] "2017-11-11T23:18:16.123456789+00:00"
[4] "2017-11-11T23:18:17.123456789+00:00"

diff(dt[[2]]) # here 1e9 nanoseconds between them

integer64
[1] 1000000000 1000000000 1000000000

With that we’re done for this piece. Happy nanotime-ing from C++!

Rcpp Gallery

Extending R with C++ and Fortran

Benchmarking Rcpp code with RcppClock

A simple example

Multiple replicates

Simulation Smoother using RcppArmadillo

Precision sampler

Three implementations in RcppArmadillo

Three implementations in R

Performance comparison

Conclusion

Constructing a Sparse Matrix Class in Rcpp

Introduction

Structure of a dgCMatrix

A dgCMatrix reference class for Rcpp

Using the Rcpp::dgCMatrix class

Passing by reference versus value

Sparse iterator class

Extending the class

Using iconv via the R exported header

Introduction

Using Scheme from R via Guile and Rcpp

Introduction

The Guile Scheme Script script.scm

The C wrapper testguile.c

The Rcpp Caller

The System Integration

Package

Guile Demo

Handling R6 objects in C++

Introduction

Creating a dummy R6 class

Creating R6 objects in C++

Calling the methods of an instance

Passing down R6 objects

Mixing Rcpp modules and Rcpp attributes

Introduction

The RCPP_EXPOSED* macros

Recreating rmultinom and rpois with Rcpp

Nullable Optional Arguments in Rcpp functions

Introduction

Nullability of Vector, Matrix or Logical Vector

Nullability of DataFrame and List

Nullability of RcppGSL::Matrix

Summary

Creating a data.table from C++

Introduction

Very Short Answer

Detailed Answer

Documenting Rcpp functions and classes in R packages

Introduction

Create a package

Add C++ functionality and Roxygen2 comment blocks

Approach 1: Standalone functions

Approach 2: Nested field structure

Approach 3: Stand-alone pages for individual class methods

Approach 4: Rcpp module docstrings

Compile, install and test

Going Further

Using RcppParallel to aggregate to a vector

The General Problem

The parallelReduce Worker

Demonstration

Performance considerations with sparse matrices in Armadillo

Introduction

General considerations: sparsity, row vs column access

Matrix multiplication

Row vs column access

The R–C++ interface

Iterators vs elementwise access

Suppressing Call Stack Info in Rcpp-Generated Errors and Warnings

Introduction

Error Messages

Warnings

A C++11 Implementation

Introducing RcppArrayFire

Introduction

Installation

A first example

Arrays as parameters

Structure of a `dgCMatrix`

The Guile Scheme Script `script.scm`

The C wrapper `testguile.c`

The `RCPP_EXPOSED*` macros

Nullability of `RcppGSL::Matrix`