1

I am trying to run the following recursive row-by-row operation to obtain norm_vec:

#include <Rcpp.h>
#ifdef _OPENMP
#include <omp.h> // OpenMP header
#endif

using namespace Rcpp;
using namespace std;

// [[Rcpp::export]]
NumericVector sim1(NumericVector X) {

  int T = X.length();
  NumericVector norm_vec = NumericVector(T);

  norm_vec[0] = 0;

  int i = 1;

  while (i <= T) {

    NumericVector sim_vec = NumericVector(i);
    NumericVector weight_vec = NumericVector(i);
    NumericVector norm = NumericVector(i);

    int j = 0;
    double sim_vec_sum = 0;

    while (j <= i) {

      sim_vec[j] = exp(-abs(X[i] - X[i-j-1]));
      sim_vec_sum += sim_vec[j];
      j++;

    }

    j = 0;
    double norm_sum = 0;

    while (j <= i) {

      weight_vec[j] = sim_vec[j]/sim_vec_sum;
      norm[j] = X[j]*weight_vec[j];
      norm_sum += norm[j];
      j++;

    }

    norm_vec[i] = norm_sum;

    i++;

  }

  return norm_vec;

}

Calling the code above with sourceCpp("sim1.cpp"), I'd like to get sim1(rnorm(n)). Whereas the code works just fine with smaller n, R completely shuts down once n gets over 19.

I'm not quite comfortable with cpp - probably I made a basic mistake. Thank you so much for your help in advance!

jayc
  • 329
  • 1
  • 8

1 Answers1

2

Loops in R go from 1 to T, but in C/C++ they go from 0 to T-1.

In your loops i goes up to T, and j goes up to i. This way i can reach T, which is more than the maximum of T-1.

Andrey Shabalin
  • 4,389
  • 1
  • 19
  • 18