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

[Christian Gunning]

> Indeed, changing the sample size from 5 to 500 results in;

For completeness, I've included an iterator example, tested over a
range of sample sizes, and plotted the results with the code below.
>From where I sit, it looks like the relative timings flatten out (more
or less) at 1e4.

I'm not sure the norm_rand is really an apples-to-apples comparison,
since this function call may be incrementally faster due to lack of
arguments.

Honestly, this is a testament to how well the sugar code works, if you
accept that it does in fact work :)
-Christian

library(inline)
library(Rcpp)
library(rbenchmark)
library(plyr)
library(lattice)

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

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

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

fun4 <- cxxfunction(signature(N="integer"), body = '
                int NN = as<int>(N);
                RNGScope scope;
                return rnorm(NN, 0.0, 1.0);
                ', plugin = "Rcpp")


fun5 <- cxxfunction(signature(N="integer"), body = '
                int NN = as<int>(N);
                RNGScope scope;
                NumericVector rn(NN);
                std::generate( rn.begin(), rn.end(), RandomNumber);
                return rn; ',
                includes='
                    inline double RandomNumber () { return
Rf_rnorm(0.0, 1.0); }',
                plugin = "Rcpp")


test.N = 4
set.seed(1)
print(fun1(test.N))

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

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

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

set.seed(1)
print(fun5(test.N))

results <- ldply(c(10^(1:5)), function(N) {
  ret <- benchmark(fun1(N), fun2(N), fun3(N), fun4(N), fun5(N),
            order = "relative", replications = 1e3L)
  ret$N <- N
  ret
}, .progress='text')

results$test = factor(results$test, labels=c("Loop, as<>", "Loop,
Rf_rnorm", "Loop, norm_rand", "Sugar", "Iterator"))

plot(
    xyplot( relative ~ N, groups=test, auto.key=list(space='right'),
            type=c('l','g', 'p'), data=results,
            scales=list(x=list(log=T)))
)

-- 
A man, a plan, a cat, a ham, a yak, a yam, a hat, a canal – Panama!