Re: [eigen] On a flexible API for submatrices, slicing, indexing, masking, etc.

[Perry de Valpine <pdevalpine@xxxxxxxxxxxx>]

Hi all,

I’m very interested in this topic because we make use of Eigen via a domain-specific language (DSL)  embedded in R (http://r-nimble.org), and the limitations especially of indexing and masking have been an obstacle.  I’ve tried to catch up on the extended discussion here, but apologies if I miss or repeat something.

More context: Specifically, NIMBLE generates C++ from code written in our nimbleFunction format within R, so we are compiling from a subset of R to C++.  For linear algebra, we generate code for Eigen, both the syntax for actual computations and also declaration of the Eigen objects and necessary initializations.  This means both that our use of Eigen is limited in scope and sometimes (often?) sub-optimal compared to hand-written code, but nevertheless it is incredibly useful to us.  Because we use our own objects to allocate memory and do other stuff, we always use Eigen::Maps.

Recently I have drafted code specific to NIMBLE’s needs that overlaps with some of the discussion here and also covers some different needs that we have.  I did most of it by defining new NullaryOps — thanks for those!  Some of the implemented features and their R counterparts that we aim to support include:

- concatenation of 2, 3 or 4 vectors (c() in R)

- construction of vectors by repeating a vector in sequence or element-wise (rep() in R)

- construction of sequences in a variety of ways (seq() and `:` in R) — already under discussion on this thread

- conversion of a boolean vector to (1-based) indices of its TRUE elements (which() in R)

- indexing with arbitrary integer or boolean vectors (A[ booleanVector, integerVector] ) — already under discussion on this thread.  (This got awkward with NullaryOps when I wanted to assign to arbitrarily indexed matrix elements, so I had to implement that separately).

- interacting with, extracting, or constructing a matrix by, the diagonal elements of a matrix (diag() in R).

- vectorized use of distribution functions (or any function defined natively to take scalar arguments that one wants to vectorize) to imitate R’s “recycling rule”.  This is a system in R that recycles elements of vector arguments as much as needed to match the longest argument.  Example: dnorm(x, mean, sd) where x is a vector of 5 numbers, mean is a scalar and sd is a vector of 2 numbers will return the vector of 5 normal distribution density function evaluations for the 5 values of x, the mean recycled 5 times, and the sd recycled 2 and a half times (yielding 5 values for all arguments).

It would be fantastic if any pieces of the above functionality could become native to Eigen, although some of the needs are clearly more general than others. 

The way I set these up using NullaryOps is limited.  The NullaryOp provided a fantastic hook to make these features work semi-natively with other Eigen operations, but on the other hand I suspect some of the efficiencies of Eigen are lost in the way I’ve done it (and I simply haven’t dealt with some issues yet, like defining meaningful functor_traits<…>::Cost values or figuring out packets).  Also the NullaryOp provided a feasible degree of learning curve whereas diving deeper into Eigen’s design would have been much harder.

In case it is of interest, the code (still in a development branch of NIMBLE) is mostly here: https://github.com/nimble-dev/nimble/blob/devel/packages/nimble/inst/include/nimble/nimbleEigen.h.  It may not be clear how we use it since it is only used by NIMBLE-generated C++, and some of the C++ we wrote is designed to simplify the C++ we need to generate.  Example: To concatenate two vectors of doubles, we generate nimCd(expression1, expression2). This finds a templated implementation at compile time, which constructs a CwiseNullaryOp containing a functor that provides operator() for one or two arguments and evaluates elements of expression1 and expression2 when needed.

In case it is of further interest, we also use a wiki for various design issues, and the page on this topic is here: https://github.com/nimble-dev/nimble/wiki/Summary-of-core-R-features-implemented-Fall-2016

Perry