[Rcpp-devel] Loops, iterators and accumulate

Thu Nov 15 15:24:04 CET 2012

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

Just to add one more test here:

Unit: microseconds
      expr    min      lq median      uq     max
1  sum1(x) 19.835 20.0445 20.184 20.3585  43.651
2 sum2a(x) 19.835 20.0095 20.184 20.3240  48.958
3  sum2(x) 93.797 94.0760 94.146 94.4250 204.355
4  sum3(x) 98.546 98.7205 98.825 99.4540 255.760
5  sum4(x) 19.835 19.9750 20.114 20.5685  43.720
6   sum(x) 16.064 16.3430 16.413 16.7270  50.565

Cheers,

Rainer

On 15/11/12 15:11, Romain Francois wrote:
> I'm not getting the same results, so I might not make the same conclusions:
> 
> Unit: microseconds expr    min     lq  median      uq    max 1   sum(x) 16.072 16.161 16.2090
> 16.3110 21.721 2  sum1(x) 10.771 10.954 11.0445 11.1910 16.163 3  sum2(x) 48.438 48.567 48.6650
> 48.8425 61.116 4 sum2a(x) 10.806 10.965 11.0695 11.2080 17.401 5  sum3(x) 10.835 10.967 11.0490
> 11.1875 23.943 6  sum4(x) 10.847 10.989 11.0690 11.2160 16.331
> 
> For me that makes more sense.
> 
> The code that gets called by the sugar version is this:
> 
> // RTYPE = REALSXP template <bool NA, typename T> class Sum<REALSXP,NA,T> : public Lazy< double
> , Sum<REALSXP,NA,T> > { public: typedef typename Rcpp::VectorBase<REALSXP,NA,T> VEC_TYPE ; 
> typedef typename Rcpp::traits::Extractor< REALSXP, NA, T>::type VEC_EXT ;
> 
> Sum( const VEC_TYPE& object_ ) : object(object_.get_ref()){}
> 
> double get() const { double result = 0 ; int n = object.size() ; for( int i=0; i<n; i++){ 
> result += object[i] ; } return result ; } private: const VEC_EXT& object ; } ;
> 
> so about the same logic as sum1.
> 
> 
> 
> I would have expected accumulate to be better since it works directly with pointers :
> 
>> evalCpp( " DEMANGLE( NumericVector::iterator )" )
> [1] "double*"
> 
> 
> 
> 
> Le 15/11/12 14:24, Hadley Wickham a écrit :
>> Hi all,
>> 
>> I'm trying to get a better handle on writing efficient numerical code with C++.  I've
>> attached five implementations of a simple sum function (mostly ignoring NAs) and timed them
>> below:
>> 
>>> source("sum.r")
>> Unit: microseconds expr   min    lq median    uq  max neval sum(x) 11.67 11.72  11.76 11.82
>> 20.6   100 sum1(x)  9.58  9.68   9.75  9.82 19.7   100 sum2(x) 37.44 37.55  37.59 37.63 48.1
>> 100 sum2a(x)  9.59  9.66   9.73  9.81 14.0   100 sum3(x) 31.84 31.92  31.98 32.03 35.6   100 
>> sum4(x)  9.60  9.69   9.77  9.86 30.5   100
>> 
>> * sum: R's built in sum * sum1: loop * sum2: iterator * sum2a: iterator, but only compute
>> x.end() once * sum3: use accumulate * sum4: use sugar sum
>> 
>> My questions:
>> 
>> * the R cpp introduction mentions that iterators are usually faster than using [ directly,
>> but I don't see that here (or in a number of other places I've tried it). Why?
> 
> We made progress since we wrote this.
> 
> operator[] still is conceptually more complex:
> 
> but with proper inlining the compiler should be able to pretty much generate the same code.
> 
>> * The introduction vignette also uses iterators in a (as far as I can tell) non-standard way
>> - is that a mistake or some C++ magic I don't understand?  (I've used them in standard style
>> in the attached code)
> 
> What do you mean ? Can you share the code you think is not standard.
> 
>> * why is accumulate so slow?
> 
> no idea. Not what I would have expected. I would have guessed it would be better since there
> could have been things like loop unrolling under the hood etc ...
> 
> maybe the implementation of the STL you use is not good. What platform are you running this on.
> The results above come from OS X Lion (not Mountain Lion), with
> 
> $ g++ --version i686-apple-darwin11-llvm-g++-4.2 (GCC) 4.2.1 (Based on Apple Inc. build 5658)
> (LLVM build 2335.15.00) Copyright (C) 2007 Free Software Foundation, Inc. This is free
> software; see the source for copying conditions.  There is NO warranty; not even for
> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
> 
>> * I'm slightly surprised that sugar::sum is about 20% faster than R's built-in sum. I would
>> have expected base sum to be fairly well optimised and close to the metal. I guess this may
>> be because R's sum deals with a few more cases of mixed arg types in ...
>> 
>> Hadley
> 
> Hmm. Not easy to find my way through the R forest of code. Once I navigate in do_summary ...
> you end up here:
> 
> static Rboolean rsum(double *x, R_xlen_t n, double *value, Rboolean narm) { LDOUBLE s = 0.0; 
> R_xlen_t i; Rboolean updated = FALSE;
> 
> for (i = 0; i < n; i++) { if (!narm || !ISNAN(x[i])) { if(!updated) updated = TRUE; s += x[i]; 
> } } *value = (double) s;
> 
> return(updated); }
> 
> So it takes care of the na.rm argument and NaN values. Hmm, I guess we should too. But we get 
> consistent results:
> 
>> x <- c(1, 2, Inf ) sum(x)
> [1] Inf
>> sum1(x)
> [1] Inf
>> sum2(x)
> [1] Inf
>> sum2a(x)
> [1] Inf
>> sum3(x)
> [1] Inf
>> sum4(x)
> [1] Inf
> 
> 
>> x <- c(1, 2, NA ) sum(x)
> [1] NA
>> sum1(x)
> [1] NA
>> sum2(x)
> [1] NA
>> sum2a(x)
> [1] NA
>> sum3(x)
> [1] NA
>> sum4(x)
> [1] NA
> 
> So, not sure why we should pay the price of this if test in the R version.
> 
> Romain
> 

-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.4.11 (GNU/Linux)
Comment: Using GnuPG with Mozilla - http://www.enigmail.net/

iEYEARECAAYFAlCk+wQACgkQoYgNqgF2egrGOwCfUPy2o22OjXXON94bNqxYMdB0
3K0AnicG44UVgg8Tst3Za015V/ZTZ75g
=kN09
-----END PGP SIGNATURE-----