Rcpp function for adding elements of a vector

Question

I have a very long vector of parameters (approximately 4^10 elements) and a vector of indices. My aim is to add together all of the values of the parameters that are indexed in the indices vector.

For instance, if I had paras = [1,2,3,4,5,5,5] and indices = [3,3,1,6] then I would want to find the cumulative sum of the third value (3) twice, the first value (1) and the sixth (5), to get 12. There is additionally the option of warping the parameter values according to their location.

I am trying to speed up an R implementation, as I am calling it millions of times.

My current code always returns NA, and I can't see where it is going wrong

Here's the Rcpp function:

double dot_prod_c(NumericVector indices, NumericVector paras, 
                   NumericVector warp = NA_REAL) {
int len = indices.size();
LogicalVector indices_ok;
for (int i = 0; i < len; i++){
    indices_ok.push_back(R_IsNA(indices[i]));
}
if(is_true(any(indices_ok))){
    return NA_REAL;
}
double counter = 0;
if(NumericVector::is_na(warp[1])){
    for (int i = 0; i < len; i++){
        counter += paras[indices[i]];
    }
} else {
    for (int i = 0; i < len; i++){
        counter += paras[indices[i]] * warp[i];
    }
}
return counter;
}

And here is the working R version:

dot_prod <- function(indices, paras, warp = NA){
    if(is.na(warp[1])){
        return(sum(sapply(indices, function(ind) paras[ind + 1])))
    } else {
        return(sum(sapply(1:length(indices), function(i){
            ind <- indices[i]
            paras[ind + 1] * warp[i]
        })))
    }
}

Here is some code for testing, and benchmarking using the microbenchmark package:

# testing
library(Rcpp)
library(microbenchmark)

parameters <- list()
indices <- list()
indices_trad <- list()

set.seed(2)
for (i in 4:12){
    size <- 4^i
    window_size <- 100
    parameters[[i-3]] <- runif(size)
    indices[[i-3]] <- floor(runif(window_size)*size)
    temp <- rep(0, size)
    for (j in 1:window_size){
        temp[indices[[i-3]][j] + 1] <- temp[indices[[i-3]][j] + 1] + 1
    }
    indices_trad[[i-3]] <- temp
}

microbenchmark(
    x <- sapply(1:9, function(i) dot_prod(indices[[i]], parameters[[i]])),
    x_c <- sapply(1:9, function(i) dot_prod_c(indices[[i]], parameters[[i]])),
    x_base <- sapply(1:9, function(i) indices_trad[[i]] %*% parameters[[i]])
)
all.equal(x, x_base) # is true, does work
all.equal(x_c, x_base) # not true - C++ version returns only NAs

Answer 1

I was having a little trouble trying to interpret your overall goal through your code, so I'm just going to go with this explanation

For instance, if I had paras = [1,2,3,4,5,5,5] and indices = [3,3,1,6] then I would want to find the cumulative sum of the third value (3) twice, the first value (1) and the sixth (5), to get 12. There is additionally the option of warping the parameter values according to their location.

since it was most clear to me.

---

There are some issues with your C++ code. To start, instead of doing this - NumericVector warp = NA_REAL - use the Rcpp::Nullable<> template (shown below). This will solve a few problems:

1. It's more readable. If you're not familiar with the Nullable class, it's pretty much exactly what it sounds like - an object that may or may not be null.
2. You won't have to make any awkward initializations, such as NumericVector warp = NA_REAL. Frankly I was surprised that the compiler accepted this.
3. You won't have to worry about accidentally forgetting that C++ uses zero-based indexing, unlike R, as in this line: if(NumericVector::is_na(warp[1])){. That has undefined behavior written all over it.

---

Here's a revised version, going off of your quoted description of the problem above:

#include <Rcpp.h>

typedef Rcpp::Nullable<Rcpp::NumericVector> nullable_t;
// [[Rcpp::export]]
double DotProd(Rcpp::NumericVector indices, Rcpp::NumericVector params, nullable_t warp_ = R_NilValue) {
  R_xlen_t i = 0, n = indices.size();
  double result = 0.0;

  if (warp_.isNull()) {
    for ( ; i < n; i++) {
      result += params[indices[i]];
    }    
  } else {
    Rcpp::NumericVector warp(warp_);
    for ( ; i < n; i++) {
      result += params[indices[i]] * warp[i];
    }  
  }

  return result;
}

---

You had some elaborate code to generate sample data. I didn't take the time to go through this because it wasn't necessary, nor was the benchmarking. You stated yourself that the C++ version wasn't producing the correct results. Your first priority should be to get your code working on simple data. Then feed it some more complex data. Then benchmark. The revised version above works on simple data:

---

args <- list(
  indices = c(3, 3, 1, 6),
  params = c(1, 2, 3, 4, 5, 5, 5),
  warp = c(.25, .75, 1.25, 1.75)
)

all.equal(
  DotProd(args[[1]], args[[2]]), 
  dot_prod(args[[1]], args[[2]]))
#[1] TRUE

all.equal(
  DotProd(args[[1]], args[[2]], args[[3]]), 
  dot_prod(args[[1]], args[[2]], args[[3]]))
#[1] TRUE

---

It's also faster than the R version on this sample data. I have no reason to believe it wouldn't be for larger, more complex data either - there's nothing magical or particularly efficient about the *apply functions; they are just more idiomatic / readable R.

---

microbenchmark::microbenchmark(
  "Rcpp" = DotProd(args[[1]], args[[2]]), 
  "R" = dot_prod(args[[1]], args[[2]]))
#Unit: microseconds
#expr    min      lq     mean  median      uq    max neval
#Rcpp  2.463  2.8815  3.52907  3.3265  3.8445 18.823   100
#R    18.869 20.0285 21.60490 20.4400 21.0745 66.531   100
#
microbenchmark::microbenchmark(
  "Rcpp" = DotProd(args[[1]], args[[2]], args[[3]]), 
  "R" = dot_prod(args[[1]], args[[2]], args[[3]]))
#Unit: microseconds
#expr    min      lq     mean median      uq    max neval
#Rcpp  2.680  3.0430  3.84796  3.701  4.1360 12.304   100
#R    21.587 22.6855 23.79194 23.342 23.8565 68.473   100

---

I omitted the NA checks from the example above, but that too can be revised into something more idiomatic by using a little Rcpp sugar. Previously, you were doing this:

LogicalVector indices_ok;
for (int i = 0; i < len; i++){
  indices_ok.push_back(R_IsNA(indices[i]));
}
if(is_true(any(indices_ok))){
  return NA_REAL;
}

It's a little aggressive - you are testing a whole vector of values (with R_IsNA), and then applying is_true(any(indices_ok)) - when you could just break prematurely and return NA_REAL on the first instance of R_IsNA(indices[i]) resulting in true. Also, the use of push_back will slow down your function quite a bit - you would have been better off initializing indices_ok to the known size and filling it by index access in your loop. Nevertheless, here's one way to condense the operation:

if (Rcpp::na_omit(indices).size() != indices.size()) return NA_REAL; 

---

For completeness, here's a fully sugar-ized version which allows you to avoid loops entirely:

#include <Rcpp.h> 

typedef Rcpp::Nullable<Rcpp::NumericVector> nullable_t;
// [[Rcpp::export]]
double DotProd3(Rcpp::NumericVector indices, Rcpp::NumericVector params, nullable_t warp_ = R_NilValue) {
  if (Rcpp::na_omit(indices).size() != indices.size()) return NA_REAL; 

  if (warp_.isNull()) {
    Rcpp::NumericVector tmp = params[indices];
    return Rcpp::sum(tmp);    
  } else {
    Rcpp::NumericVector warp(warp_), tmp = params[indices];
    return Rcpp::sum(tmp * warp); 
  }
}

/*** R

all.equal(
  DotProd3(args[[1]], args[[2]]), 
  dot_prod(args[[1]], args[[2]]))
#[1] TRUE

all.equal(
  DotProd3(args[[1]], args[[2]], args[[3]]), 
  dot_prod(args[[1]], args[[2]], args[[3]]))
#[1] TRUE

*/