Overview

This post is a mind dump of working with the Armadillo library inside of R.

• Software Requirements for RcppArmadillo
• Using RcppArmadillo from a C++ script with sourceCpp()
• Examples of R syntax and conceptually corresponding Armadillo syntax
• Row and Column binds
• Subsetting operations

Software Requirements for RcppArmadillo

Unlike the majority of R packages that exist, the RcppArmadillo package requires a development toolchain to use. The toolchain varies from platform to platform. Please see one of the following guides for help:

• Windows: Installing RTools for Compiled Code via Rcpp
• macOS: R Compiler Tools for Rcpp on macOS
• Linux:
  • Debian (Ubuntu): Install the r-base-dev package via apt-get.
  • RHEL (CentOS): The development tools should be present.

Finally, please make sure to install the RcppArmadillo package via:

install.packages("RcppArmadillo")

Verify it works by running:

Rcpp::evalCpp("2 + 2", depends = "RcppArmadillo")
#> [1] 4

Common errors:

On macOS, if you receive an error such as -lgfortran and so on, then please make sure you installed the fortran binaries.

Using RcppArmadillo from a C++ script with sourceCpp()

The preferred method of compiling C++ code is to place it in a separate file and use sourceCpp("path/to/file.cpp") to compile and import C++ functions.

When writing C++ scripts that use armadillo, every C++ code file must have at the top:

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

1. The #include<RcppArmadillo.h> statement provides the Rcpp.h and armadillo.h headers with the appropriate casting magic.
2. The // [[Rcpp::depends(RcppArmadillo)]] statement adds onto the compiler’s search path the RcppArmadillo package directory where the armadillo.h header files are found. This approach takes advantage of Rcpp Attributes.

Optional After the Rcpp depends statement, you are able to use:

using namespace arma;

This avoids having code prefixed by a namespace. The downside to this approach is sometimes function names can overlap making for a fun evening of debugging.

Example

add_two.cpp

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

arma::vec add_two(arma::vec x){
  return x + 2;
}

add_two.R

# Load package
library("Rcpp")

# Compile C++ code and import function into R
sourceCpp("path/to/add_two.cpp")

# Example C++ function call from R
x = c(1,2,3)
add_two(1,2,3)

Examples of R syntax and conceptually corresponding Armadillo syntax

Based off of Armadillo’s syntax table with the left column containing R commands instead of Matlab/octave.

Row and Column Binds

When working with C++ objects, often there is a need to combine objects together. Traditionally, in R, there is an unlimited amount of times objects can be bound together. The same – unfortunately – cannot be said for C++. In particular, only two objects may be combined simultaneously. The objects can be joined together using either a row or a column matching scheme:

• rbind: join_vert(A,B) or join_cols(A,B)
• cbind: join_horiz(A,B) or join_rows(A,B)

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

// Performs the operation of rbind
// [[Rcpp::export]]
arma::mat bind_row_vectors(arma::vec A, arma::vec B){
  arma::mat out = join_vert(A, B); // same as join_cols
  return out;
}

// Performs the operation of cbind
// [[Rcpp::export]]
arma::mat bind_col_vectors(arma::vec A, arma::vec B){
  arma::mat out = join_horiz(A, B); // same as join_rows
  return out;
}

# Define two vectors
a = 1:3
b = 4:6

# Perform a row bind operation, e.g. rbind
bind_row_vectors(a, b)

##      [,1]
## [1,]    1
## [2,]    2
## [3,]    3
## [4,]    4
## [5,]    5
## [6,]    6

# Perform a column bind operation, e.g. cbind
bind_col_vectors(a, b)

##      [,1] [,2]
## [1,]    1    4
## [2,]    2    5
## [3,]    3    6

Subsetting with conditions

Accessing a portion of data or modifying values in place is possible by using one of the find*() variants:

• which(X > 2): find(X > 2)
• is.finite(X): find_finite(X)
  • Values excluding: -Inf, +Inf and NaN.
  • NB NA values in the integer case, e.g. ivec/imat, transform to INT_MIN. This method misses detecting this value. Use find(X == INT_MIND) to capture such values.
• is.infinite(x) & is.nan(x) or !is.finite(x): find_nonfinite(X)
  • Values that contain: -Inf, +Inf and NaN.
• seq_along(X)[!duplicated(X)]: find_unique(X)
  • Element index of the first instance of a value in a vector.

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

// Obtains the subset of values matching condition
// [[Rcpp::export]]
arma::vec find_equal(arma::vec A, double b){
  arma::uvec idx = find(A == b); // Substitute == with >, <, >=, <=, !=
  arma::vec out = A.elem(idx);   // Retrieve elements from positional index
  return out;
}

// [[Rcpp::export]]
arma::vec find_and_replace(arma::vec A,
                           double find_val = 3.0,
                           double replace_val = 1.5) {
  arma::uvec idx = find(A >= find_val); // Substitute == with >, <, >=, <=, !=
  A.elem(idx).fill(replace_val);        // Retrieve elements from positional index
                                        // Fill with the appropriate value
  return A;
}

# Define a vector
A = c(1,4,3,3,4,2,3,2)
b = 3

# Locate elements
find_equal(A, b)

##      [,1]
## [1,]    3
## [2,]    3
## [3,]    3

find_and_replace(A, b)

##      [,1]
## [1,]  1.0
## [2,]  1.5
## [3,]  1.5
## [4,]  1.5
## [5,]  1.5
## [6,]  2.0
## [7,]  1.5
## [8,]  2.0