Modernizing RcppGSL examples

Part of #56's push to switch examples away from the cxxfunction to Rcpp attributes.

Author: coatless

inst/examples/RcppGibbs/RcppGibbs.R

@@ -83,8 +83,12 @@ RCgibbs <- cmpfun(Rgibbs)
 ## mat <- Rgibbs(10000,10); dim(mat)
 ## would give: [1] 10000     2
 
-
 ## Now for the Rcpp version -- Notice how easy it is to code up!
+
+## NOTE: This is the old way to compile Rcpp code inline.
+## The code here has left as a historical artifact and tribute to the old way.
+## Please use the code under the "new" inline compilation section.
+
 gibbscode <- '
 
   using namespace Rcpp;   // inline does that for us already
@@ -114,7 +118,7 @@ gibbscode <- '
 '
 
 # Compile and Load
-RcppGibbs <- cxxfunction(signature(n="int", thin = "int"),
+RcppGibbs_old <- cxxfunction(signature(n="int", thin = "int"),
                          gibbscode, plugin="Rcpp")
 
 
@@ -146,7 +150,7 @@ gslgibbscode <- '
 '
 
 ## Compile and Load
-GSLGibbs <- cxxfunction(signature(ns="int", thns = "int"),
+GSLGibbs_old <- cxxfunction(signature(ns="int", thns = "int"),
                         body=gslgibbscode, includes=gslgibbsincl,
                         plugin="RcppGSL")
 
@@ -158,6 +162,69 @@ GSLGibbs <- cxxfunction(signature(ns="int", thns = "int"),
 #                      libargs="-lgsl -lgslcblas")
 
 
+## NOTE: Within this section, the new way to compile Rcpp code inline has been
+## written. Please use the code next as a template for your own project.
+
+## Use of the cppFunction() gives the ability to immediately compile embed C++
+## without having to worry about header specification or Rcpp attributes.
+
+cppFunction('
+NumericMatrix RcppGibbs(int N, int thn){
+    // Note: n and thin are SEXPs which the Rcpp automatically converts to ints
+
+    // Setup storage matrix
+    NumericMatrix mat(N, 2);
+    
+    // The rest of the code follows the R version
+    double x = 0, y = 0;
+    
+    for (int i = 0; i < N; i++) {
+        for (int j = 0; j < thn; j++) {
+            x = R::rgamma(3.0,1.0/(y*y+4));
+            y = R::rnorm(1.0/(x+1),1.0/sqrt(2*x+2));
+        }
+        mat(i,0) = x;
+        mat(i,1) = y;
+    }
+    
+    return mat;             // Return to R
+}')
+
+
+## Use of the sourceCpp() is preferred for users who wish to source external
+## files or specify their headers and Rcpp attributes within their code.
+## Code here is able to easily be extracted and placed into its own C++ file. 
+
+## Compile and Load
+sourceCpp(code="
+#include <RcppGSL.h>
+#include <gsl/gsl_rng.h>
+#include <gsl/gsl_randist.h>
+
+using namespace Rcpp;  // just to be explicit
+
+// [[Rcpp::depends(RcppGSL)]]
+
+// [[Rcpp::export]]
+NumericMatrix GSLGibbs(int N, int thin){
+    gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937);
+    double x = 0, y = 0;
+    NumericMatrix mat(N, 2);
+    for (int i = 0; i < N; i++) {
+        for (int j = 0; j < thin; j++) {
+            x = gsl_ran_gamma(r,3.0,1.0/(y*y+4));
+            y = 1.0/(x+1)+gsl_ran_gaussian(r,1.0/sqrt(2*x+2));
+        }
+        mat(i,0) = x;
+        mat(i,1) = y;
+    }
+    gsl_rng_free(r);
+    
+    return mat;           // Return to R
+}")
+
+
+
 ## Now for some tests
 ## You can try other values if you like
 ## Note that the total number of interations are N*thin!

inst/examples/RcppGibbs/timeRNGs.R

@@ -3,7 +3,11 @@ suppressMessages(library(Rcpp))
 suppressMessages(library(inline))
 suppressMessages(library(rbenchmark))
 
-rcppGamma <- cxxfunction(signature(xs="numeric"), plugin="Rcpp", body='
+## NOTE: This is the old way to compile Rcpp code inline.
+## The code here has left as a historical artifact and tribute to the old way.
+## Please use the code under the "new" inline compilation section.
+
+rcppGamma_old <- cxxfunction(signature(xs="numeric"), plugin="Rcpp", body='
   NumericVector x(xs);
   int n   = x.size();
 
@@ -20,7 +24,7 @@ rcppGamma <- cxxfunction(signature(xs="numeric"), plugin="Rcpp", body='
 ')
 
 
-gslGamma <- cxxfunction(signature(xs="numeric"), plugin="RcppGSL",
+gslGamma_old <- cxxfunction(signature(xs="numeric"), plugin="RcppGSL",
                         include='#include <gsl/gsl_rng.h>
                                  #include <gsl/gsl_randist.h>',
                         body='
@@ -39,7 +43,7 @@ gslGamma <- cxxfunction(signature(xs="numeric"), plugin="RcppGSL",
 ')
 
 
-rcppNormal <- cxxfunction(signature(xs="numeric"), plugin="Rcpp", body='
+rcppNormal_old <- cxxfunction(signature(xs="numeric"), plugin="Rcpp", body='
   NumericVector x(xs);
   int n   = x.size();
 
@@ -56,7 +60,7 @@ rcppNormal <- cxxfunction(signature(xs="numeric"), plugin="Rcpp", body='
 ')
 
 
-gslNormal <- cxxfunction(signature(xs="numeric"), plugin="RcppGSL",
+gslNormal_old <- cxxfunction(signature(xs="numeric"), plugin="RcppGSL",
                         include='#include <gsl/gsl_rng.h>
                                  #include <gsl/gsl_randist.h>',
                         body='
@@ -75,6 +79,83 @@ gslNormal <- cxxfunction(signature(xs="numeric"), plugin="RcppGSL",
 ')
 
 
+## NOTE: Within this section, the new way to compile Rcpp code inline has been
+## written. Please use the code next as a template for your own project.
+
+cppFunction('
+NumericVector rcppGamma(NumericVector x){
+    int n   = x.size();
+    
+    const double y = 1.234;
+    for (int i=0; i<n; i++) {
+        x[i] = R::rgamma(3.0, 1.0/(y*y+4));
+    }
+    
+    // Return to R
+    return x;
+}')
+
+## This approach is a bit sloppy. Generally, you will want to use 
+## sourceCpp() if there are additional includes that are required.
+cppFunction('
+NumericVector gslGamma(NumericVector x){
+    int n   = x.size();
+    
+    gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937);
+    const double y = 1.234;
+    for (int i=0; i<n; i++) {
+        x[i] = gsl_ran_gamma(r,3.0,1.0/(y*y+4));
+    }
+    gsl_rng_free(r);
+    
+    // Return to R
+    return x;
+}', includes = '#include <gsl/gsl_rng.h>
+                #include <gsl/gsl_randist.h>',
+    depends = "RcppGSL")
+
+
+cppFunction('
+NumericVector rcppNormal(NumericVector x){
+    int n   = x.size();
+    
+    const double y = 1.234;
+    for (int i=0; i<n; i++) {
+        x[i] = R::rnorm(1.0/(y+1),1.0/sqrt(2*y+2));
+    }
+    
+    // Return to R
+    return x;
+}')
+
+
+## Here we demonstrate the use of sourceCpp() to show the continuity 
+## of the code artifact.
+
+sourceCpp(code = '
+#include <RcppGSL.h>
+#include <gsl/gsl_rng.h>
+#include <gsl/gsl_randist.h>
+
+using namespace Rcpp;
+
+// [[Rcpp::depends("RcppGSL")]]
+
+// [[Rcpp::export]]
+NumericVector gslNormal(NumericVector x){
+    int n   = x.size();
+    
+    gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937);
+    const double y = 1.234;
+    for (int i=0; i<n; i++) {
+        x[i] = 1.0/(y+1)+gsl_ran_gaussian(r,1.0/sqrt(2*y+2));
+    }
+    gsl_rng_free(r);
+    
+    // Return to R
+    return x;
+}')
+
 x <- rep(NA, 1e6)
 res <- benchmark(rcppGamma(x),
                  gslGamma(x),