#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
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 viaapt-get
. - RHEL (CentOS): The development tools should be present.
- Debian (Ubuntu): Install the
Finally, please make sure to install the RcppArmadillo package via:
install.packages("RcppArmadillo")
Verify it works by running:
::evalCpp("2 + 2", depends = "RcppArmadillo")
Rcpp#> [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:
- The
#include<RcppArmadillo.h>
statement provides theRcpp.h
andarmadillo.h
headers with the appropriate casting magic. - The
// [[Rcpp::depends(RcppArmadillo)]]
statement adds onto the compiler’s search path theRcppArmadillo
package directory where thearmadillo.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)]]
::vec add_two(arma::vec x){
armareturn 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
= c(1,2,3)
x 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.
R | Armadillo | Notes |
---|---|---|
A[1, 1] | A(0, 0) | indexing in Armadillo starts at 0 |
A[n, p] | A(n-1, p-1) | |
nrow(A) | A.n_rows | read only |
ncol(A,2) | A.n_cols | |
dim(Q)[3] | Q.n_slices | Q is a cube (3D array) |
length(A) | A.n_elem | |
A[, k] | A.col(k) | this is a conceptual example only; exact conversion from R to Armadillo syntax will require taking into account that indexing starts at 0 |
A[k, ] | A.row(k) | |
A[, p:q] | A.cols(p, q) | |
A[p:q,] | A.rows(p, q) | |
Q[,,k] | Q.slice(k) | Q is a cube (3D array) |
Q[,, t:u] | Q.slices(t, u) | |
t(A) | A.t() or trans(A) | matrix transpose / Hermitian transpose (for complex matrices, the conjugate of each element is taken) |
A = numeric(length(A)) | A.zeros() | |
A = rep(1, length(A)) | A.ones() | |
A = matrix(0, k, k) | A = zeros |
|
matrix(1, k, k) | A = ones |
|
A * B | A % B | element-wise multiplication |
A / B | A / B | element-wise division |
solve(A, B) | solve(A,B) | conceptually similar to inv(A)*B, but more efficient |
A = A + 1 | A++ | |
A = A - 1 | A– | |
matrix(c(1,2,3,4), nrow = 2) | A << 1 << 2 << endr << 3 << 4 << endr; | element initialisation, with special element endr indicating end of row |
X = c(A) | X = vectorise(A) | |
X = cbind(A, B) | X = join_horiz(A,B) | |
X = rbind(A, B) | X = join_vert(A,B) | |
print(A) | Rcpp::Rcout << A << std::endl; |
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]]
::mat bind_row_vectors(arma::vec A, arma::vec B){
arma::mat out = join_vert(A, B); // same as join_cols
armareturn out;
}
// Performs the operation of cbind
// [[Rcpp::export]]
::mat bind_col_vectors(arma::vec A, arma::vec B){
arma::mat out = join_horiz(A, B); // same as join_rows
armareturn out;
}
# Define two vectors
= 1:3
a = 4:6
b
# 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
andNaN
. - NB
NA
values in the integer case, e.g.ivec
/imat
, transform toINT_MIN
. This method misses detecting this value. Usefind(X == INT_MIND)
to capture such values.
- Values excluding:
is.infinite(x) & is.nan(x)
or!is.finite(x)
:find_nonfinite(X)
- Values that contain:
-Inf
,+Inf
andNaN
.
- Values that contain:
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]]
::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
armareturn out;
}
// [[Rcpp::export]]
::vec find_and_replace(arma::vec A,
armadouble find_val = 3.0,
double replace_val = 1.5) {
::uvec idx = find(A >= find_val); // Substitute == with >, <, >=, <=, !=
arma.elem(idx).fill(replace_val); // Retrieve elements from positional index
A// Fill with the appropriate value
return A;
}
# Define a vector
= c(1,4,3,3,4,2,3,2)
A = 3
b
# 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