[Rcpp-devel] NumericVector Double mismatch when indexing an array

[Jeffrey Pollock]

I was interested to see if there was any real speed difference between the
different methods suggested, and it looks like... there isn't...

library(inline)
library(Rcpp)
library(rbenchmark)

fun1 <- cxxfunction(body = '
                RNGScope scope;
                NumericVector rn(5);
                for (int i = 0; i < 5; i++) {
                    rn[i] = as<double>(rnorm(1, 0.0, 1.0));
                }
                return rn;
                ', plugin = "Rcpp")

fun2 <- cxxfunction(body = '
                RNGScope scope;
                NumericVector rn(5);
                for (int i = 0; i < 5; i++) {
                    rn[i] = Rf_rnorm(0.0, 1.0);
                }
                return rn;
                ', plugin = "Rcpp")

fun3 <- cxxfunction(body = '
                RNGScope scope;
                NumericVector rn(5);
                for (int i = 0; i < 5; i++) {
                    rn[i] = norm_rand();
                }
                return rn;
                ', plugin = "Rcpp")

fun4 <- cxxfunction(body = '
                RNGScope scope;
                return rnorm(5, 0.0, 1.0);
                ', plugin = "Rcpp")

set.seed(1)
print(fun1())

set.seed(1)
print(fun2())

set.seed(1)
print(fun3())

set.seed(1)
print(fun4())

benchmark(fun1(), fun2(), fun3(), fun4(), order = "relative", replications
= 1e6L)

output;

> set.seed(1)
> print(fun1())
[1] -0.6264538  0.1836433 -0.8356286  1.5952808  0.3295078

> set.seed(1)
> print(fun2())
[1] -0.6264538  0.1836433 -0.8356286  1.5952808  0.3295078

> set.seed(1)
> print(fun3())
[1] -0.6264538  0.1836433 -0.8356286  1.5952808  0.3295078

> set.seed(1)
> print(fun4())
[1] -0.6264538  0.1836433 -0.8356286  1.5952808  0.3295078

> benchmark(fun1(), fun2(), fun3(), fun4(), order = "relative",
replications = 1e6L)
    test replications elapsed relative user.self sys.self user.child
sys.child
3 fun3()      1000000    7.15    1.000      7.01     0.00         NA
NA
4 fun4()      1000000    7.33    1.025      7.27     0.01         NA
NA
2 fun2()      1000000    7.53    1.053      7.41     0.00         NA
NA
1 fun1()      1000000    7.99    1.117      7.91     0.00         NA
NA


On Tue, Sep 25, 2012 at 8:19 PM, Goldfeld, Keith
<Keith.Goldfeld at nyumc.org>wrote:

> Yes - I agree that iterating is not the best way to go, but I am using
> this for an agent based simulation where iteration is really my only
> option.  I won't be generating more than a few thousand records, so the
> time shouldn't be much of a factor.  Of course, in R it is painfully slow -
> so that is why I am going the Rcpp route.
>
> -----Original Message-----
> From: dmbates at gmail.com [mailto:dmbates at gmail.com] On Behalf Of Douglas
> Bates
> Sent: Tuesday, September 25, 2012 3:13 PM
> To: Rodney Sparapani
> Cc: Goldfeld, Keith; rcpp-devel at lists.r-forge.r-project.org
> Subject: Re: [Rcpp-devel] NumericVector Double mismatch when indexing an
> array
>
> On Tue, Sep 25, 2012 at 1:53 PM, Rodney Sparapani <rsparapa at mcw.edu>
> wrote:
> > On 09/25/2012 01:37 PM, Goldfeld, Keith wrote:
> >
> >>>  code<- 'Rcpp::RNGScope scope;
> >>
> >>
> >> +                    Rcpp::NumericVector rn(5);
> >>
> >> +                    for (int i=0; i<  5; i++) {
> >>
> >> +                        rn(i) = rnorm(1,0,1);
> >>
> >> +                    }
> >>
> >> +                    return Rcpp::wrap(rn);
> >>
> >> +  '
> >>
> >>>
> >>
> >>>  fun<- cxxfunction(body=code, plugin="Rcpp")
> >
> >
> > You just need an explicit type conversion like so:
> >
> > funk2 <- cxxfunction(body='RNGScope scope;
> >           NumericVector rn(5);
> >           for (int i=0; i < 5; i++) rn[i] = as<double>(rnorm(1,0,1));
> >           return wrap(rn);',
> >           includes="using namespace Rcpp;\n", plugin="Rcpp")
> >
> > funk2()
>
> That works but you wouldn't want to use it as a general model because you
> are creating vectors of length 1 then dereferencing the first element in
> each of these vectors and copying it to another vector.
> Obviously with n=5 the difference in execution time will be minimal but
> with n=5 million it won't be.  Using the Rcpp sugar function rnorm will be
> the easiest approach unless, like me, you get a little queasy when using
> Rcpp sugar functions just because it is so hard to pick them apart and see
> what is actually happening.
>
> I would probably end up going back to the R API and calling norm_rand or
> Rf_rnorm.  That latter returns a double from two double arguments, mu and
> sigma.
> _______________________________________________
> Rcpp-devel mailing list
> Rcpp-devel at lists.r-forge.r-project.org
> https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20120925/1769fa62/attachment.html>