A b solve_system(A,b) [,1] [1,] 3.8 [2,] -0.2 > solve(A) %*% b [,1] [1,] 3.8 [2,] -0.2 References Cpp11 official documentation R-lib-cpp11 tag on Stackoverflow 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. " />

Cpp11 (R package) vendoring

[This article was first published on Pachá, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

R/Shiny Training: Should you find this blog to be of interest, kindly 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.

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)
invert_matrix <- function(A) {
  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)
solve_system <- function(A,b) {
  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)
    {
        stop("X must be a square matrix");
    }

    // Copy the matrix

    writable::doubles_matrix<> Acopy(N, N);

    for (int i = 0; i < N; i++)
    {        
        for (int j = 0; j < N; j++)
        {
            Acopy(i, j) = A(i, j);
        }
    }

    // Create the identity matrix as a starting point for Gauss-Jordan

    writable::doubles_matrix<> Ainv(N, N);

    for (int i = 0; i < N; i++)
    {
        for (int j = 0; j < N; j++)
        {
            if (i == j)
            {
                Ainv(i, j) = 1.0;
            }
            else
            {
                Ainv(i, j) = 0.0;
            }
        }
    }

    // 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++)
        {
            Acopy(i, j) /= a;
            Ainv(i, j) /= a;
        }

        for (int j = 0; j < M; j++)
        {
            if (i != j)
            {
                a = Acopy(j, i);

                for (int k = 0; k < M; k++)
                {
                    Acopy(j, k) -= Acopy(i, k) * a;
                    Ainv(j, k) -= Ainv(i, k) * a;
                }
            }
        }
    }

    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)
    {
        stop("A must be a square matrix");
    }

    int N2 = b.nrow();
    int M2 = b.ncol();

    if (N1 != N2)
    {
        stop("b must have the same number of rows as A");
    }

    if (M2 != 1)
    {
        stop("b must be a column vector");
    }

    // 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
    
    doubles_matrix<> Ainv = invert_matrix_(A);

    // Multiply Ainv by b

    writable::doubles_matrix<> x(N1, 1);

    for (int i = 0; i < N1; i++)
    {
        x(i, 0) = 0.0;

        for (int j = 0; j < N1; j++)
        {
            x(i, 0) += Ainv(i, j) * b(j, 0);
        }
    }

    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

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.

To leave a comment for the author, please follow the link and comment on their blog: Pachá.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)