16 - Worked Examples

1 Motivation

The previous package skeleton left out some essential details, such as testing for memory leaks and debugging. This vignette will cover these aspects in more detail. The references for this vignette are Vaughan 2024, Padgham 2022, Vaughan 2019, and Wickham 2019.

2 Instructional examples

2.1 Solving a matrix using the Gauss-Jordan method

This example package covers the following topics:

• Integers
• Doubles
• Matrices
• Conditionals
• Loops (for)
• Vendoring

See the cpp4rgaussjordan package.

2.1.1 Details

This implementation is a naive approach, but it can be used, for example, to obtain the Ordinary Least Squares (OLS) estimator as shown in the next section.

Vendoring means that the dependency code, the cpp4r C++ headers, are copied the project’s source tree. This ensures the dependency code is fixed and stable until it is updated.

The advantage of vendoring is that changes to the cpp4r package could never break the package’s code. The disadvantage is that fixes and new features will be available after vendoring the code again.

2.1.2 Vendoring

You can start with create_package("~/github/cpp4rgaussjordan"). You can use VSCode, but all the steps also apply to RStudio.

After opening ~/github/cpp4rgaussjordan, you can run use_cpp4r() to have a readily available skeleton in your project.

You can run use_apache_licence() to have a LICENSE file and indicate in DESCRIPTION that your package is distributed under the Apache License.

Then you can run cpp_vendor() to copy the C++ headers into inst/include.

2.1.3 Building and testing

You can use devtools to build and test the package:

# build

devtools::clean_dll()
cpp4r::cpp_register()
devtools::document()
devtools::load_all()

# test

A <- matrix(c(2,1,3,-1), nrow = 2, ncol = 2)
invert_matrix(A)

> invert_matrix(A)
     [,1] [,2]
[1,]  0.2  0.6
[2,]  0.2 -0.4

2.2 Ordinary Least Squares (OLS) estimator

This example package covers the following topics:

• Integers
• Doubles
• Doubles matrices
• Conditionals

See the cpp4rols package.

This implementation is extremely naive, quite similar to the Gauss-Jordan example with extra steps, and it is enough to show how to use C++ code within R.

One approach is to create one function per step, which means creating one function to obtain \(X^tX\), another for \((X^tX)^{-1}\) (e.g, implementing the Gauss-Jordan method to invert a matrix), another for \(X^tY\) and then calling each of those functions to obtain \(\hat{\beta} = (X^tX)^{-1}(X^tY)\).

A good challenge would be to implement the QR decomposition used by the lm() function in R and use it to obtain the OLS estimator in C++. Drury 2016 provides a good starting point, this is not trivial to implement.

It is hard to beat the performance of the lm() function in R, which calls compiled C and FORTRAN functions, and these functions are fast and robust lm().

2.3 Linear programming (Simplex phase 2)

This example package covers the following topics:

• Integers
• Doubles
• Doubles matrices
• Conditionals
• Loops (for and while)
• Messages

See the cpp4rsimplex package.

2.3.1 Algorithm

The simplex algorithm is well described in Introduction to Linear Optimization and there is efficient software to solve this, including BerkelaarCsardi 2023.

A problem written in canonical form is represented by a table such as:

$$ \begin{array}{ccc|c} x_1 & \cdots & x_n & \ \hline c_1 & \cdots & c_n & -z \ a_{11} & \cdots & a_{1n} & b_1 \ \vdots & \ddots & \vdots & \vdots \ a_{m1} & \cdots & a_{mn} & b_m \end{array} $$

where \(c_1, \ldots, c_n\) are the coefficients of the objective function (i.e., costs), \(a_{11}, \ldots, a_{mn}\) are the coefficients of the constraints, and \(b_1, \ldots, b_m\) are the right-hand side of the constraints.

The simplex algorithm to solve the problem consists in the next steps:

1. If \(c_j \geq 0\) for all \(j\), then the current solution is optimal. Basic variables are equal to \(b_i\) and non-basic variables are equal to 0.
2. If \(c_j < 0\) for some \(j\), choose it to enter the base. Choose the variable with the most negative \(c_j\), and assume that it is \(j = s\).
3. If \(a_{is} \leq 0\) for all \(i\), then the problem is unbounded.
4. If \(a_{is} > 0\) for some \(i\), choose \(i = r\) such that \(b_r / a_{rs} = \min(b_i / a_{is},\: a_{is} > 0)\) and pivot on \(a_{rs}\), to then go back to step 1.

The coefficients are updated according to:

1. \(a_{ij} \leftarrow a_{ij} - (a_{is} a_{rj}) / a_{rs}\) for \(j \neq s\)
2. \(a_{rj} \leftarrow a_{rj} / a_{rs}\)
3. \(b_i \leftarrow b_i - (a_{is} b_r) / a_{rs}\) for \(i \neq r\)
4. \(b_r \leftarrow b_r / a_{rs}\)
5. \(c_j \leftarrow c_j - (c_s a_{rj}) / a_{rs}\)
6. \(-z \leftarrow -z - (c_s b_r) / a_{rs}\)

This algorithm is equivalent to Gauss method to solve linear systems.

2.3.2 Numerical example

A simple example is the following minimization problem:

$$ \begin{aligned} \text{min} \quad & -x_1 - 3x_2 \ \text{subject to} \quad & x_1 + x_2 \geq 3 \ & -3x_1 + x_2 \geq 2 \ & x_1, x_2 \geq 0 \end{aligned} $$

In canonical form, this problem is:

$$ \begin{aligned} \text{min} \quad & -x_1 - 3x_2 + 0x_3 + 0x_4 \ \text{subject to} \quad & x_1 + x_2 + x_3 = 3 \ & -3x_1 + x_2 + x_4 = 2 \ & x_1, x_2,x_3,x_4 \geq 0 \end{aligned} $$

The initial tableau for the problem is:

$$ \begin{array}{cccc|c} x_1 & x_2 & x_3 & x_4 & -z \ \hline -1 & -3 & 0 & 0 & 0 \ 1 & 1 & 1 & 0 & 3 \ -3 & 1 & 0 & 1 & 2 \end{array} $$

The first row is the cost row, the last column is the right-hand side, and the rest is the matrix \(A\).

The first step is to pivot on row 2 and column 2:

$$ \begin{array}{cccc|c} x_1 & x_2 & x_3 & x_4 & -z \ \hline -10 & 0 & 0 & 3 & 6 \ 4 & 0 & 1 & -1 & 1 \ -3 & 1 & 0 & 1 & 2 \end{array} $$

The second step is to pivot on row 2 and column 1:

$$ \begin{array}{cccc|c} x_1 & x_2 & x_3 & x_4 & -z \ \hline 0 & 0 & 5/2 & 1/2 & 17/2 \ 1 & 0 & 1/4 & -1/4 & 1/4 \ 0 & 1 & 3/4 & 1/4 & 11/4 \end{array} $$

Here a stopping criterion is reached: the minimum cost is non-negative, therefore the solution is optimal and is \(x^* = (1/4, 11/4, 0 , 0)\) with an optimal value \(z^* = -17/2\).

2.3.3 Building and testing

You can use devtools to build and test the package:

# build

devtools::clean_dll()
cpp4r::cpp_register()
devtools::document()
devtools::load_all()

# test

c <- c(-1, -3)
b <- c(3, 2)

A <- matrix(
    c(1, -3, 1, 1),
    nrow = 2,
    ncol = 2,
    byrow = FALSE
)

cpp4r_simplex_phase2(c, b, A)

The result should be:

Initial tableau:
-1 -3  0  0  0 
 1  1  1  0  3 
-3  1  0  1  2 
Minimum cost: -3
Pivot row: 2
Pivot column: 2
====
New tableau:
-10  0  0  3  6 
 4  0  1 -1  1 
-3  1  0  1  2 
Minimum cost: -10
Pivot row: 1
Pivot column: 1
====
New tableau:
 0  0  2.5  0.5  8.5 
 1  0  0.25 -0.25  0.25 
 0  1  0.75  0.25  2.75 
Minimum cost: 0
Optimal solution found in 2 steps !

2.4 Using OMP (parallelization)

This example package covers the following topics:

• Integers
• Doubles
• Lists (unnamed and named)
• Loops (for)
• OpenMP parallelization

See the cpp4romp package.

2.4.1 Motivation

One common phrase that you might find when you need to search how to do something with cpp4r is: “cpp4r does not offer OpenMP support.” This is a myth. cpp4r does offer OpenMP support. The requirements are: A processor and C++ compiler that support OpenMP.

This has been tested on Windows, where you need to install Rtools, and Ubuntu where you do not need anything special because the gcc compiler comes with the operating system and it just works.

If you are using macOS, you need to install libomp via Homebrew in order to extend the clang compiler, and this is explained in the OpenBox documentation .

2.4.2 Enabling OpenMP

You can create an R package called cpp4romp with the following code:

usethis::create_project("cpp4romp")
usethis::use_cpp4r()

Then, you can create the file R/cpp4romp-package.R with the following contents:

## usethis namespace: start
#' @useDynLib cpp4romp, .registration = TRUE
## usethis namespace: end
NULL

In order to get the #pragma instruction to work, you need to add the following to src/Makevars:

PKG_CXXFLAGS = $(SHLIB_OPENMP_CXXFLAGS) -DARMA_OPENMP_THREADS=1
PKG_LIBS = $(SHLIB_OPENMP_CXXFLAGS)
CXX_STD = CXX11

Not adding this means that the pragma instruction will be ignored and the functions will run in a single thread.

2.4.3 Unnamed list

You can add a function called squared_unnamed_ in src/code.cpp that will square each element in a vector of doubles, so the file content corresponds to the following:

#include <cpp4r.hpp>
#include <omp.h>

using namespace cpp4r;

[[cpp4r::register]] list squared_unnamed_(doubles x) {
  // create vectors y = x^2 and z = thread number
  int n = x.size();
  writable::doubles y(n);
  writable::doubles z(n);
  
  #pragma omp parallel for
  for (int i = 0; i < n; ++i) {
    y[i] = x[i] * x[i];
    z[i] = omp_get_thread_num();
  }

  //create a list containing y and z
  writable::list out;
  out.push_back(y);
  out.push_back(z);
  return out;
}

The previous function returns an unnamed list with two elements: the squared vector and the thread number. The function is registered with [[cpp4r::register]] so that it can be called from R.

C++ is strict with types, so you need to create a wrapper function that will convert the integers to doubles to avoid accidental errors, it will go in R/cpp4romp-package.R:

#' Unnamed list with squared numbers and the threads used
#' @param x A vector of doubles
#' @export
squared_unnamed <- function(x) {
  squared_unnamed_(as.double(x))
}

2.4.4 Named list

You can add a function called squared_named_ in src/code.cpp that does the same but returns a named list. The additional content corresponds to the following:

[[cpp4r::register]] list squared_named_(doubles x) {
  // create vectors y = x^2 and z = thread number
  int n = x.size();
  writable::doubles y(n);
  writable::doubles z(n);
  
  #pragma omp parallel for
  for (int i = 0; i < n; ++i) {
    y[i] = x[i] * x[i];
    z[i] = omp_get_thread_num();
  }

  //create a list containing y and z
  writable::list out;
  out.push_back({"x^2"_nm = y});
  out.push_back({"thread"_nm = z});
  return out;
}

As in the previous part, you can add a wrapper and documentation:

#' Named list with squared numbers and the threads used
#' @param x A vector of doubles
#' @export
squared_named <- function(x) {
  squared_named_(as.double(x))
}

2.4.5 Building and testing

You can use devtools to build and test the package:

cpp4r::cpp_register()
devtools::document()
devtools::install()

Then, you can test the package from a new R session:

> library(cpp4romp)
> squared_unnamed(1:10)
[[1]]
 [1]   1   4   9  16  25  36  49  64  81 100

[[2]]
 [1] 0 0 1 1 2 3 4 5 6 7

> squared_named(1:10)
$`x^2`
 [1]   1   4   9  16  25  36  49  64  81 100

$thread
 [1] 0 0 1 1 2 3 4 5 6 7

2.5 Implementing a rejection and bootstrap sampler

This example package covers the following topics:

• Integers
• Doubles
• Matrices
• Push back
• Random seed state

3 Extensive R/C++ tests

See the cpp4rtest package. It covers over 1,000 unit tests for cpp4r features, including:

• Vectors (integers, doubles, logicals, strings)
• Matrices (integers, doubles, logicals, strings)
• Lists (named and unnamed)
• Data frames
• Error handling (stop, warning, message)
• C++ side unit tests with testthat
• Roxygen documentation

4 Useful R/C++ Integration Examples

Here are some examples of C++ code integration with R using the cpp11 package:

• arrow: An R package that provides bindings to the Arrow C++ library. Arrow is a
  columnar in-memory analytics format that is extremely fast and efficient.
• cpp4rarmadillo: An R package that provides bindings to the Armadillo C++ library. Armadillo is a high-quality linear algebra library with a syntax similar to MATLAB.
• cpp11eigen: An R package that provides bindings to the Eigen C++ library. Eigen is a high-performance linear algebra library with a permissive license.
• haven: A package that reads and writes SPSS, Stata, and SAS files in R.
• kendallknight: Implements the Kendall’s correlation coefficient in C++, achieving speedup by using an algorithm with a complexity of \(O(n \log(n))\) instead of \(O(n^2)\) in base R.
• mice: A package that imputes missing data using multivariate chained equations.
• redatam: C++ implementation of the Redatam file format, callable from both R and Python.
• RPostgres: A C++ interface to PostgreSQL using cpp4r.
• tidyr: A package that uses C++ functions to reshape data frames.

5 References