Rcpp Version 0.12.14
timeRNGs.R
Go to the documentation of this file.
1 
2 suppressMessages(library(Rcpp))
3 suppressMessages(library(inline))
4 suppressMessages(library(rbenchmark))
5 
6 ## NOTE: This is the old way to compile Rcpp code inline.
7 ## The code here has left as a historical artifact and tribute to the old way.
8 ## Please use the code under the "new" inline compilation section.
9 
10 rcppGamma_old <- cxxfunction(signature(xs="numeric"), plugin="Rcpp", body='
11  NumericVector x(xs);
12  int n = x.size();
13 
14  // Initialize Random number generator
15  RNGScope scope;
16 
17  const double y = 1.234;
18  for (int i=0; i<n; i++) {
19  x[i] = ::Rf_rgamma(3.0, 1.0/(y*y+4));
20  }
21 
22  // Return to R
23  return x;
24 ')
25 
26 
27 gslGamma_old <- cxxfunction(signature(xs="numeric"), plugin="RcppGSL",
28  include='#include <gsl/gsl_rng.h>
29  #include <gsl/gsl_randist.h>',
30  body='
31  NumericVector x(xs);
32  int n = x.size();
33 
34  gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937);
35  const double y = 1.234;
36  for (int i=0; i<n; i++) {
37  x[i] = gsl_ran_gamma(r,3.0,1.0/(y*y+4));
38  }
39  gsl_rng_free(r);
40 
41  // Return to R
42  return x;
43 ')
44 
45 
46 rcppNormal_old <- cxxfunction(signature(xs="numeric"), plugin="Rcpp", body='
47  NumericVector x(xs);
48  int n = x.size();
49 
50  // Initialize Random number generator
51  RNGScope scope;
52 
53  const double y = 1.234;
54  for (int i=0; i<n; i++) {
55  x[i] = ::Rf_rnorm(1.0/(y+1),1.0/sqrt(2*y+2));
56  }
57 
58  // Return to R
59  return x;
60 ')
61 
62 
63 gslNormal_old <- cxxfunction(signature(xs="numeric"), plugin="RcppGSL",
64  include='#include <gsl/gsl_rng.h>
65  #include <gsl/gsl_randist.h>',
66  body='
67  NumericVector x(xs);
68  int n = x.size();
69 
70  gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937);
71  const double y = 1.234;
72  for (int i=0; i<n; i++) {
73  x[i] = 1.0/(y+1)+gsl_ran_gaussian(r,1.0/sqrt(2*y+2));
74  }
75  gsl_rng_free(r);
76 
77  // Return to R
78  return x;
79 ')
80 
81 
82 ## NOTE: Within this section, the new way to compile Rcpp code inline has been
83 ## written. Please use the code next as a template for your own project.
84 
85 cppFunction('
86 NumericVector rcppGamma(NumericVector x){
87  int n = x.size();
88 
89  const double y = 1.234;
90  for (int i=0; i<n; i++) {
91  x[i] = R::rgamma(3.0, 1.0/(y*y+4));
92  }
93 
94  // Return to R
95  return x;
96 }')
97 
98 ## This approach is a bit sloppy. Generally, you will want to use
99 ## sourceCpp() if there are additional includes that are required.
100 cppFunction('
101 NumericVector gslGamma(NumericVector x){
102  int n = x.size();
103 
104  gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937);
105  const double y = 1.234;
106  for (int i=0; i<n; i++) {
107  x[i] = gsl_ran_gamma(r,3.0,1.0/(y*y+4));
108  }
109  gsl_rng_free(r);
110 
111  // Return to R
112  return x;
113 }', includes = '#include <gsl/gsl_rng.h>
114  #include <gsl/gsl_randist.h>',
115  depends = "RcppGSL")
116 
117 
118 cppFunction('
119 NumericVector rcppNormal(NumericVector x){
120  int n = x.size();
121 
122  const double y = 1.234;
123  for (int i=0; i<n; i++) {
124  x[i] = R::rnorm(1.0/(y+1),1.0/sqrt(2*y+2));
125  }
126 
127  // Return to R
128  return x;
129 }')
130 
131 
132 ## Here we demonstrate the use of sourceCpp() to show the continuity
133 ## of the code artifact.
134 
135 sourceCpp(code = '
136 #include <RcppGSL.h>
137 #include <gsl/gsl_rng.h>
138 #include <gsl/gsl_randist.h>
139 
140 using namespace Rcpp;
141 
142 // [[Rcpp::depends("RcppGSL")]]
143 
144 // [[Rcpp::export]]
145 NumericVector gslNormal(NumericVector x){
146  int n = x.size();
147 
148  gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937);
149  const double y = 1.234;
150  for (int i=0; i<n; i++) {
151  x[i] = 1.0/(y+1)+gsl_ran_gaussian(r,1.0/sqrt(2*y+2));
152  }
153  gsl_rng_free(r);
154 
155  // Return to R
156  return x;
157 }')
158 
159 x <- rep(NA, 1e6)
160 res <- benchmark(rcppGamma(x),
161  gslGamma(x),
162  rcppNormal(x),
163  gslNormal(x),
164  columns=c("test", "replications", "elapsed", "relative", "user.self", "sys.self"),
165  order="relative",
166  replications=20)
167 print(res)
168 
169