Worked Examples
Motivation
The previous package skeleton left out some essential details, such as testing for memory leaks and debugging. This chapter will cover these aspects in more detail. The references for this chapter are Vaughan, Hester, and Francois (2024), Padgham (2022), Vaughan (2019), and Wickham and Bryan (2023).
Instructional examples
Solving a matrix using the Gauss-Jordan method
This example package covers the following topics:
- Integers
- Doubles
- Matrices
- Conditionals
- Loops (
for) - Vendoring
See the cpp11gaussjordan package.
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 cpp11 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 cpp11 package could never break the package’s code. The disadvantage is that fixes and new features will be available after vendoring the code again.
Vendoring
I started with create_package("~/github/cpp11gaussjordan"). I used VSCode but all my steps also apply to RStudio.
After opening ~/github/cpp11gaussjordan I ran use_cpp11() to have a readily availably skeleton in my project.
I ran use_apache_licence() to have a LICENSE file and indicate in DESCRIPTION that my package is distributed under the Apache License.
Then I ran cpp_vendor() to copy the C++ headers into inst/include.
Building and testing
I used devtools to build and test the package:
# build
devtools::clean_dll()
cpp11::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.4Naive Ordinary Least Squares (OLS) estimator
This example package covers the following topics:
- Integers
- Doubles
- Doubles matrices
- Conditionals
See the cpp11ols 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.
My approach was to create one function per step, which means to create 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 call 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().
Linear programming (Simplex phase 2)
This example package covers the following topics:
- Integers
- Doubles
- Doubles matrices
- Conditionals
- Loops (
forandwhile) - Messages
See the cpp11simplex package.
Algorithm
The simplex algorithm is well described in Introduction to Linear Optimization and there is efficient software to solve this, including Berkelaar and Csardi (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:
- 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.
- If \(c_j < 0\) for some \(j\), we choose it to enter the base. We chose the variable with the most negative \(c_j\), let’s say that it is \(j = s\).
- If \(a_{is} \leq 0\) for all \(i\), then the problem is unbounded.
- If \(a_{is} > 0\) for some \(i\), we choose \(i = r\) such that \(\frac{b_r}{a_{rs}} = \min(\frac{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:
- \(a_{ij} \leftarrow a_{ij} - \frac{a_{is} a_{rj}}{a_{rs}}\) for \(j \neq s\)
- \(a_{rj} \leftarrow \frac{a_{rj}}{a_{rs}}\)
- \(b_i \leftarrow b_i - \frac{a_{is} b_r}{a_{rs}}\) for \(i \neq r\)
- \(b_r \leftarrow \frac{b_r}{a_{rs}}\)
- \(c_j \leftarrow c_j - \frac{c_s a_{rj}}{a_{rs}}\)
- \(-z \leftarrow -z - \frac{c_s b_r}{a_{rs}}\)
This algorithm is equivalent to Gauss method to solve linear systems.
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 we reached a stopping criterion: the minimum cost is non-negative, therefore the solution is optimal and is \(x^* = (\frac{1}{4}, \frac{11}{4}, 0 , 0)\) with an optimal value \(z^* = -\frac{17}{2}\).
Building and testing
I used devtools to build and test the package:
## build
devtools::clean_dll()
cpp11::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
)
cpp11_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 !Using OMP (parallelization)
This example package covers the following topics:
- Integers
- Doubles
- Lists (unnamed and named)
- Loops (
for) - OpenMP parallelization
See the cpp11omp package.
Motivation
One common phrase that I find when I need to Google how to do something with cpp11 is: “cpp11 does not offer OpenMP support.” This is a myth. cpp11 does offer OpenMP support. The requirements are: A processor and C++ compiler that support OpenMP.
I tested this on Windows, where you need to install Rtools, and Ubuntu where I didn’t 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 (“Enable OpenMP on macOS For Building Scikit-Learn - OpenBox Documentation” n.d.).
Enabling OpenMP
I created an R package called cpp11omp with the following code:
usethis::create_project("cpp11omp")
usethis::use_cpp11()Then, I created the file R/cpp11omp-package.R with the following contents:
### usethis namespace: start
#' @useDynLib cpp11omp, .registration = TRUE
### usethis namespace: end
NULLIn order to get the #pragma instruction to work, I needed to add the following to src/Makevars:
PKG_CXXFLAGS = $(SHLIB_OPENMP_CXXFLAGS) -DARMA_OPENMP_THREADS=1
PKG_LIBS = $(SHLIB_OPENMP_CXXFLAGS)
CXX_STD = CXX11Not adding this means that the pragma instruction will be ignored and the functions will run in a single thread.
Unnamed list
I added 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 <cpp11.hpp>
#include <omp.h>
using namespace cpp11;
[[cpp11::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 [[cpp11::register]] so that it can be called from R.
C++ is strict with types, so I need to create a wrapper function that will convert the integers to doubles to avoid accidental errors, it will go in R/cpp11omp-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))
}Named list
I added 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:
[[cpp11::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, I added 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))
}Building and testing
I used devtools to build and test the package:
cpp11::cpp_register()
devtools::document()
devtools::install()Then, I tested the package from a new R session:
> library(cpp11omp)
> 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 7Useful 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. - cpp11armadillo: An R package that provides bindings to the Armadillo C++ library. Armadiilo 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 permisive license.
- cpp11tesseract: An R package that provides bindings to the Tesseract OCR
C++ engine. This package allows to extract text from images. - 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
cpp11. - tidyr: A package that uses C++ functions to reshape data frames.