[Rcpp-devel] Minimal Working Example for RcppParallel gives nonsensical results

[steven pav]

(Based on this Q: http://stackoverflow.com/q/36206524/164611 , and
available from this gist:
https://gist.github.com/shabbychef/9c2ed6b2ee92f104a8b1 )

I am trying to write some code using the parallelReduce framework from
RcppParallel. Somehow I must be violating thread safety. I reduced my real
code down to a simple MWE that does not work. I'll include it below. This
very silly example just computes the length of an input vector. I thought
it was working correctly, but noticed that if I ran it many times quickly
(via replicate in R), I would get different answers. Per Dirk's advice, I
tried to remove any hint of Rcpp functions to maintain thread safety, but I
am not sure what is left. Any help appreciated.

yunoparallel.cpp:

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::depends(RcppParallel)]]
#include <RcppParallel.h>
using namespace RcppParallel;

// accumulate # of elements seen on the input iterator.
void accum_nel_int(RVector<double>::const_iterator vecit,
                   RVector<double>::const_iterator end,
                   RVector<int>::iterator xret) {
    for (;vecit != end;++vecit) { *xret += 1; }
}
void join_em_int(RVector<int> lhs, const RVector<int> rhs) {
    lhs[0] += rhs[0];
}

struct NEL_int : public Worker
{
   // source vector
   const RVector<double> input;
   // accumulated value
   RVector<int> output;

   // constructors
   NEL_int(const NumericVector input, IntegerVector output) : input(input),
output(output) {}
   NEL_int(const NEL_int& othr, Split) : input(othr.input),
output(othr.output) {}

   // accumulate just the element of the range I've been asked to
   void operator()(std::size_t begin, std::size_t end) {
       accum_nel_int(input.begin() + begin, input.begin() + end,
output.begin());
   }

   // join my value with that of another NEL_int
   void join(const NEL_int& rhs) {
       join_em_int(output,rhs.output);
   }
};

// dumb example counting # of elements on input:
// [[Rcpp::export]]
IntegerVector dumbCount_int(NumericVector x) {
   // allocate the output
   IntegerVector output(1);
   // declare the instance
   NEL_int wrkr(x,output);
   // call parallel_reduce to start the work
   parallelReduce(0, x.length(), wrkr, 5000);
   // return the computed number of non-na elements
   return output;
}

//for vim modeline: (do not edit)
//
vim:ts=4:sw=4:tw=79:fdm=marker:fmr=FOLDUP,UNFOLD:cms=//%s:syn=cpp:ft=cpp:ai:si:et:cin:nu:fo=croql:cino=p0t0c5(0:

In R:

library(Rcpp)
source('yunoparallel.cpp')
x <- rnorm(1e5)
dumbCount_int(x)
dumbCount_int(x)
# fine, right? no.
replicate(10,dumbCount_int(x))
sd(replicate(100,dumbCount_int(x)))

sessionInfo:

> sessionInfo()
R version 3.2.4 Revised (2016-03-16 r70336)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 14.04.4 LTS

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8
LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                  LC_ADDRESS=C
          LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8
LC_IDENTIFICATION=C

attached base packages:
[1] graphics  grDevices datasets  stats     utils     methods   base

other attached packages:
[1] Rcpp_0.12.3     colorout_1.1-2  fortunes_1.5-2  Quandl_2.7.0
 xts_0.9-7       zoo_1.7-12      devtools_1.10.0 drat_0.1.0

loaded via a namespace (and not attached):
[1] httr_1.1.0          R6_2.1.2            tools_3.2.4
memoise_1.0.0       grid_3.2.4          jsonlite_0.9.19     digest_0.6.9
     RcppParallel_4.3.15
[9] lattice_0.20-33


-- 

--sep
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20160325/80a55af3/attachment.html>