Re: [Rcpp-devel] Optimising 2d convolution
On 7 January 2012 at 11:40, baptiste auguie wrote:
| Hi,
|
| (Rcpp)Armadillo has a 1D convolution function <
| http://arma.sourceforge.net/docs.html#conv >. It might be a very
| inefficient way, to use it sequentially for 2D; I haven't tried.
Good point. Off-list, I already suggested to Hadley to look at what the image
processing folks have done. And as Conrad is now proud author of several
published papers in the area (kudos!), he may have a pointers. Conrad?
Dirk
| HTH,
|
| baptiste
|
|
|
| On 7 January 2012 07:43, Hadley Wickham wrote:
| > Hi all,
| >
| > The rcpp-devel list has been very helpful to me so far, so I hope you
| > don't mind another question! I'm trying to speed up my 2d convolution
| > function:
| >
| >
| >
| > library(inline)
| > # 2d convolution
-
| >
| > convolve_2d <- cxxfunction(signature(sampleS = "numeric", kernelS =
| > "numeric"), plugin = "Rcpp", '
| > Rcpp::NumericMatrix sample(sampleS), kernel(kernelS);
| > int x_s = sample.nrow(), x_k = kernel.nrow();
| > int y_s = sample.ncol(), y_k = kernel.ncol();
| >
| > Rcpp::NumericMatrix output(x_s + x_k - 1, y_s + y_k - 1);
| > for (int row = 0; row < x_s; row++) {
| > for (int col = 0; col < y_s; col++) {
| > for (int i = 0; i < x_k; i++) {
| > for (int j = 0; j < y_k; j++) {
| > output(row + i, col + j) += sample(row, col) * kernel(i, j);
| > }
| > }
| > }
| > }
| > return output;
| > ')
| >
| >
| > x <- diag(1000)
| > k <- matrix(runif(20* 20), ncol = 20)
| > system.time(convolve_2d(x, k))
| > # user system elapsed
| > # 14.759 0.028 15.524
| >
| > I have a vague idea that to get better performance I need to get
| > closer to bare pointers, and I need to be careful about the order of
| > my loops to ensure that I'm working over contiguous chunks of memory
| > as often as possible, but otherwise I've basically exhausted my
| > C++/Rcpp knowledge. Where should I start looking to improve the
| > performance of this function?
| >
| > The example data basically matches the real problem - x is not usually
| > going to be much bigger than 1000 x 1000 and k typically will be much
| > smaller. (And hence, based on what I've read, this technique should
| > be faster than doing it via a discrete fft)
| >
| > 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
| ___
| 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
--
"Outside of a dog, a book is a man's best friend. Inside of a dog, it is too
dark to read." -- Groucho Marx
___
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
Re: [Rcpp-devel] Optimising 2d convolution
On 7 January 2012 11:40, baptiste auguie wrote:
> Hi,
>
> (Rcpp)Armadillo has a 1D convolution function <
> http://arma.sourceforge.net/docs.html#conv >. It might be a very
> inefficient way, to use it sequentially for 2D; I haven't tried.
Also, it only works for separable kernels.
baptiste
>
> HTH,
>
> baptiste
>
>
>
> On 7 January 2012 07:43, Hadley Wickham wrote:
>> Hi all,
>>
>> The rcpp-devel list has been very helpful to me so far, so I hope you
>> don't mind another question! I'm trying to speed up my 2d convolution
>> function:
>>
>>
>>
>> library(inline)
>> # 2d convolution
>> -
>>
>> convolve_2d <- cxxfunction(signature(sampleS = "numeric", kernelS =
>> "numeric"), plugin = "Rcpp", '
>> Rcpp::NumericMatrix sample(sampleS), kernel(kernelS);
>> int x_s = sample.nrow(), x_k = kernel.nrow();
>> int y_s = sample.ncol(), y_k = kernel.ncol();
>>
>> Rcpp::NumericMatrix output(x_s + x_k - 1, y_s + y_k - 1);
>> for (int row = 0; row < x_s; row++) {
>> for (int col = 0; col < y_s; col++) {
>> for (int i = 0; i < x_k; i++) {
>> for (int j = 0; j < y_k; j++) {
>> output(row + i, col + j) += sample(row, col) * kernel(i, j);
>> }
>> }
>> }
>> }
>> return output;
>> ')
>>
>>
>> x <- diag(1000)
>> k <- matrix(runif(20* 20), ncol = 20)
>> system.time(convolve_2d(x, k))
>> # user system elapsed
>> # 14.759 0.028 15.524
>>
>> I have a vague idea that to get better performance I need to get
>> closer to bare pointers, and I need to be careful about the order of
>> my loops to ensure that I'm working over contiguous chunks of memory
>> as often as possible, but otherwise I've basically exhausted my
>> C++/Rcpp knowledge. Where should I start looking to improve the
>> performance of this function?
>>
>> The example data basically matches the real problem - x is not usually
>> going to be much bigger than 1000 x 1000 and k typically will be much
>> smaller. (And hence, based on what I've read, this technique should
>> be faster than doing it via a discrete fft)
>>
>> 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
___
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
Re: [Rcpp-devel] Optimising 2d convolution
Hi,
(Rcpp)Armadillo has a 1D convolution function <
http://arma.sourceforge.net/docs.html#conv >. It might be a very
inefficient way, to use it sequentially for 2D; I haven't tried.
HTH,
baptiste
On 7 January 2012 07:43, Hadley Wickham wrote:
> Hi all,
>
> The rcpp-devel list has been very helpful to me so far, so I hope you
> don't mind another question! I'm trying to speed up my 2d convolution
> function:
>
>
>
> library(inline)
> # 2d convolution -
>
> convolve_2d <- cxxfunction(signature(sampleS = "numeric", kernelS =
> "numeric"), plugin = "Rcpp", '
> Rcpp::NumericMatrix sample(sampleS), kernel(kernelS);
> int x_s = sample.nrow(), x_k = kernel.nrow();
> int y_s = sample.ncol(), y_k = kernel.ncol();
>
> Rcpp::NumericMatrix output(x_s + x_k - 1, y_s + y_k - 1);
> for (int row = 0; row < x_s; row++) {
> for (int col = 0; col < y_s; col++) {
> for (int i = 0; i < x_k; i++) {
> for (int j = 0; j < y_k; j++) {
> output(row + i, col + j) += sample(row, col) * kernel(i, j);
> }
> }
> }
> }
> return output;
> ')
>
>
> x <- diag(1000)
> k <- matrix(runif(20* 20), ncol = 20)
> system.time(convolve_2d(x, k))
> # user system elapsed
> # 14.759 0.028 15.524
>
> I have a vague idea that to get better performance I need to get
> closer to bare pointers, and I need to be careful about the order of
> my loops to ensure that I'm working over contiguous chunks of memory
> as often as possible, but otherwise I've basically exhausted my
> C++/Rcpp knowledge. Where should I start looking to improve the
> performance of this function?
>
> The example data basically matches the real problem - x is not usually
> going to be much bigger than 1000 x 1000 and k typically will be much
> smaller. (And hence, based on what I've read, this technique should
> be faster than doing it via a discrete fft)
>
> 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
___
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
Re: [Rcpp-devel] Optimising 2d convolution
On Fri, Jan 6, 2012 at 2:58 PM, Douglas Bates wrote:
> On Fri, Jan 6, 2012 at 2:02 PM, Hadley Wickham wrote:
>>> What are you doing the timing on? On a modest desktop (2.6 GHz Athlon
>>> X4) I get less than a second for this
>> ...
k <- matrix(runif(20* 20), ncol = 20)
system.time(convolve_2d(x, k))
>>> user system elapsed
>>> 0.864 0.000 0.862
>>
>> Ooops, I must have been using a old k - I now get
>>
>>> system.time(convolve_2d(x, k))
>> user system elapsed
>> 2.265 0.006 2.298
>>
>> I'm on a 2nd gen macbook air, so the processor isn't the fastest (2.13
>> GHz). Maybe my intuition is wrong but it seems like you should be able
>> to do it faster than that.
>>
>> But let me think - I'm doing 1000 * 1000 * 20 * 20 = 4e08 multiplies +
>> adds in ~2s. So that's 5e-9 s = 5 ns per operation. Looking at it
>> that way it does seem pretty fast already.
>
> You can make this about 40% faster by using block operations in Eigen
> through RcppEigen. In this code fragment I did remember to use a
> typedef and a const modifier. I allocated the result as an
> Rcpp::NumericMatrix then mapped it to avoid an extra allocation and
> copy but that is a minor speedup for matrices of this size.
I just saw a superfluous statement after I changed to allocating an
Rcpp::NumericMatrix. The .setZero() is no longer needed because Rcpp
zeros the storage.
R version 2.14.1 (2011-12-22)
Copyright (C) 2011 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
Natural language support but running in an English locale
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> library(inline)
> # 2d convolution -
>
> convolve_2d <- cxxfunction(signature(sampleS = "numeric", kernelS =
+ "numeric"), plugin = "Rcpp", '
+Rcpp::NumericMatrix sample(sampleS), kernel(kernelS);
+int x_s = sample.nrow(), x_k = kernel.nrow();
+int y_s = sample.ncol(), y_k = kernel.ncol();
+
+Rcpp::NumericMatrix output(x_s + x_k - 1, y_s + y_k - 1);
+for (int row = 0; row < x_s; row++) {
+ for (int col = 0; col < y_s; col++) {
+for (int i = 0; i < x_k; i++) {
+ for (int j = 0; j < y_k; j++) {
+output(row + i, col + j) += sample(row, col) * kernel(i, j);
+ }
+}
+ }
+}
+return output;
+ ')
>
> convolve_2d_Eigen <- cxxfunction(signature(sampleS = "numeric", kernelS =
+ "numeric"), plugin = "RcppEigen", '
+typedef Eigen::Map MMat;
+
+const MMat sample(as(sampleS)), kernel(as(kernelS));
+int x_s = sample.rows(), x_k = kernel.rows();
+int y_s = sample.cols(), y_k = kernel.cols();
+int x_o = x_s + x_k - 1, y_o = y_s + y_k - 1;
+
+Rcpp::NumericMatrix out(x_o, y_o);
+MMat output(out.begin(), x_o, y_o);
+
+for (int row = 0; row < x_s; row++) {
+ for (int col = 0; col < y_s; col++) {
+output.block(row, col, x_k, y_k) += sample(row, col) * kernel;
+ }
+}
+return out;
+ ')
Loading required package: lattice
Attaching package: ‘Matrix’
The following object(s) are masked from ‘package:base’:
det
>
> x <- diag(1000)
> set.seed(1)
> k <- matrix(runif(20* 20), ncol = 20)
> system.time(loop <- convolve_2d(x, k))
user system elapsed
0.840 0.004 0.847
> system.time(block <- convolve_2d_Eigen(x, k))
user system elapsed
0.580 0.004 0.584
> all.equal(loop, block)
[1] TRUE
>
> library(microbenchmark)
> microbenchmark(loop = convolve_2d(x, k),
+block = convolve_2d_Eigen(x, k), times=30)
Unit: milliseconds
expr min lq median uq max
1 block 473.0885 476.9951 478.4295 481.5343 492.4735
2 loop 830.6481 833.8049 835.1170 841.9029 854.5755
>
> proc.time()
user system elapsed
51.726 1.180 52.486
___
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
Re: [Rcpp-devel] Optimising 2d convolution
On Fri, Jan 6, 2012 at 2:02 PM, Hadley Wickham wrote:
>> What are you doing the timing on? On a modest desktop (2.6 GHz Athlon
>> X4) I get less than a second for this
> ...
>>> k <- matrix(runif(20* 20), ncol = 20)
>>> system.time(convolve_2d(x, k))
>> user system elapsed
>> 0.864 0.000 0.862
>
> Ooops, I must have been using a old k - I now get
>
>> system.time(convolve_2d(x, k))
> user system elapsed
> 2.265 0.006 2.298
>
> I'm on a 2nd gen macbook air, so the processor isn't the fastest (2.13
> GHz). Maybe my intuition is wrong but it seems like you should be able
> to do it faster than that.
>
> But let me think - I'm doing 1000 * 1000 * 20 * 20 = 4e08 multiplies +
> adds in ~2s. So that's 5e-9 s = 5 ns per operation. Looking at it
> that way it does seem pretty fast already.
You can make this about 40% faster by using block operations in Eigen
through RcppEigen. In this code fragment I did remember to use a
typedef and a const modifier. I allocated the result as an
Rcpp::NumericMatrix then mapped it to avoid an extra allocation and
copy but that is a minor speedup for matrices of this size.
R version 2.14.1 (2011-12-22)
Copyright (C) 2011 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
Natural language support but running in an English locale
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> library(inline)
> # 2d convolution -
>
> convolve_2d <- cxxfunction(signature(sampleS = "numeric", kernelS =
+ "numeric"), plugin = "Rcpp", '
+Rcpp::NumericMatrix sample(sampleS), kernel(kernelS);
+int x_s = sample.nrow(), x_k = kernel.nrow();
+int y_s = sample.ncol(), y_k = kernel.ncol();
+
+Rcpp::NumericMatrix output(x_s + x_k - 1, y_s + y_k - 1);
+for (int row = 0; row < x_s; row++) {
+ for (int col = 0; col < y_s; col++) {
+for (int i = 0; i < x_k; i++) {
+ for (int j = 0; j < y_k; j++) {
+output(row + i, col + j) += sample(row, col) * kernel(i, j);
+ }
+}
+ }
+}
+return output;
+ ')
>
> convolve_2d_Eigen <- cxxfunction(signature(sampleS = "numeric", kernelS =
+ "numeric"), plugin = "RcppEigen", '
+typedef Eigen::Map MMat;
+
+const MMat sample(as(sampleS)), kernel(as(kernelS));
+int x_s = sample.rows(), x_k = kernel.rows();
+int y_s = sample.cols(), y_k = kernel.cols();
+int x_o = x_s + x_k - 1, y_o = y_s + y_k - 1;
+
+Rcpp::NumericMatrix out(x_o, y_o);
+MMat output(out.begin(), x_o, y_o);
+output.setZero(); // Must explicitly zero out Eigen matrices
+for (int row = 0; row < x_s; row++) {
+ for (int col = 0; col < y_s; col++) {
+output.block(row, col, x_k, y_k) += sample(row, col) * kernel;
+ }
+}
+return out;
+ ')
Loading required package: lattice
Attaching package: ‘Matrix’
The following object(s) are masked from ‘package:base’:
det
>
> x <- diag(1000)
> set.seed(1)
> k <- matrix(runif(20* 20), ncol = 20)
> system.time(loop <- convolve_2d(x, k))
user system elapsed
0.840.000.84
> system.time(block <- convolve_2d_Eigen(x, k))
user system elapsed
0.588 0.000 0.587
> all.equal(loop, block)
[1] TRUE
>
> library(microbenchmark)
> microbenchmark(loop = convolve_2d(x, k),
+block = convolve_2d_Eigen(x, k), times=30)
Unit: milliseconds
expr min lq median uq max
1 block 479.5808 481.5312 482.7564 492.0526 495.4081
2 loop 839.4510 844.4272 846.7940 854.1737 859.0746
>
> proc.time()
user system elapsed
52.242 1.180 53.027
___
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
Re: [Rcpp-devel] Idiom for accessing scalars
On 6 January 2012 at 13:46, Douglas Bates wrote: | I didn't make myself clear. What I meant was that it is not possible | to use asInteger in Rcpp and count on the name being remapped to | Rf_asInteger. Right. I had to enforce that a long time ago because the non-prefixed variants too often clash with things commonly used in other places. One example was length(), and I think I even had issues with error(). Hence the use of the remap toggle and the need for ugliness via Rf_*. Dirk -- "Outside of a dog, a book is a man's best friend. Inside of a dog, it is too dark to read." -- Groucho Marx ___ 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
Re: [Rcpp-devel] Optimising 2d convolution
> What are you doing the timing on? On a modest desktop (2.6 GHz Athlon > X4) I get less than a second for this ... >> k <- matrix(runif(20* 20), ncol = 20) >> system.time(convolve_2d(x, k)) > user system elapsed > 0.864 0.000 0.862 Ooops, I must have been using a old k - I now get > system.time(convolve_2d(x, k)) user system elapsed 2.265 0.006 2.298 I'm on a 2nd gen macbook air, so the processor isn't the fastest (2.13 GHz). Maybe my intuition is wrong but it seems like you should be able to do it faster than that. But let me think - I'm doing 1000 * 1000 * 20 * 20 = 4e08 multiplies + adds in ~2s. So that's 5e-9 s = 5 ns per operation. Looking at it that way it does seem pretty fast already. 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
Re: [Rcpp-devel] Optimising 2d convolution
> Fancy indexing tricks aside, another avenue you might want to look > into is doing the convolution in fourier space ... I have looked into that, and I'm pretty sure for my case (k much smaller than x), that it's less efficient (and less numerically stable). In my experience, it also make the code much harder to understand - convolutions in fourier space are circular, which means you have to add a whole lot of extra padding and it's non-obvious how to debug. I spent several months going down the fft route, and I've already made more progress in a week of using loops. 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
Re: [Rcpp-devel] When does using iterators for subscripting help?
> 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. But that's half the fun of programming in R! ;) Thanks for the other explanations. 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
Re: [Rcpp-devel] Idiom for accessing scalars
On Fri, Jan 6, 2012 at 1:12 PM, Dirk Eddelbuettel wrote: > > On 6 January 2012 at 12:59, Douglas Bates wrote: > | On Fri, Jan 6, 2012 at 12:39 PM, John Chambers > wrote: > | > The "Rf_" part of the API in particular is ugly and somewhat of an add-on > | > forced in a few examples by the use of some common names in the macro > files. > | > | But, as it stands, that is a requirement when using Rcpp. > > Where? I can think of one propagated use, which is at the bottom of the > try/catch structure where we use ::Rf_error. But the commonly used macros > hide it, and we could/should obviously wrap this. > > Otherwise, and especially since the 'Rcpp sugar' initiative took off, I don't > really touch any ::Rf_* myself anymore. Inside the Rcpp code base, sure. But > not really in user-facing stuff and Rcpp applications. I didn't make myself clear. What I meant was that it is not possible to use asInteger in Rcpp and count on the name being remapped to Rf_asInteger. > | I think of the Rf_ part as more due to the fact that C doesn't have a > | concept of namespaces so anything in the R API is at the top level > | namespace leading to some conflicts. > > Agreed. But speaking stylistically, for the same reason that we prefer C++ > versions of C header files (eg cstdint over stdint.h, cstdio over stdio.h, > ...) I am with John on the preference for C++ idioms when given a choice. I suppose I could have just checked whether Rcpp::as calls Rf_asInteger. If so, everything is cool. Unfortunately, I haven't been able to find that specialization. ___ 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
Re: [Rcpp-devel] When does using iterators for subscripting help?
2012/1/6 Hadley Wickham : >> 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 x_(as >(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 > x_(as >(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 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 x_(as >(x)); to emphasize that it is a read-only vector. > * I find it hard to visually parse >(x). > Is there a reason you don't write it as > x_(as>(x)? Maybe because that would then > form the >> operator? Yes, exactly. In production code I usually write something like typedef Eigen::Map MVec; const MVec x_(as(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
Re: [Rcpp-devel] Idiom for accessing scalars
On 6 January 2012 at 12:59, Douglas Bates wrote: | On Fri, Jan 6, 2012 at 12:39 PM, John Chambers wrote: | > The "Rf_" part of the API in particular is ugly and somewhat of an add-on | > forced in a few examples by the use of some common names in the macro files. | | But, as it stands, that is a requirement when using Rcpp. Where? I can think of one propagated use, which is at the bottom of the try/catch structure where we use ::Rf_error. But the commonly used macros hide it, and we could/should obviously wrap this. Otherwise, and especially since the 'Rcpp sugar' initiative took off, I don't really touch any ::Rf_* myself anymore. Inside the Rcpp code base, sure. But not really in user-facing stuff and Rcpp applications. | I think of the Rf_ part as more due to the fact that C doesn't have a | concept of namespaces so anything in the R API is at the top level | namespace leading to some conflicts. Agreed. But speaking stylistically, for the same reason that we prefer C++ versions of C header files (eg cstdint over stdint.h, cstdio over stdio.h, ...) I am with John on the preference for C++ idioms when given a choice. Dirk -- "Outside of a dog, a book is a man's best friend. Inside of a dog, it is too dark to read." -- Groucho Marx ___ 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
Re: [Rcpp-devel] Optimising 2d convolution
Hi,
Fancy indexing tricks aside, another avenue you might want to look
into is doing the convolution in fourier space ...
I'll stop there so I minimize the damage I'm doing as far as spreading
misinformation in a public forum goes. I'm also interested in
investigating that avenue further, as I wrote a 1d convolution for
something else in Rcpp and want to drive the "fourier trick" via C
code ... I just haven't had the time to do the follow up reading.
Anyway, I "know" it works in the 1d case (that's how the base R
convolve works), I guess there are particulars to sort out in the 2d
case ...
And if you *really* want it to go fast, I guess you can go the CUDA route :-)
http://developer.download.nvidia.com/compute/cuda/2_2/sdk/website/projects/convolutionFFT2D/doc/convolutionFFT2D.pdf
HTH (but I'm sure it didn't :-)
-steve
On Fri, Jan 6, 2012 at 1:43 PM, Hadley Wickham wrote:
> Hi all,
>
> The rcpp-devel list has been very helpful to me so far, so I hope you
> don't mind another question! I'm trying to speed up my 2d convolution
> function:
>
>
>
> library(inline)
> # 2d convolution -
>
> convolve_2d <- cxxfunction(signature(sampleS = "numeric", kernelS =
> "numeric"), plugin = "Rcpp", '
> Rcpp::NumericMatrix sample(sampleS), kernel(kernelS);
> int x_s = sample.nrow(), x_k = kernel.nrow();
> int y_s = sample.ncol(), y_k = kernel.ncol();
>
> Rcpp::NumericMatrix output(x_s + x_k - 1, y_s + y_k - 1);
> for (int row = 0; row < x_s; row++) {
> for (int col = 0; col < y_s; col++) {
> for (int i = 0; i < x_k; i++) {
> for (int j = 0; j < y_k; j++) {
> output(row + i, col + j) += sample(row, col) * kernel(i, j);
> }
> }
> }
> }
> return output;
> ')
>
>
> x <- diag(1000)
> k <- matrix(runif(20* 20), ncol = 20)
> system.time(convolve_2d(x, k))
> # user system elapsed
> # 14.759 0.028 15.524
>
> I have a vague idea that to get better performance I need to get
> closer to bare pointers, and I need to be careful about the order of
> my loops to ensure that I'm working over contiguous chunks of memory
> as often as possible, but otherwise I've basically exhausted my
> C++/Rcpp knowledge. Where should I start looking to improve the
> performance of this function?
>
> The example data basically matches the real problem - x is not usually
> going to be much bigger than 1000 x 1000 and k typically will be much
> smaller. (And hence, based on what I've read, this technique should
> be faster than doing it via a discrete fft)
>
> 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
--
Steve Lianoglou
Graduate Student: Computational Systems Biology
| Memorial Sloan-Kettering Cancer Center
| Weill Medical College of Cornell University
Contact Info: http://cbio.mskcc.org/~lianos/contact
___
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
Re: [Rcpp-devel] Optimising 2d convolution
On Fri, Jan 6, 2012 at 12:43 PM, Hadley Wickham wrote:
> Hi all,
>
> The rcpp-devel list has been very helpful to me so far, so I hope you
> don't mind another question! I'm trying to speed up my 2d convolution
> function:
>
>
>
> library(inline)
> # 2d convolution -
>
> convolve_2d <- cxxfunction(signature(sampleS = "numeric", kernelS =
> "numeric"), plugin = "Rcpp", '
> Rcpp::NumericMatrix sample(sampleS), kernel(kernelS);
> int x_s = sample.nrow(), x_k = kernel.nrow();
> int y_s = sample.ncol(), y_k = kernel.ncol();
>
> Rcpp::NumericMatrix output(x_s + x_k - 1, y_s + y_k - 1);
> for (int row = 0; row < x_s; row++) {
> for (int col = 0; col < y_s; col++) {
> for (int i = 0; i < x_k; i++) {
> for (int j = 0; j < y_k; j++) {
> output(row + i, col + j) += sample(row, col) * kernel(i, j);
> }
> }
> }
> }
> return output;
> ')
>
>
> x <- diag(1000)
> k <- matrix(runif(20* 20), ncol = 20)
> system.time(convolve_2d(x, k))
> # user system elapsed
> # 14.759 0.028 15.524
>
> I have a vague idea that to get better performance I need to get
> closer to bare pointers, and I need to be careful about the order of
> my loops to ensure that I'm working over contiguous chunks of memory
> as often as possible, but otherwise I've basically exhausted my
> C++/Rcpp knowledge. Where should I start looking to improve the
> performance of this function?
>
> The example data basically matches the real problem - x is not usually
> going to be much bigger than 1000 x 1000 and k typically will be much
> smaller. (And hence, based on what I've read, this technique should
> be faster than doing it via a discrete fft)
What are you doing the timing on? On a modest desktop (2.6 GHz Athlon
X4) I get less than a second for this
> library(inline)
> convolve_2d <- cxxfunction(signature(sampleS = "numeric", kernelS =
+ "numeric"), plugin = "Rcpp", '
+Rcpp::NumericMatrix sample(sampleS), kernel(kernelS);
+int x_s = sample.nrow(), x_k = kernel.nrow();
+int y_s = sample.ncol(), y_k = kernel.ncol();
+Rcpp::NumericMatrix output(x_s + x_k - 1, y_s + y_k - 1);
+for (int row = 0; row < x_s; row++) {
+ for (int col = 0; col < y_s; col++) {
+for (int i = 0; i < x_k; i++) {
+ for (int j = 0; j < y_k; j++) {
+output(row + i, col + j) += sample(row, col) * kernel(i, j);
+ }
+}
+ }
+}
+return output;
+ ')
> x <- diag(1000)
> k <- matrix(runif(20* 20), ncol = 20)
> system.time(convolve_2d(x, k))
user system elapsed
0.864 0.000 0.862
___
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
Re: [Rcpp-devel] Idiom for accessing scalars
On Fri, Jan 6, 2012 at 12:39 PM, John Chambers wrote: > At the risk of derailing this thread into the quicksand of programming style > > > Except when required by the specific application, my preference would be to > stay with the Rcpp idiom. Mixing in the older C API seems to risk more > programming error. Of course, each application is different and we have to > make decisions about efficiency (machine) vs possible inefficiency (human). I would argue that special-purpose functions and macros for handling the funky C structures are the best way of handling the funky C structures. Relative to working with C++ classes and methods they are clunky but, short of converting to something like CXXR, we are stuck with the C structures. > The "Rf_" part of the API in particular is ugly and somewhat of an add-on > forced in a few examples by the use of some common names in the macro files. But, as it stands, that is a requirement when using Rcpp. I think of the Rf_ part as more due to the fact that C doesn't have a concept of namespaces so anything in the R API is at the top level namespace leading to some conflicts. > This year in our Stanford grad course, I'm planning to push Rcpp, so in a > few months I may have a different view. :-) > > Cheers, > John > > > On 1/6/12 10:15 AM, Dirk Eddelbuettel wrote: >> >> On 6 January 2012 at 13:00, Steve Lianoglou wrote: >> | Hi, >> | >> | 2012/1/6 Douglas Bates: >> | [snip] >> |> As I mentioned in another thread, I prefer the idiom >> |> >> |> int n2 = ::Rf_asInteger(n); >> |> >> |> because asInteger is part of the R API (in Rcpp it must be called as >> |> Rf_asInteger - the :: is a hint to the compiler that it will be found >> |> in the global namespace). Functions like asInteger and asReal are >> |> used in thousands of places in the base R code and have been optimized >> |> to death as well as being as general as they can possibly be. >> | >> | Cool ... I actually missed that tip, thanks for pointing it out again. >> >> And if one overcomes the "ick" factor of mixing APIs, it saves a >> little: >> >> R> library(inline) >> R> library(microbenchmark) >> R> >> R> fr<- cxxfunction(signature(xs="integer"), plugin="Rcpp", body=' >> + int x = ::Rf_asInteger(xs); >> + return Rcpp::wrap(x); >> + ') >> R> >> R> frcpp<- cxxfunction(signature(xs="integer"), plugin="Rcpp", body=' >> + int x = Rcpp::as(xs); >> + return Rcpp::wrap(x); >> + ') >> R> >> R> microbenchmark(fr(123), frcpp(123)) >> Unit: nanoseconds >> expr min lq median uq max >> 1 fr(123) 841 878.5 908.5 976 15684 >> 2 frcpp(123) 976 999.5 1017.0 1085 6589 >> R> >> >> A good chunk here is "fixed" cost of calling a function, returning a >> value, ... so that the "variable" gain of ::Rf_asInteger() looks pretty >> good. >> >> Dirk, somewhat wondering why we're debating 100ns gains in the context of >> R >> > ___ > 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 ___ 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
Re: [Rcpp-devel] Idiom for accessing scalars
On 6 January 2012 at 10:39, John Chambers wrote: | At the risk of derailing this thread into the quicksand of programming | style Au contraire: thanks for joining in. | Except when required by the specific application, my preference would be | to stay with the Rcpp idiom. Mixing in the older C API seems to risk | more programming error. Of course, each application is different and we | have to make decisions about efficiency (machine) vs possible | inefficiency (human). Couldn't agree more. | The "Rf_" part of the API in particular is ugly and somewhat of an | add-on forced in a few examples by the use of some common names in the | macro files. | | This year in our Stanford grad course, I'm planning to push Rcpp, so in | a few months I may have a different view. :-) Let the bug reports roll in :) Dirk | Cheers, |John | | On 1/6/12 10:15 AM, Dirk Eddelbuettel wrote: | > On 6 January 2012 at 13:00, Steve Lianoglou wrote: | > | Hi, | > | | > | 2012/1/6 Douglas Bates: | > | [snip] | > |> As I mentioned in another thread, I prefer the idiom | > |> | > |> int n2 = ::Rf_asInteger(n); | > |> | > |> because asInteger is part of the R API (in Rcpp it must be called as | > |> Rf_asInteger - the :: is a hint to the compiler that it will be found | > |> in the global namespace). Functions like asInteger and asReal are | > |> used in thousands of places in the base R code and have been optimized | > |> to death as well as being as general as they can possibly be. | > | | > | Cool ... I actually missed that tip, thanks for pointing it out again. | > | > And if one overcomes the "ick" factor of mixing APIs, it saves a | > little: | > | > R> library(inline) | > R> library(microbenchmark) | > R> | > R> fr<- cxxfunction(signature(xs="integer"), plugin="Rcpp", body=' | > +int x = ::Rf_asInteger(xs); | > +return Rcpp::wrap(x); | > + ') | > R> | > R> frcpp<- cxxfunction(signature(xs="integer"), plugin="Rcpp", body=' | > +int x = Rcpp::as(xs); | > +return Rcpp::wrap(x); | > + ') | > R> | > R> microbenchmark(fr(123), frcpp(123)) | > Unit: nanoseconds | > expr minlq median uq max | > 1fr(123) 841 878.5 908.5 976 15684 | > 2 frcpp(123) 976 999.5 1017.0 1085 6589 | > R> | > | > A good chunk here is "fixed" cost of calling a function, returning a | > value, ... so that the "variable" gain of ::Rf_asInteger() looks pretty good. | > | > Dirk, somewhat wondering why we're debating 100ns gains in the context of R | > -- "Outside of a dog, a book is a man's best friend. Inside of a dog, it is too dark to read." -- Groucho Marx ___ 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
[Rcpp-devel] Optimising 2d convolution
Hi all,
The rcpp-devel list has been very helpful to me so far, so I hope you
don't mind another question! I'm trying to speed up my 2d convolution
function:
library(inline)
# 2d convolution -
convolve_2d <- cxxfunction(signature(sampleS = "numeric", kernelS =
"numeric"), plugin = "Rcpp", '
Rcpp::NumericMatrix sample(sampleS), kernel(kernelS);
int x_s = sample.nrow(), x_k = kernel.nrow();
int y_s = sample.ncol(), y_k = kernel.ncol();
Rcpp::NumericMatrix output(x_s + x_k - 1, y_s + y_k - 1);
for (int row = 0; row < x_s; row++) {
for (int col = 0; col < y_s; col++) {
for (int i = 0; i < x_k; i++) {
for (int j = 0; j < y_k; j++) {
output(row + i, col + j) += sample(row, col) * kernel(i, j);
}
}
}
}
return output;
')
x <- diag(1000)
k <- matrix(runif(20* 20), ncol = 20)
system.time(convolve_2d(x, k))
#user system elapsed
# 14.759 0.028 15.524
I have a vague idea that to get better performance I need to get
closer to bare pointers, and I need to be careful about the order of
my loops to ensure that I'm working over contiguous chunks of memory
as often as possible, but otherwise I've basically exhausted my
C++/Rcpp knowledge. Where should I start looking to improve the
performance of this function?
The example data basically matches the real problem - x is not usually
going to be much bigger than 1000 x 1000 and k typically will be much
smaller. (And hence, based on what I've read, this technique should
be faster than doing it via a discrete fft)
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
Re: [Rcpp-devel] Idiom for accessing scalars
At the risk of derailing this thread into the quicksand of programming style Except when required by the specific application, my preference would be to stay with the Rcpp idiom. Mixing in the older C API seems to risk more programming error. Of course, each application is different and we have to make decisions about efficiency (machine) vs possible inefficiency (human). The "Rf_" part of the API in particular is ugly and somewhat of an add-on forced in a few examples by the use of some common names in the macro files. This year in our Stanford grad course, I'm planning to push Rcpp, so in a few months I may have a different view. :-) Cheers, John On 1/6/12 10:15 AM, Dirk Eddelbuettel wrote: On 6 January 2012 at 13:00, Steve Lianoglou wrote: | Hi, | | 2012/1/6 Douglas Bates: | [snip] |> As I mentioned in another thread, I prefer the idiom |> |> int n2 = ::Rf_asInteger(n); |> |> because asInteger is part of the R API (in Rcpp it must be called as |> Rf_asInteger - the :: is a hint to the compiler that it will be found |> in the global namespace). Functions like asInteger and asReal are |> used in thousands of places in the base R code and have been optimized |> to death as well as being as general as they can possibly be. | | Cool ... I actually missed that tip, thanks for pointing it out again. And if one overcomes the "ick" factor of mixing APIs, it saves a little: R> library(inline) R> library(microbenchmark) R> R> fr<- cxxfunction(signature(xs="integer"), plugin="Rcpp", body=' +int x = ::Rf_asInteger(xs); +return Rcpp::wrap(x); + ') R> R> frcpp<- cxxfunction(signature(xs="integer"), plugin="Rcpp", body=' +int x = Rcpp::as(xs); +return Rcpp::wrap(x); + ') R> R> microbenchmark(fr(123), frcpp(123)) Unit: nanoseconds expr minlq median uq max 1fr(123) 841 878.5 908.5 976 15684 2 frcpp(123) 976 999.5 1017.0 1085 6589 R> A good chunk here is "fixed" cost of calling a function, returning a value, ... so that the "variable" gain of ::Rf_asInteger() looks pretty good. Dirk, somewhat wondering why we're debating 100ns gains in the context of R ___ 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
Re: [Rcpp-devel] When does using iterators for subscripting help?
> 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 x_(as >(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 x_(as >(x)) - you're creating a new Map (templated to contain VectorXds) called x_ and initialising with a converted x? * I find it hard to visually parse >(x). Is there a reason you don't write it as x_(as>(x)? Maybe because that would then form the >> operator? * 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
Re: [Rcpp-devel] Idiom for accessing scalars
On 6 January 2012 at 13:00, Steve Lianoglou wrote: | Hi, | | 2012/1/6 Douglas Bates : | [snip] | > As I mentioned in another thread, I prefer the idiom | > | > int n2 = ::Rf_asInteger(n); | > | > because asInteger is part of the R API (in Rcpp it must be called as | > Rf_asInteger - the :: is a hint to the compiler that it will be found | > in the global namespace). Functions like asInteger and asReal are | > used in thousands of places in the base R code and have been optimized | > to death as well as being as general as they can possibly be. | | Cool ... I actually missed that tip, thanks for pointing it out again. And if one overcomes the "ick" factor of mixing APIs , it saves a little: R> library(inline) R> library(microbenchmark) R> R> fr <- cxxfunction(signature(xs="integer"), plugin="Rcpp", body=' +int x = ::Rf_asInteger(xs); +return Rcpp::wrap(x); + ') R> R> frcpp <- cxxfunction(signature(xs="integer"), plugin="Rcpp", body=' +int x = Rcpp::as(xs); +return Rcpp::wrap(x); + ') R> R> microbenchmark(fr(123), frcpp(123)) Unit: nanoseconds expr minlq median uq max 1fr(123) 841 878.5 908.5 976 15684 2 frcpp(123) 976 999.5 1017.0 1085 6589 R> A good chunk here is "fixed" cost of calling a function, returning a value, ... so that the "variable" gain of ::Rf_asInteger() looks pretty good. Dirk, somewhat wondering why we're debating 100ns gains in the context of R -- "Outside of a dog, a book is a man's best friend. Inside of a dog, it is too dark to read." -- Groucho Marx ___ 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
Re: [Rcpp-devel] Idiom for accessing scalars
Hi, 2012/1/6 Douglas Bates : [snip] > As I mentioned in another thread, I prefer the idiom > > int n2 = ::Rf_asInteger(n); > > because asInteger is part of the R API (in Rcpp it must be called as > Rf_asInteger - the :: is a hint to the compiler that it will be found > in the global namespace). Functions like asInteger and asReal are > used in thousands of places in the base R code and have been optimized > to death as well as being as general as they can possibly be. Cool ... I actually missed that tip, thanks for pointing it out again. Cheers, -steve -- Steve Lianoglou Graduate Student: Computational Systems Biology | Memorial Sloan-Kettering Cancer Center | Weill Medical College of Cornell University Contact Info: http://cbio.mskcc.org/~lianos/contact ___ 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
Re: [Rcpp-devel] Idiom for accessing scalars
2012/1/4 Romain François : > Le 04/01/12 14:29, Hadley Wickham a écrit : > >> Hi all, >> >> I'm still just getting my feet wet in Rcpp, so please excuse the >> naivety of my question, but is this the appropriate idiom for treating >> an input as a C++ scalar? >> >> f<- cxxfunction(signature(x = "integer"), plugin = "Rcpp", ' >> Rcpp::IntegerVector x1(x); >> int x2 = x1[0]; >> >> return(Rcpp::NumericVector(x2)); >> ') >> >> i.e. is there any way I can turn x in x2 more succinctly? (This >> function is just for illustrative purposes - I'm trying to do >> something at least a little more complicated) >> >> Thanks! >> >> Hadley > > Yep, > > We have "as" for this, and wrap for the other way around: > > > f<- cxxfunction(signature(x = "integer"), plugin = "Rcpp", ' > int x2 = as(x) ; > > return wrap(x2) ; > ') As I mentioned in another thread, I prefer the idiom int n2 = ::Rf_asInteger(n); because asInteger is part of the R API (in Rcpp it must be called as Rf_asInteger - the :: is a hint to the compiler that it will be found in the global namespace). Functions like asInteger and asReal are used in thousands of places in the base R code and have been optimized to death as well as being as general as they can possibly be. ___ 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
Re: [Rcpp-devel] When does using iterators for subscripting help?
2012/1/4 Hadley Wickham : >> NumericVector:::iterator is actually alias to double*. Here is a trick >> (probably does not work on not gcc compilers): > > Ah, interesting - thanks! 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 x_(as >(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. With this version I get a minor speedup > microbenchmark(operator = count_bin(x, binwidth, origin, nbins = n), +iterator = count_bini(x, binwidth, origin, nbins = n), +eigen= count_eigen(x, binwidth, origin, nbins = n) + ) Unit: milliseconds expr min lq median uq max 1eigen 145.6456 145.7102 145.7536 145.8210 151.0456 2 iterator 153.3059 153.3603 153.4182 153.6056 155.3246 3 operator 156.2418 156.7063 156.7637 156.8982 159.8635 ___ 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
