In dirkschumacher/armacmp: Translate R Linear Algebra Code to Armadillo C++

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
library(armacmp)

This documents describes all language elements that can be compiled to Armadillo flavored C++ code.

Function signature

armacmp always compiles functions. Every function needs to have a return statement with an optional type argument.

my_fun <- compile(function(X, y = type_colvec())) {
  return(X %*% y, type = type_colvec())
}

You can define your inputs using the standard function syntax. By default paramters are of type type_matrix. But they can have other types such as:

• type_colvec - a column vector
• type_rowvec - a row vector
• type_scalar_integer - a single int value
• type_scalar_numeric - a single double value
• type_scalar_bool - a single boolean value

All functions need to return a value using the return function. The return function can have an optional return type argument with return(value, type = type_rowvec())

Type system

The package supports only a very limited number of types. For scalar values these are bool, int and double. All matrices and vectors are of type arma::mat/colvec/rowvec for double matrices, arma::imat/icolvec/irowvec for integer matrices and arma::umat/ucolvec/urowvec for unsigned integer matrices (used for indicies).

All operations are performed either on Armadillo types or on scalar types and not on R's data structures.

The package tries to figure out the correct types of variables and functions, but you can give type hints. If the package cannot figure out the type, the C++ type auto is used which can have suboptimal performance when used together with Armadillo operations.

Operators

Assignments

Example:

compile(function(X) {
  X2 <- X + 1
  X2 <- X2 + 1
  return(X2)
})

Multiplication

• %*% - matrix multiplication translates to the * operator on both arma::mat types and regular variables/expressions.
• * - elementwise multiplication gets translated to % for arma::mat types and to * otherwise.

The rest

The following operators are translated to their C++ counter-part.

• /
• +
• -
• !
• <
• >
• >=
• <=

Functions

General math functions

• x ^ y
• exp
• abs
• log
• sum
• min
• max
• sqrt
• cos
• acos
• cosh
• acosh
• sin
• asin
• sinh
• asinh
• tan
• atan
• tanh
• atanh

Linear algebra

• t
• chol
• diag
• QR
• svd
• solve
• ncol
• nrow
• crossprod
• tcrossprod
• det
• [i, j]
• [i, ]
• [, j]
• [i]
• length

Rest

• cumsum
• sort
• unique
• cumsum
• cumprod
• pnorm
• seq_len returns a arma::colvec (note that is corresponds to a sequence of numerics)
• a:b returns an arma::colvec from a to b (note that is corresponds to a sequence of numerics)
• rep.int returns a arma::colvec (note that is corresponds to a sequence of numerics)

Control Flow

Currently if/else, for and while loops are supported.

if/else

if_clause <- compile(function(X) {
  # infers that test needs to be bool (auto)
  test <- sum(log(X)) < 10
  X_new <- X
  if (test) {
    X_new <- (t(X) %*% X) + 10
  } else if (sum(X) < 10) {
    X_new <- t(X) %*% X
  } else {
    X_new <- (t(X) %*% X) + 10
  }
  return(X_new)
})

Just make sure all return statements return the same type.

Also if/else cannot be used as an expression:

# this is not possible
x <- if (y) 1 else 2

For loops

For loops are supported:

# currently only seq_len is available to create sequences of ints
x <- seq_len(n)
for (i in x) {
  ...
}

Example:

for_loop <- compile(function(X, offset = type_scalar_numeric()) {
  X_new <- X
  # only seq_len is currently supported
  for (i in seq_len(20)) {
    X_new <- log(t(X_new) %*% X_new + i + offset)
  }
  return(X_new)
})

While loops

x <- 0
while (x < 10) {
  x <- x + 1
}

Functions

You can also define functions within your functions (lambdas).

fun <- compile(function(X) {
  square <- function(y) {
    return(y^2)
  }
  Y <- square(X)
  return(Y)
})
fun(matrix(1:10))

Using lambdas it is possible to define recursive functions as well.

is_even <- compile(function(x = type_scalar_numeric()) {
  is_even_internal <- function(y = type_scalar_numeric()) {
    if (y == 0) return(TRUE)
    if (y < 0) return(FALSE)
    return(is_even_internal(y - 2), type = type_scalar_logical())
  }
  return(is_even_internal(x))
})
is_even(5)
is_even(4)

dirkschumacher/armacmp documentation built on Oct. 22, 2021, 7:10 p.m.