[Rcpp-devel] When does using iterators for subscripting help?

[Douglas Bates]

2012/1/6 Hadley Wickham <hadley at rice.edu>:
>> I'm coming late to the party but ...  I was able to squeeze a couple
>> of milliseconds from the computation by expressing the counts as
>> integers and ensuring that I did not copy x by using an Eigen::Map.
>> It may well be that an Rcpp::NumericVector will not be a copy in this
>> case but I have found it difficult to determine exactly when an Rcpp
>> vector is going to be copied.  With Eigen I can ensure that the
>> original storage in R will be used for the vector.
>>
>> count_eigen <- cxxfunction(signature(x="numeric", binwidth="numeric",
>>                                     origin="numeric", nbins="integer"), '
>> double binwidth_ = ::Rf_asReal(binwidth), origin_ = ::Rf_asReal(origin);
>>
>> Eigen::VectorXi counts(::Rf_asInteger(nbins));
>> Eigen::Map<Eigen::VectorXd> x_(as<Eigen::Map<Eigen::VectorXd> >(x));
>>
>> int n = x_.size();
>>
>> counts.setZero();
>> for (int i = 0; i < n; i++) counts[((x_[i] - origin_) / binwidth_)]++;
>>
>> return wrap(counts);
>> ', plugin="RcppEigen")
>>
>>
>> As you see, I use ::Rf_asReal() and ::Rf_asInteger() for conversion to
>> scalar doubles or scalar integers.  Those functions are part of the R
>> API and are fast and general.
>
> Thanks.  Would you mind explaining a couple of the C++ idioms that you used?
>
> * What's going on with Eigen::Map<Eigen::VectorXd>
> x_(as<Eigen::Map<Eigen::VectorXd> >(x)) - you're creating a new Map
> (templated to contain VectorXds) called x_ and initialising with a
> converted x?

Dirk and I wrote an introductory vignette for the RcppEigen package
which, if you have not already done so, would be worth looking at.

An Eigen::VectorXd is a double precision vector with storage allocated
from the C++ heap.  An Eigen::Map<Eigen::Vector> behaves similarly but
maps the storage from some other source, in this case the storage from
R.  Had I followed the recommendations in the vignette mentioned above
I would have declared it as

const Eigen::Map<Eigen::VectorXd>  x_(as<Eigen::Map<Eigen::VectorXd> >(x));

to emphasize that it is a read-only vector.

> * I find it hard to visually parse  <Eigen::Map<Eigen::VectorXd> >(x).
>  Is there a reason you don't write it as
> x_(as<Eigen::Map<Eigen::VectorXd>>(x)?  Maybe because that would then
> form the >> operator?

Yes, exactly.  In production code I usually write something like

typedef Eigen::Map<Eigen::VectorXd>  MVec;

const MVec x_(as<MVec>(x));

As I said earlier, this construction guarantees that the storage is
not copied.  An Rcpp::NumericVector also uses the storage from R but
it can be difficult to determine exactly when a copy constructor,
which actually allocates new storage for an SEXPREC, is invoked.
> * What does ::Rf_asReal mean?  How is it different to Rf_asReal?
>
> Thanks!
>
> Hadley
>
> --
> Assistant Professor / Dobelman Family Junior Chair
> Department of Statistics / Rice University
> http://had.co.nz/