Rcpp Function - SegFault Error & Memory management

Question

I need to loop through columns in a matrix and sum all values in each row where column index is higher.

I've done it fine using a for loop and rowSums without issue as I'm familiar with basic R code. I wrote a C++ function to speed up runtime - code below. When I run in RStudio I encounter a fatal error. When I run in R I get a segfault error as memory is not mapped.

Error doesn't occur when I run matrix of say 10 rows - but ideally want to run 10k.

Do I need to allocate memory somewhere?

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
Rcpp::NumericVector PVCalc_cpp(Rcpp::NumericMatrix x, int y) {  
int nrow = x.nrow(), ncol = x.ncol();
RNGScope scope;

Rcpp::NumericVector out(nrow);
double total = 0;

      for (int i = 0; i < nrow; i++) {
        total = 0;
        for (int j = y; j < ncol; j++){
          total += x(i, j+1);
        }
        out(i) =  floor((100. * total)+.5)/100;
      }
  return out;
}

Consider what happens when `j = (ncol - 1)`. You will remove the error by changing `x(i, j+1)` to `x(i, j)`, but then I think you won't get the calculation you want. Can you show your working R code and/or a better mathematical description of your quantity of interest and/or example input and desired output? So that we can see exactly what you're trying to do; I didn't quite understand your description in the question. Then it will be easier to fix your out of bounds error. — duckmayr, Jun 02 '20 at 18:19
Also, an aside: Why do you declare `RNGScope scope;`? You don't use it — duckmayr, Jun 02 '20 at 18:22
So I want df_x[r,10] to equal the sum(df_y[r,11-1000]) essentially each column represents a time period - and I want the sum over all future time periods. ``` for (j in 1:(int_ncols-2)){ df_y[,j] <- rowSums(df_x[,c((j+1):int_ncols)]) } df_x[,(int_ncols-1)] <- df_x[,int_ncols]``` — Barry King, Jun 02 '20 at 18:59

duckmayr · Answer 1 · 2020-06-03T22:21:37.573

The problem is that when j = ncol - 1, trying to access x(i, j+1) in the line

total += x(i, j+1);

is the same as trying to access x(i, ncol), and you are out of bounds. If I understand your problem correctly, you can just change j+1 to j, assuming y is passed correctly. So, we could change your code to the following:

#include <Rcpp.h>

// [[Rcpp::export]]
Rcpp::NumericVector PVCalc_cpp(const Rcpp::NumericMatrix& x, int y) {  
    int nrow = x.nrow(), ncol = x.ncol();
    Rcpp::NumericVector out(nrow);
    double total = 0;
    for (int i = 0; i < nrow; i++) {
        total = 0;
        for (int j = y; j < ncol; j++){
            total += x(i, j);
        }
        /* The following will always be equal to total but rounded;
           to type less R code for the R comparison,
           I omit the rounding here.
           Notice also the use of [] access operator instead of ();
           this is more efficient.
         */
        // out[i] = floor((100. * total)+.5)/100;
        out[i] = total;
    }
    return out;
}

We can then verify that this both works, and is faster than rowSums():

## Load required packages
library(Rcpp)
library(microbenchmark)
## Source C++ code
sourceCpp("so.cpp")
## Generate example data
x <- matrix(1:9, nrow = 3)
x

     [,1] [,2] [,3]
[1,]    1    4    7
[2,]    2    5    8
[3,]    3    6    9

## Check results are the same
rowSums(x[ , -1])

[1] 11 13 15

PVCalc_cpp(x, 1)

[1] 11 13 15

## Benchmark small example
library(microbenchmark)
microbenchmark(base = rowSums(x[ , -1]),
               rcpp = PVCalc_cpp(x, 1))

Unit: microseconds
 expr   min     lq     mean median    uq      max neval
 base 5.591 5.9210  8.61073  6.475 6.786  137.125   100
 rcpp 2.337 2.5795 19.90118  3.035 3.222 1651.094   100

## Check larger example
set.seed(123)
x <- matrix(rnorm(1e6), nrow = 1e3)
y <- sample(seq_len(ncol(x)), size = 1)
all.equal(rowSums(x[ , -(1:y)]), PVCalc_cpp(x, y))

[1] TRUE

microbenchmark(base = rowSums(x[ , -(1:y)]),
               rcpp = PVCalc_cpp(x, y))

Unit: milliseconds
 expr      min       lq     mean   median       uq       max neval
 base 5.377342 6.052347 6.954338 6.482907 7.834190 11.580706   100
 rcpp 1.447596 1.909504 2.085185 2.023343 2.158256  3.159366   100

Thanks - this sped things up a lot for my smaller input matrix of about 500 * 720 records. However, when I jumped up to 1000 records I'm still getting the same error: *** caught segfault *** address 0x11583a000, cause 'memory not mapped' Traceback: 1: PVCalc_cpp(as.matrix(cf_pv_exp_temp[, c(int_matchcol:int_ncols)])) 2: data.frame(PVCalc_cpp(as.matrix(cf_pv_exp_temp[, c(int_matchcol:int_ncols)]))) 3: eval(ei, envir) 4: eval(ei, envir) 5: withVisible(eval(ei, envir)) 6: source("/.../CF Projection Code V2.R") — Barry King, Jun 02 '20 at 20:53
@BarryKing Such error should only occur if you pass `y` equal to `ncol(input_matrix)` — duckmayr, Jun 02 '20 at 20:55
I'm not sure I understand why it would work for a smaller matrix - but not for a similar matrix with a larger number of rows if it wasn't a memory issue? — Barry King, Jun 02 '20 at 21:22
@BarryKing See [Accessing an array out of bounds gives no error, why?](https://stackoverflow.com/q/1239938/8386140); accessing an element out of bounds is undefined behavior that may *or may not* cause a segfault — duckmayr, Jun 02 '20 at 21:55
Nice job as always. Minor suggestion: would it be better if you used the square brackets as opposed to the parentheses when assigning to the `out` vector? I haven't tested it, but I think the square bracket should be more efficient as it doesn't check bounds. — Joseph Wood, Jun 03 '20 at 18:53
@JosephWood That's correct. I typically do that, but didn't spend a lot of time optimizing here. Will put in an edit later — duckmayr, Jun 03 '20 at 19:04

Rcpp Function - SegFault Error & Memory management

1 Answers1