[Rcpp-devel] Loops, iterators and accumulate

Thu Nov 15 15:11:22 CET 2012

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

-- 
Romain Francois
Professional R Enthusiast
+33(0) 6 28 91 30 30

R Graph Gallery: http://gallery.r-enthusiasts.com
`- http://bit.ly/SweN1Z : SuperStorm Sandy

blog:            http://romainfrancois.blog.free.fr
|- http://bit.ly/RE6sYH : OOP with Rcpp modules
`- http://bit.ly/Thw7IK : Rcpp modules more flexible