translating R's %in% to C++

Question

I am trying to translate a piece of R code to C++ (via Rcpp) that uses the %in% operator. I know Rcpp has a sugar in(), but I don't understand how it relates. There is another SO question regarding %in% but I couldn't sort out how that one relates too. The Rcpp gallery is wonderful but when I search in() I don't find any posts.

Here is a simple example of the problem. Suppose you want to create an identity matrix but constrain some elements on the diagonal to be zero. (This is obviously no longer an identity matrix but you understand what I mean.) You can do it in R with a function like this:

r_eye <- function(n, v) {
  A <- matrix(data = 0, nrow = n, ncol = n)
  for(i in 1:n) {
    for(j in 1:n) {
      if(i==j & !(i %in% v)) {
        A[i,j] <- 1
      }
    }
  }
  return(A)
}

where n is an integer that determines the size of the square matrix and v is a integer vector of "constraints". Starting with a matrix of zeros, you loop over rows and columns. If row equals column, and if row is not among the constraints, assign the value one to the matrix entry. For instance, if you run r_eye(n=6, v=c(1,2,3)), you get an "identity" matrix where the first three diagonal elements are zero and the remaining are one:

> r_eye(n = 6, v = c(1,2,3))
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    0    0    0    0    0    0
[2,]    0    0    0    0    0    0
[3,]    0    0    0    0    0    0
[4,]    0    0    0    1    0    0
[5,]    0    0    0    0    1    0
[6,]    0    0    0    0    0    1

The bottleneck to translating this to C++ is the second part of the "if" condition: !(i %in% v). I can create an identity matrix with

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
NumericMatrix cpp_eye(int n) {
  NumericMatrix A(n,n);
  for(int i = 0; i <= n; i++) {
    for(int j = 0; j <= n; j++) {
      if(i==j) {
          A(i,j) = 1;
        } 
      }
    }
  return A;
}

but is it possible to squeeze in that second "if" condition like the R function does? Or do I have to consider a different approach?

@F.Privé the point of the R function was to make it look like the cpp function — invictus, Jul 26 '17 at 13:23

bjhend · Answer 1 · 2017-07-25T23:17:07.383

0

As far as I could find out, an R vector is passed as Rcpp::NumericVector, which provides begin() and end() like a std::vector. So I would try this:

#include <algorithm>

bool in(int value, const NumericVector & vec) {
    return vec.end() != std::find(vec.begin(), vec.end(), value);
}


NumericMatrix cpp_eye(int n, NumericVector v) {
    NumericMatrix A(n,n);
    for(int i = 0; i <= n; i++) {
        for(int j = 0; j <= n; j++) {
            if((i==j) && in(i, v)) {
               A(i,j) = 1;
            } 
        }
    }
    return A;
}

edited Jul 25 '17 at 23:17

answered Jul 25 '17 at 22:27

bjhend

1,538
11
25

can you please explain the custom `in()` function here? – invictus Jul 25 '17 at 22:37
It does the same as ``i %in% v`` does in R. It applies ``std::find`` from STL's algorithm lib to search for an element ``value`` by iterating over ``vec`` from ``vec.begin()`` until ``vec.end()``. ``std::find`` returns either an iterator pointing to the found element or ``vec.end()``. So, by comparing the result to ``vec.end()`` you know if ``value`` is contained in ``vec`` or not. Note, that ``vec.end()`` points _behind_ the last element of ``vec``. – bjhend Jul 25 '17 at 22:44
cool. do you know how fast `std::find()` is? From the related SO post I linked to they say it is slow. If so C++ may not speed things up. – invictus Jul 25 '17 at 22:50
``std::find`` goes through the elements successively. Thus, if ``n`` is very big it might pay off to sort ``v`` first and then use ``std::binary_search`` instead (it already returns a ``bool`` so comparing with ``vec.end()`` should be removed). Another alternative ist to copy ``v`` into a ``std::set`` and use its ``find`` function, which also does a binary search. – bjhend Jul 25 '17 at 23:06

score 0 · Answer 2 · answered Jul 26 '17 at 08:06

I don't understand why you want to use %in% in this situation. Even for a more complex problem, you should go all over v rather than go all over the indices and test if they are part of v.

For this example, you could do:

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
NumericMatrix cpp_eye(int n, const IntegerVector& v) {

  int i, j;

  NumericMatrix A(n, n);
  for (i = 0; i < n; i++) A(i, i) = 1;

  for (j = 0; j < v.size(); j++) {
    i = v[j] - 1; // we are in C++
    A(i, i) = 0;
  }

  return A;
}

"I don't understand..." because it didn't occur to me, mon frere — invictus, Jul 26 '17 at 13:24

translating R's %in% to C++

2 Answers2