RcppParallel and C++. Inconsistent results

Question

I've been playing around with RcppParallel and coded up a fairly simple example to figure out how things work. The code is displayed below.

The function float pdf(double x, double sigma) calculates a scaled version of a Gaussian distribution with mean 0 and standard deviation sigma.

Struct_1 is a struct that creates a worker to perform some calculations. I populate a matrix to figure out why certain things are not working correctly.

void Struct_check() performs the calculations.

The function seems to work but every now and again it does not work as expected. I think that it has to do with the types used to perform the calculations in the function pdf!

An example run is displayed below the code.

I would appreciate any help help!

#include <RcppParallel.h>
#include <RcppArmadillo.h>
#include <RcppArmadilloExtensions/sample.h>
#include <math.h>

#define pi           3.14159265358979323846  /* pi */
using namespace arma;
using namespace Rcpp;
using namespace R;
using namespace sugar;
using namespace std;
using namespace RcppParallel;

// Enable C++11 via this plugin (Rcpp 0.10.3 or later)
// [[Rcpp::plugins(cpp11)]]
// [[Rcpp::depends(RcppParallel)]]


// Returns the probability of x, given the distribution described by mu and sigma.
float pdf(double x,  double sigma)
{
     return exp( -1 * x * x / (2 * sigma * sigma)) / sigma;
}

struct Struct_1 : public Worker
{
     arma::vec wr;
     arma::vec sr; 
     NumericVector w2;

     // source matrix
     const RVector<double> input;

     // destination matrix
     RMatrix<double> output;

     // initialize with source and destination
     Struct_1(const NumericMatrix input, NumericMatrix output) 
          : input(input), output(output) {}

     //what is done. 
     void operator()(std::size_t begin, std::size_t end) {

          for (std::size_t i=begin; i<end; i++){ //the processor loop!

               NumericVector w2(3); 

               for (int comp_j=0; comp_j<3; ++comp_j){
                    w2(comp_j) = wr(comp_j) * pdf( input[i], sr(comp_j) ) ;
               }

               double sw1 = sum(w2);

               output(i,0) = w2(0); 
               output(i,1) = w2(1); 
               output(i,2) = w2(2); 
               output(i,3) = sw1; 

               w2 = w2/sw1;

               output(i,4) = w2(0); 
               output(i,5) = w2(1); 
               output(i,6) = w2(2); 

               double sw2 = sum(w2);

               output(i,7) = sw2;

          }//end of i loop
     }//end of operator
};


// [[Rcpp::depends("RcppArmadillo")]]
// [[Rcpp::export]]
void Struct_check(){

     //Some vecs defined
     arma::vec wr = {0.2522, 0.58523, 0.16257};
     arma::vec s2r = {1.2131, 2.9955, 7.5458}; 
     arma::vec sr = sqrt(s2r);

     //an arma mat that will be used in the struct 
     arma::mat arb_mat;
     arb_mat.randn(20);
     Rcout<<"Arb_mat=\n"<<arb_mat<<endl;

     NumericMatrix r_i_x_NM = as<NumericMatrix>(wrap( arb_mat )); //convert to NumericMatrix
     NumericMatrix output( r_i_x_NM.nrow() , 8 ); //define the output matrix

     Struct_1 struct_1( r_i_x_NM , output);
     struct_1.wr = wr;
     struct_1.sr = sr;

     Rcout<<"nrow output = "<<output.nrow()<<endl;
     Rcout<<"ncol output = "<<output.ncol()<<endl;

     parallelFor(0, r_i_x_NM.length(), struct_1);
     Rcout<<"completed Parallell calculations"<<endl;
     Rcout<<"output = \n"<<output<<endl;

}

Run from within Rstudio. I am running OS X El Capitan if that matter.

Struct_check()
Arb_mat=
  -0.4539
   0.7915
   0.2581
   1.5917
   0.3718
   0.4452
   0.1230
  -1.4719
   0.0024
   2.6166
  -0.4839
  -1.2865
   2.0492
  -1.5980
  -0.7531
  -0.7312
  -1.4482
   0.0202
   0.4434
  -0.0224

nrow output = 20
ncol output = 8
completed Parallell calculations
output = 
0.210336 0.326704 0.0583792 0.595419 0.353256 0.548696 0.0980473 1.00000
0.176872 0.304564 0.0567753 0.538211 0.328629 0.565882 0.105489 1.00000
0.222778 0.334398 0.0589211 0.616097 0.361596 0.542768 0.0956361 1.00000
0.0805904 0.221529 0.0500356 0.352155 0.228849 0.629067 0.142084 1.00000
0.216296 0.330423 0.0586421 0.605361 0.357301 0.545827 0.0968712 1.00000
0.211018 0.327133 0.0584096 0.596561 0.353724 0.548365 0.0979106 1.00000
0.227556 0.337284 0.0591224 0.623962 0.364695 0.540551 0.0947533 1.00000
0.0937487 0.235521 0.0512670 0.380537 0.246359 0.618918 0.134723 1.00000
0.228979 0.338136 0.0591817 0.626297 0.365608 0.539897 0.0944947 1.00000
0.0136216 0.107837 0.0375975 0.159056 0.0856401 0.677981 0.236379 1.00000
0.207911 0.325174 0.0582705 0.591355 0.351584 0.549879 0.0985372 1.00000
0.115751 0.256513 0.0530344 0.425298 0.272164 0.603137 0.124699 1.00000
0.0405607 0.167755 0.0448066 0.253123 0.160241 0.662743 0.177015 1.00000
0.0799309 0.220793 0.0499695 0.350694 0.227922 0.629590 0.142488 1.00000
0.181248 0.307594 0.0569989 0.545841 0.332053 0.563523 0.104424 1.00000
0.183689 0.309265 0.0571216 0.550075 0.333934 0.562222 0.103843 1.00000
**0.228941 0.338113 0.0591801 0.618557 0.591026 0.872861 0.152777 1.61666**
0.228941 0.338113 0.0591801 0.626234 0.365583 0.539915 0.0945016 1.61666
0.211153 0.327218 0.0584156 0.596786 0.353816 0.548300 0.0978837 1.00000
0.228932 0.338108 0.0591798 0.626220 0.365578 0.539919 0.0945032 1.00000

The error occurs when -1.4482 is evaluated to produce the following line 0.228941 0.338113 0.0591801 0.618557 0.591026 0.872861 0.152777 1.61666

In R - checking I get :

wr <- c(0.2522, 0.58523, 0.16257)
s2r <- c(1.2131, 2.9955, 7.5458)
sr <- sqrt(s2r)

w<-NULL
for (i in 1:3){
     w[i] = wr[i]*exp( -0.5*((-1.4482/sr[i])^ 2))/(sr[i])
}

w
[1] 0.09646706 0.23826346 0.05150315

sum(w)
[1] 0.3862337

w = w/sum(w)

w
[1] 0.2497635 0.6168894 0.1333471