Rcpp function crashes

Question

My problem:

I am using R.3.0.1 together with RStudio 0.97.551 on a 64bit Windows7 PC and I have begun to outsource a function to C/C++ using Rcpp. The function compiles, but evaluating it within an R function produces a runtime error. I am not able to find out why and how to fix this.

Details

Below is my cpp-file... let's say it's called "vector.cpp"

#include <Rcpp.h>

using namespace Rcpp;

// [[Rcpp::export]]
NumericVector l5(double z, double k, double s, double K, double theta, double x, double h, NumericVector m){
    int n = m.size();
    NumericVector a(n);
    NumericVector bu(n);
    NumericVector b(n);
    NumericVector c(n);
    for(int i=0; i<n+1; i++){
        a[i] = pow(z,m[i]) * (pow((x*pow(h,m[i])/K), theta) - 1) * (K/theta);
        for (int j=0; j<i; j++){
            bu[i] += pow(z,j) * (1 - z) * fmax(((pow((s/K), theta) - 1) * (K/theta)), ((pow((x*pow(h, j)/K), theta) - 1) * (K/theta)));
        }
        b[i] = k *bu[i];
        c[i] = k * pow(z, m[i]) * fmax(((pow((s/K), theta) - 1) * (K/theta)), ((pow((x*pow(h, m[i])/K), theta) - 1) * (K/theta)));
    }

    return wrap(a-b-c);
}

which I compile in R (or RStudio) using the command

sourceCpp(<path to file>/vector.cpp)

It compiles - so far so good. However, when I go on to use the function l5 within other R functions it often leads R to crash (both in RStudio and the plain R GUI). In fact, also evaluating it itself isn't any more stable. To reproduce this, e.g. try evaluating l6 multiple times

l6 <- function(zs, ks, ss, Ks, thetas, xs, hs){
    z=zs
    k=ks
    s=ss
    K=Ks
    theta=thetas
    x=xs
    h=hs
    m=0:30
    res <- l5(z, k, s, K, theta,x, h, m)
return(res)
}

and run

l6(0.9, 0.1, 67, 40, 0.5, 44, 1.06)

Specifically, it produces the following runtime error

This application has requested Runtime to terminate it in an unusual way. 

So what is wrong with my function?

Solution

As Dirk suggested below there is an elementary mistake in the for loop, where i runs from 0 to n and thus has n+1 elements, but I intialized only vectors of length n. To avoid this mistake I now implemented the function using iterators

#include <Rcpp.h>

using namespace Rcpp;

// [[Rcpp::export]]
NumericVector l5(double z, double k, double s, double K, double theta, double x, double h, NumericVector m){
    int n = m.size();
    NumericVector a(n);
    NumericVector bu(n);
    NumericVector b(n);
    NumericVector c(n);
    for(NumericVector::iterator i = m.begin(); i != m.end(); ++i){
        a[*i] = pow(z, m[*i]) * (pow((x*pow(h, m[*i])/K), theta) - 1) * (K/theta);
        for(int j=0; j<*i; j++){
            bu[*i] += pow(z,j) * (1 - z) * fmax(((pow((s/K), theta) - 1) * (K/theta)), ((pow((x*pow(h, j)/K), theta) - 1) * (K/theta)));
        }
        b[*i] = k *bu[*i];
        c[*i] = k * pow(z, m[*i]) * fmax(((pow((s/K), theta) - 1) * (K/theta)), ((pow((x*pow(h, m[*i])/K), theta) - 1) * (K/theta)));
    }

    return wrap(a-b-c);
}

Many thanks again!

Answer 1

You are making an elementary C/C++ error:

for(int i=0; i<n+1; i++)

will be accessed n+1 times, but you allocated n spaces.