[Rcpp-devel] Loops, iterators and accumulate
Rainer M Krug
r.m.krug at gmail.com
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-----
More information about the Rcpp-devel
mailing list