<div dir="ltr">(Based on this Q: <a href="http://stackoverflow.com/q/36206524/164611">http://stackoverflow.com/q/36206524/164611</a> , and available from this gist: <a href="https://gist.github.com/shabbychef/9c2ed6b2ee92f104a8b1">https://gist.github.com/shabbychef/9c2ed6b2ee92f104a8b1</a> )<div><br></div><div>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.</div><div><br></div><div>yunoparallel.cpp:</div><div><br></div><div><div>#include <Rcpp.h></div><div>using namespace Rcpp;</div><div><br></div><div>// [[Rcpp::depends(RcppParallel)]]</div><div>#include <RcppParallel.h></div><div>using namespace RcppParallel;</div><div><br></div><div>// accumulate # of elements seen on the input iterator.</div><div>void accum_nel_int(RVector<double>::const_iterator vecit, </div><div> RVector<double>::const_iterator end, </div><div> RVector<int>::iterator xret) {</div><div> for (;vecit != end;++vecit) { *xret += 1; }</div><div>}</div><div>void join_em_int(RVector<int> lhs, const RVector<int> rhs) {</div><div> lhs[0] += rhs[0];</div><div>}</div><div><br></div><div>struct NEL_int : public Worker</div><div>{ </div><div> // source vector</div><div> const RVector<double> input;</div><div> // accumulated value</div><div> RVector<int> output;</div><div> </div><div> // constructors</div><div> NEL_int(const NumericVector input, IntegerVector output) : input(input), output(output) {}</div><div> NEL_int(const NEL_int& othr, Split) : input(othr.input), output(othr.output) {}</div><div> </div><div> // accumulate just the element of the range I've been asked to</div><div> void operator()(std::size_t begin, std::size_t end) {</div><div> accum_nel_int(input.begin() + begin, input.begin() + end, output.begin());</div><div> }</div><div> </div><div> // join my value with that of another NEL_int</div><div> void join(const NEL_int& rhs) { </div><div> join_em_int(output,rhs.output);</div><div> }</div><div>};</div><div><br></div><div>// dumb example counting # of elements on input:</div><div>// [[Rcpp::export]]</div><div>IntegerVector dumbCount_int(NumericVector x) {</div><div> // allocate the output</div><div> IntegerVector output(1);</div><div> // declare the instance</div><div> NEL_int wrkr(x,output);</div><div> // call parallel_reduce to start the work</div><div> parallelReduce(0, x.length(), wrkr, 5000);</div><div> // return the computed number of non-na elements</div><div> return output;</div><div>}</div><div><br></div><div>//for vim modeline: (do not edit)</div><div>// 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:</div><div><br></div><div>In R:</div><div><br></div><div><div>library(Rcpp)</div><div>source('yunoparallel.cpp')</div><div>x <- rnorm(1e5)</div><div>dumbCount_int(x)</div><div>dumbCount_int(x)</div><div># fine, right? no.</div><div>replicate(10,dumbCount_int(x))</div><div>sd(replicate(100,dumbCount_int(x)))</div></div><div><br></div><div>sessionInfo:</div><div><br></div><div><div>> sessionInfo() </div><div>R version 3.2.4 Revised (2016-03-16 r70336)</div><div>Platform: x86_64-pc-linux-gnu (64-bit)</div><div>Running under: Ubuntu 14.04.4 LTS</div><div><br></div><div>locale:</div><div> [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 </div><div> [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 </div><div><br></div><div>attached base packages:</div><div>[1] graphics grDevices datasets stats utils methods base </div><div><br></div><div>other attached packages:</div><div>[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 </div><div><br></div><div>loaded via a namespace (and not attached):</div><div>[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</div><div>[9] lattice_0.20-33 </div></div><div><br></div><div><br></div>-- <br><div class="gmail_signature"><br><div>--sep</div></div>
</div></div>