RCPP and R discrepancy

Question

I'm new to C++ programming and apologize if my solution is in plain sight. I am attempting to use RCPP to speed up a slow R function. I think I've narrowed down the issue to a nested for loop. I've simplified the function and provided one R and one RCPP version for comparison. Will someone please explain why my RCPP function yields different results? Thanks!

## Data ##

set.seed(666)
input <- rmultinom(10,2,c(.4,.5,.6)) + 1

## R ##

testR <- \(input){
  M1 <- matrix(c(0.5,0.4,0.0,0.3,0.5,0.0,0.2,0.1,1.0),3,3)
  M2 <- matrix(c(0.75,0.0,0.0,0.0,0.6,0.0,0.25,0.4,1.0),3,3)
  
  Mrows <- nrow(M1)
  tmsteps <- ncol(input)
  N <- nrow(input)
  
  alphas <- NULL; tmp <- NULL; out <- NULL
  for(i in 1:N){
    alphas = c(0,-1e6,-1e6)
    for(j in 1:tmsteps){
      for(k in 1:Mrows){
        tmp[k] = sum(alphas + M1[,k] + M2[k, input[i,j] ])
      }
      alphas <- tmp
    }
    out[i] <- sum(alphas)
  }
  sum(out)
  }

## RCPP ##

cppFunction('double testRCPP(IntegerMatrix input){
    NumericVector v1 = {0.5,0.4,0.0,0.3,0.5,0.0,0.2,0.1,1.0};
    v1.attr("dim") = Dimension(3, 3);
    NumericMatrix M1 = as<NumericMatrix>(v1);
    
    NumericVector v2 = {0.75,0.0,0.0,0.0,0.6,0.0,0.25,0.4,1.0};
    v2.attr("dim") = Dimension(3, 3);
    NumericMatrix M2 = as<NumericMatrix>(v2);
    
    int Mrows = M1.nrow();
    int tmsteps = input.ncol();
    int N = input.nrow();
    
    NumericVector alphas(3);
    NumericVector tmp(3);
    NumericVector out(N);
    
    for(int i=0; i<N; i++){
    alphas = {0,-1e6,-1e6};
      for(int j=0; j<tmsteps; j++){
        for(int k=0; k<Mrows; k++){
        tmp[k] = sum(alphas + M1(_,k) + M2(k, (input(i,j) - 1) ));
        }
      alphas = tmp;
    }
      out += alphas;
    }
    return(sum(out));
  }')

> testRCPP(input)
[1] -2.273726e+14
> testR(input)
[1] -354293536945

Answer 1

I have figured out how to get the Rcpp to behave like the R function. I think my issue has to do with C++ variable scoping.

I had previously been initializing the tmp variable outside the nested for loop.

NumericVector tmp(3);

    for(int i=0; i<N; i++){
    alphas = {0,-1e6,-1e6};
    ...

All is good when I declare the tmp variable inside the loop, although I don't understand why yet.

    for(int i=0; i<N; i++){
    alphas = {0,-1e6,-1e6};
      for(int j=0; j<tmsteps; j++){
      NumericVector tmp(3);
        for(int k=0; k<Mrows; k++){
        tmp[k] = sum(alphas + M1(_,k) + M2(k, (input(i,j) - 1) ));
        }
      alphas = tmp;
    }
...