How can I use RcppArmadillo?

[sangwon-hyun]

I see from templates that when I'm looking to use Rcpp, I can use this at the beginning of an Rcpp code block:

#include <Rcpp.h>
using namespace Rcpp;

What should I do when I'm looking to use libraries like Armadillo or Eigen (or numeric)? My guess would be to do:

#include <numeric>
#include <RcppArmadillo.h>
#include <RcppEigen.h>
#include <Rcpp.h>

using namespace Rcpp;
using namespace arma;
using namespace Eigen;

But that doesn't work :)

[jacobbien]

Thanks for this comment. c8e0f5e now makes it possible to use RcppArmadillo and adds a template. To get started, you can work from template:

rmarkdown::draft("create-witharmadillo.Rmd",
                 template = "make-an-r-package-with-armadillo",
                 package = "litr")
litr::render("create-witharmadillo.Rmd")

(First, install the latest development version of litr.)

[sangwon-hyun]

Ok thank you!

I'm happy to report that I'm able to use not only RcppArmadillo but also RcppEigen! Like this:

// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::depends(RcppEigen)]]
#include <RcppArmadillo.h>
#include <RcppEigen.h>
#include <numeric>
using namespace arma;
using namespace Eigen;

using Eigen::Map;                       // 'maps' rather than copies
using Eigen::MatrixXd;                  // variable size matrix, double precision

//' Solve "barebones" sylvester equation that takes upper triangular matrices as coefficients.
//' 
//' @param TA Upper-triangular matrix
//' @param TB Upper-triangular matrix
//' @param C matrix
//' @export
// [[Rcpp::export]]
Eigen::MatrixXd matrix_function_solve_triangular_sylvester_barebonesC2(const Eigen::MatrixXd & TA,
								     const Eigen::MatrixXd & TB,
								     const Eigen::MatrixXd & C){
  // Eigen::eigen_assert(TA.rows() == TA.cols());
  // Eigen::eigen_assert(TA.Eigen::isUpperTriangular());
  // Eigen::eigen_assert(TB.rows() == TB.cols());
  // Eigen::eigen_assert(TB.Eigen::isUpperTriangular());
  // Eigen::eigen_assert(C.rows() == TA.rows());
  // Eigen::eigen_assert(C.cols() == TB.rows());

  // typedef typename MatrixType::Index Index;
  // typedef typename MatrixType::Scalar Scalar;

  int m = TA.rows();
  int n = TB.rows();
  Eigen::MatrixXd X(m, n);

  for (int i = m - 1; i >= 0; --i) {
    for (int j = 0; j < n; ++j) {

      // Compute T_A X = \sum_{k=i+1}^m T_A_{ik} X_{kj}
      double TAX;
      if (i == m - 1) {
      	TAX = 0;
      } else {
	MatrixXd TAXmatrix = TA.row(i).tail(m-1-i) * X.col(j).tail(m-1-i);
      	TAX = TAXmatrix(0,0);
      }

      // Compute X T_B = \sum_{k=1}^{j-1} X_{ik} T_B_{kj}
      double XTB;
      if (j == 0) {
      	XTB = 0;
      } else {
      	MatrixXd XTBmatrix = X.row(i).head(j) * TB.col(j).head(j);
      	XTB = XTBmatrix(0,0);
      }

      X(i,j) = (C(i,j) - TAX - XTB) / (TA(i,i) + TB(j,j));
    }
  }
  return X;
}

works well! Here's some test code:

library(tidyverse)
set.seed(100)
A = runif(4) %>% matrix(ncol=2)
B = runif(4) %>% matrix(ncol=2)
C = runif(4) %>% matrix(ncol=2)
A[lower.tri(A)]=0
B[lower.tri(B)]=0
C[lower.tri(C)]=0
a <- matrix_function_solve_triangular_sylvester_barebonesC2(A, B, C)
a

I have another simple request! It seems like litr::document() is meant to be run at the end of any workflow. However, the current website showcasing the bookdown template doesn't make this clear:

and just ends with "etc." instead of showing 5conclude.Rmd.

(To elaborate, my workflow was to (1) run litr::render("index.Rmd") to create a package in ./mypackage, then (2) run install.package("mypackage", repos=NULL, type="source"). All without litr::document(). As I started to define Rcpp functions, I'd get the following error during this second step:

which, it turns out, litr::document() already anticipates! But the above error message by install.packages() is not useful for your average litr+bookdown user..)

Here's one simple suggestion: when showcasing templates, instead of ending with "etc.", maybe end with 5conclude.Rmd? Something like this:

index.Rmd, 1description.Rmd, ⋯, 5conclude.Rmd

Or an alternative is to showcase somewhere that the workflow is incomplete if you don't end with litr::document(). Maybe in the FAQ, or at the end of basic-example.html ?

[jacobbien]

Thanks, a03a34c follows both of these last two suggestions.