Cpp11 (R package) vendoring
R and Shiny Training: If you find this blog to be interesting, please note that I offer personalized and group-based training sessions that may be reserved through Buy me a Coffee. Additionally, I provide training services in the Spanish language and am available to discuss means by which I may contribute to your Shiny project.
Updated 2023-05-23: I explicitly mention to add the vendoring after creating the project instead of assuming the readers know I run a commented line in the install script. Also, if we include “cpp11.hpp”, then we don’t need to call “doubles.hpp” and “matrix.hpp” because it loads that and also “list.hpp”, “external_pointer.hpp” and other headers. Therefore, I simplified and kept just the call to doubles and matrix in order to load the minimal things required to build the package.
Motivation
In my previous post A step by step guide to write an R package that uses C++ code (Ubuntu) I explained the steps in my journey to use C++ function from R.
Now I will explain how I figured out how to use vendoring, which means that I copy the code for the dependencies into your project’s source tree, and in this case this means I copied the C++ headers provided by the cpp11
package into my own R package. 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 your existing code. The disadvantage is that you no longer get bugfixes and new features until you manually run cpp_vendor()
in your project.
Creating an R package with C++ code and vendoring
Similar to the previous post, I started with create_package("~/github/cpp11gaussjordan")
. I am using VSCode but all my steps also apply to RStudio.
After opening ~/github/cpp11gaussjordan
I run use_cpp11()
to have a readily availably skeleton in my project.
I run use_apache_licence()
to have a LICENSE
file and indicate in DESCRIPTION
that my package is distributed under the Apache License.
Then I run cpp_vendor()
to copy the C++ headers into inst/include
.
You can find the final result on my GitHub profile.
R side
I run use_r("cpp11gaussjordan-package")
and added the following contents.
#' @useDynLib cpp11gaussjordan, .registration = TRUE
NULL
#' Invert (some) square matrices
#' @export
#' @param A numeric matrix
#' @return numeric matrix
#' @examples
#' A <- matrix(c(2,1,3,-1), nrow = 2, ncol = 2)
#' invert_matrix(A)
<- function(A) {
invert_matrix invert_matrix_(A)
}
#' Solve (some) linear systems
#' @export
#' @param A numeric matrix
#' @param b numeric matrix
#' @return numeric matrix
#' @examples
#' A <- matrix(c(2,1,3,-1), nrow = 2, ncol = 2)
#' b <- matrix(c(7,4), nrow = 2, ncol = 1)
<- function(A,b) {
solve_system solve_system_(A, b)
}
Here solve_system()
is the “end-user” function, and solve_system_()
is the “development” function. The “development” functions were written in C++ and cpp11
handles adding the .Call()
when required to link the C++ code and be able to use the C++ functions from R.
C++ side
I replaced code.cpp
contents with the following lines.
#include <cpp11.hpp>
#include <cpp11/doubles.hpp>
using namespace cpp11;
[[cpp11::register]] doubles_matrix<> invert_matrix_(doubles_matrix<> A)
{
// Obtain (X'X)^-1 via Gauss-Jordan
// Check dimensions
int N = A.nrow();
int M = A.ncol();
if (N != M)
{
("X must be a square matrix");
stop}
// Copy the matrix
::doubles_matrix<> Acopy(N, N);
writable
for (int i = 0; i < N; i++)
{
for (int j = 0; j < N; j++)
{
(i, j) = A(i, j);
Acopy}
}
// Create the identity matrix as a starting point for Gauss-Jordan
::doubles_matrix<> Ainv(N, N);
writable
for (int i = 0; i < N; i++)
{
for (int j = 0; j < N; j++)
{
if (i == j)
{
(i, j) = 1.0;
Ainv}
else
{
(i, j) = 0.0;
Ainv}
}
}
// Overwrite Ainv by steps with the inverse of A
// in other words, find the echelon form of A
for (int i = 0; i < M; i++)
{
double a = Acopy(i, i);
for (int j = 0; j < M; j++)
{
(i, j) /= a;
Acopy(i, j) /= a;
Ainv}
for (int j = 0; j < M; j++)
{
if (i != j)
{
= Acopy(j, i);
a
for (int k = 0; k < M; k++)
{
(j, k) -= Acopy(i, k) * a;
Acopy(j, k) -= Ainv(i, k) * a;
Ainv}
}
}
}
return Ainv;
}
[[cpp11::register]] doubles_matrix<> solve_system_(doubles_matrix<> A, doubles_matrix<> b)
{
// Use Gauss-Jordan to solve the system Ax = b
// Check dimensions
int N1 = A.nrow();
int M1 = A.ncol();
if (N1 != M1)
{
("A must be a square matrix");
stop}
int N2 = b.nrow();
int M2 = b.ncol();
if (N1 != N2)
{
("b must have the same number of rows as A");
stop}
if (M2 != 1)
{
("b must be a column vector");
stop}
// Obtain the inverse
// writable::doubles_matrix<> Acopy(N1,N1);
// for (int i = 0; i < N1; i++)
// {
// for (int j = 0; j < N1; j++)
// {
// Acopy(i, j) = A(i, j);
// }
// }
// doubles_matrix<> Ainv = invert_matrix_(Acopy);
// I don´t need to create a copy in this case
<> Ainv = invert_matrix_(A);
doubles_matrix
// Multiply Ainv by b
::doubles_matrix<> x(N1, 1);
writable
for (int i = 0; i < N1; i++)
{
(i, 0) = 0.0;
x
for (int j = 0; j < N1; j++)
{
(i, 0) += Ainv(i, j) * b(j, 0);
x}
}
return x;
}
I also needed a Makevars
file to indicate the compiler (either clang
or gcc
) to use the vendored headers. The content is the following.
PKG_CPPFLAGS = -I../inst/include
Because of the vendoring option, I also had to remove the LinkingTo: cpp11
line from DESCRIPTION
.
Putting all together
I created a vscode-install.r
file in the root of the project, and I added it to .Rbuildignore
. It contains the following lines.
# cpp_vendor() # run only when updating C++ headers
clean_dll()
cpp_register()
document()
install()
# load_all()
To test that the functions work, after running install()
I solved a simple system.
> library(cpp11gaussjordan)
> A <- matrix(c(2,1,3,-1), nrow = 2, ncol = 2)
> b <- matrix(c(7,4), nrow = 2, ncol = 1)
> solve_system(A,b)
1]
[,1,] 3.8
[2,] -0.2
[> solve(A) %*% b
1]
[,1,] 3.8
[2,] -0.2 [
References
Ask a lot
Blog posts like this are summaries of what worked for me. I asked so much on Stackoverflow, that they created the “r-lib-cpp11” for my questions. Please ask there, there is probably something I figured how to solve after many hours on Google.