2012/1/6 Hadley Wickham <had...@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/ _______________________________________________ Rcpp-devel mailing list Rcpp-devel@lists.r-forge.r-project.org https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
