Unexpected behaviour in Rcpp

Question

Please note that this error was taken from a bigger context, which I cannot obviously report here entirely.

I have the following functions in the file fun.cpp

#include <RcppArmadilloExtensions/sample.h>

using namespace Rcpp;

// [[Rcpp::depends(RcppArmadillo)]]

arma::vec colMeans(arma::mat data){

   int n_0 = data.n_rows;

   arma::vec xbar(data.n_cols);
   for(int i = 0; i < data.n_rows; i++){
      for(int j = 0; j < data.n_cols; j++){
         xbar[j] += data(i,j) /n_0;
      }
   }

   return xbar;
}

// [[Rcpp::export]]
List PosteriorNIW(arma::mat data, arma::vec mu0, double lambda0,
              double df0, arma::mat V){

   // Compute posterior
   int n = data.n_rows;

   arma::vec xbar = colMeans(data);

   double lambdan = lambda0 + n;

   arma::vec mun = (lambda0 * mu0 + n * xbar) / lambdan;

   arma::mat S;
   S.zeros(data.n_cols, data.n_cols);
   for(int i = 0; i < n; i++){
      S += (arma::conv_to<arma::vec>::from(data.row(i)) - xbar) * arma::trans(arma::conv_to<arma::vec>::from(data.row(i)) - xbar);
   }

   arma::mat Vn = V + S + ((lambda0*n)/(lambda0 + n)) * (xbar - mu0) * arma::trans(xbar - mu0);

   return List::create(_["mun"] = mun,
                  _["Vn"] = Vn,
                  _["lambdan"] = lambdan);
 }

Calling now:

library(Rcpp); library(RcppArmadillo)
mu0 <- c(3,3)
V0 <- matrix(c(2.5,0.0,0.0,2.5), nrow = 2)
sourceCpp("fun.cpp")

data <- cbind(rep(5,15),rep(0,15))
PosteriorNIW(data,  mu0, 1, 1, V0)

gives the expected result.

$mun
    [,1]
[1,] 4.8750
[2,] 0.1875

$Vn
     [,1]    [,2]
[1,]  6.250 -5.6250
[2,] -5.625 10.9375

$lambdan
[1] 16

Now if I add to the file fun.cpp the following functions (again, these are taken from a bigger context so don't bother trying to understand but just paste them) strange things happens:

// [[Rcpp::export]]
NumericMatrix myFun(arma::mat t_dish, arma::cube data){
   int l = 0;
   for(int j = 0; j < data.n_rows; j++){
      l++;
   }
   NumericMatrix Dk(l, 2);
   return Dk;
}

// [[Rcpp::export]]
int myFun2(arma::cube n_cust){

   arma::mat temp = n_cust.subcube(arma::span(0), arma::span(), arma::span());
   int i;
   for(i = 0; i < n_cust.n_cols; i++){
      arma::rowvec temp2 = temp.row(i);
   }

   return i + 1;
}

// [[Rcpp::export]]
arma::vec myFun3(arma::mat k_tables){
   arma::vec temp(k_tables.n_cols * k_tables.n_rows);
   int l = 0;
   if(!R_IsNA(k_tables(0,0))){
      l++;
   }

   arma::vec temp2(l);
   arma::vec tmp3 = sort(temp2);
   return tmp3;
}

double myFun4(arma::vec x, double nu, arma::vec mu, arma::mat Sigma){
   arma::vec product = (arma::trans(x - mu) * arma::inv(Sigma) * (x - mu));
   double num = pow(1 + (1 / nu) * product[0], - ( nu + 2 ) / 2);
   double den = pow(sqrt(M_PI * nu),2) * sqrt(arma::det(Sigma));
   return  num / den;
}

bool myFun5(NumericVector X, double z) { 
   return std::find(X.begin(), X.end(), z)!=X.end(); 
}

calling PosteriorNIW(data, mu0, 1, 1, V0) repeatedly starts giving different results every time. Note that there is no randomness in the functions and that obviously those functions have got no impact as they are not called in the original function.

I have tried on a different machine to make sure it was not a problem of my compiler but the error keeps happening.

I know that removing those function (even just one of them) fixes the problem but clearly this is not a feasible solution when I am working with more functions.

I would like to know if other users are able to replicate this behavior and if yes if there is a fix for it.

Thank you in advance

EDIT:

The version of R is 3.3.2 and Rtools is 3.4. Both Rcpp and RcppArmadillo are up-to-date

This is *not* reproducible. I get the same result (from calling `PosteriorNIW(data, mu0, 1, 1, V0) ` twice) whether or not those extra functions are included. — Dirk Eddelbuettel, Feb 15 '17 at 18:24
I will try again on a different machine. Any possible reason why this might be happening on mine? (Updating rtools or things like that). However please not that the error happens only if `data` is not declared again — adaien, Feb 15 '17 at 19:15
Please read about [creating a good reproducible example](http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example). — Dirk Eddelbuettel, Feb 15 '17 at 21:41
I provided all is needed to run the code, but if the error depends on the machine used unfortunately there is nothing I can do about. I tried on two different ones and the error is the same. — adaien, Feb 15 '17 at 22:53
I get different results on runs even *before* I add your second block of functions to the file. Can you paste in the results you get from the first "expected" code block, and make absolute sure they are repeated? — Spacedman, Feb 16 '17 at 10:00
Thank you @Spacedman for trying. I made an edit where I show the results. I have tried again and the first block always gives the right answer. Can you please check if the answer is correct if you call the function after initializing `data` again? — adaien, Feb 16 '17 at 10:08

score 2 · Accepted Answer · answered Feb 16 '17 at 10:21

You're not zeroing xbar in your colMeans function. If I do do that:

arma::vec colMeans(arma::mat data){

   int n_0 = data.n_rows;

   arma::vec xbar;
   xbar.zeros(data.n_cols);
   for(int i = 0; i < data.n_rows; i++){
      for(int j = 0; j < data.n_cols; j++){
         xbar[j] += data(i,j) /n_0;
      }
   }

   return xbar;
}

I get this everytime:

> PosteriorNIW(data,  mu0, 1, 1.1, V0)
$mun
       [,1]
[1,] 4.8750
[2,] 0.1875

$Vn
       [,1]    [,2]
[1,]  6.250 -5.6250
[2,] -5.625 10.9375

$lambdan
[1] 16

Even when I do add your extra block of code.

I don't know if these vectors are documented to be initialised to zero by their constructor (in which case this might be a bug there) or not, in which case its your bug!

My debugging technique in these situations is to print out things as they are calculated and finding out where things are going off. That led me to the `colMeans` function and then I spotted the unzeroed vector. — Spacedman, Feb 16 '17 at 10:51
I do the same usually, but when I was trying to print out `n` in the function `PosteriorNIW`, just adding the printing made the function works and removed the error.. Cheers — adaien, Feb 16 '17 at 11:00

Unexpected behaviour in Rcpp

1 Answers1