Computing‎ > ‎

User supplied functions in RcppGSL

User supplied functions and RcppGSL

User supplied functions and RcppGSL

In this short demo, we show how to use RcppGSL together with user-defined R functions. We also show through time benchmarks that this might not necessarily be a useful way to go and that coding functions directly in C can result in much quicker execution.

To take an example, let us assume one is interested in calculating accurately for some . If is large, then this integral is highly oscillatory and traditional quadrature will be inaccurate and time-consuming. Luckily, the GNU Scientific Library has a routine that suits just this case, it is called QAWF (ported from QuadPack). The idea behind this routine is to approximate locally the function using polynomials, which can be exactly integrated against the sine function (note that in traditional quadrature the whole function would be approximated, which is much more difficult for large .

RcppGSL is a great way to interact with GSL package in C, using Rcpp objects. The “problem” with using GSL through RcppGSL is that GSL requires functions defined in a very specific C format. We can use however the magic behind Rcpp to actually do a proper implementation. In the following box we imlement the passing of an arbitrary function to GSL as well as specialized example using a C coded function:

// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::depends(RcppGSL)]]

#include <RcppGSL.h>
#include <gsl/gsl_integration.h>

// set up struct that contains an Rcpp Function, this will be used in the function provided to gsl
struct my_f_params {Rcpp::Function G;};

// R callback function
double myObjB(double x, void *p) {
  struct my_f_params *params = (struct my_f_params *)p;
  Rcpp::Function G = (params->G);
  Rcpp::NumericVector y = G(x);
  double res = y[0];
  return(res);
}

// specialized C implementation
double myObjA(double x, void *p) {
  double res = exp(-x);
  return(res);
}

// [[Rcpp::export]] 
Rcpp::NumericVector qawf_general(Rcpp::Function f, Rcpp::NumericVector omega, bool rFun = true) {
  // setup output vector
  Rcpp::NumericVector resVector = omega;

  // define integration workspace
  gsl_integration_workspace *work = gsl_integration_workspace_alloc(1e5);
  gsl_integration_workspace *cyclew = gsl_integration_workspace_alloc(1e5);
  double currOmega = omega[0];
  gsl_integration_qawo_table *table = gsl_integration_qawo_table_alloc(currOmega, 1, GSL_INTEG_SINE, 1e1);  

  // declare output varaibles
  double result = 0.0;
  double abserr = 0.0;

  gsl_function F;
  struct my_f_params params = {f};

  // if R function is to be used, then point gsl to that, otherwise use C implementation
  if (rFun) {
    F.function = &myObjB;
  } else {
    F.function = &myObjA;    
  }
  F.params = &params;

  for (int nn=0;nn<omega.size();nn++) {
    currOmega = omega[nn];
    int set = gsl_integration_qawo_table_set(table,omega[nn], 1, GSL_INTEG_SINE);
    int res = gsl_integration_qawf(&F, 0.0, 1e-8, 1e5, work, cyclew, table, &result, &abserr);
    resVector[nn] = result;
  }

  gsl_integration_qawo_table_free(table);
  gsl_integration_workspace_free(cyclew);
  gsl_integration_workspace_free(work);

  return resVector;             // return vector  
}

To show how this code works, we can run the following example:

# Demonstrate R callback speed using gsl function qasw ----

# compile and source the C++ function
Rcpp::sourceCpp("qawf_general.cpp")

# start timer
start.time <- proc.time()[3]

# demonstrate equivalence of R and C results
R.result <- qawf_general(function(x) exp(-x), seq(1, 100), TRUE)
R.time <- proc.time()[3]

C.result <- qawf_general(function(x) exp(-x), seq(1, 100), FALSE)
C.time <- proc.time()[3]

print(paste("Maximum relative difference between R and C result is:", max(abs(R.result/C.result - 
    1))))
## [1] "Maximum relative difference between R and C result is: 0"

# now calculate using non-oscillatory quadrature, we need many subdivisions
quad.result <- sapply(seq(1, 100), function(y) integrate(function(x) exp(-x) * 
    sin(y * x), 0, Inf, subdivisions = 1000, abs.tol = 1e-08)$value)
quad.time <- proc.time()[3]

print(paste("Maximum relative difference between qawf and direct quadrature is:", 
    max(abs(quad.result/C.result - 1))))
## [1] "Maximum relative difference between qawf and direct quadrature is: 3.7802436905654e-05"

# now compare timings
print(paste("R function based qawf took:", R.time - start.time, "seconds"))
## [1] "R function based qawf took: 0.215 seconds"
print(paste("C function based qawf took:", C.time - R.time, "seconds"))
## [1] "C function based qawf took: 0.00600000000000023 seconds"
print(paste("Direct (non oscillatory) quadrature took:", quad.time - C.time, 
    "seconds"))
## [1] "Direct (non oscillatory) quadrature took: 0.128 seconds"

So we see that indeed all three approaches work in this case. However, the direct quadrature would not work very well with higher \( \omega \) values and we already see it to be less accurate.

Finally, we note that although passing R functions to GSL does work this way, it is painfully slow compared to passing directly a C function. Therefore it should only be done for very complex functions (that are potentially expensive to calculate anyway), that cannot be reasonably be ported to C. However, if a C port is feasible, it is usually much-much faster.


Comments