Using gsl_function using a supplied R function

Question

I am currently attempting to use the gsl_function from GSL library through RcppGSL using a .cpp file and calling it using sourceCpp(). The idea is to perform numerical integration with gsl_integration_qags, also from GSL. My C code invokes a user defined R function (SomeRFunction in my code below) saved into the global environment. The code is:

#include <RcppGSL.h>
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_vector.h>

// [[Rcpp::depends(RcppGSL)]]
// [[Rcpp::export]]
double integfn1(double ub, double lb, double abserr, double Rsq, int pA){ 

  double result, abserror;
  Rcpp::Environment global = Rcpp::Environment::global_env();
  Rcpp::Function f1 = global["SomeRFunction"];

  struct my_f_params { double Rsq; int pA; };
  struct my_f_params params = { Rsq, pA };
  gsl_function F;
  F.function = & f1;
  F.params = &params;
  double lb1 =  lb;
  double ub1 = ub;

  gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000);

  gsl_integration_qags (&F, lb1, ub1, 1e-8, 1e-8, 1000, w, &result, &abserror);

  printf ("result          = % .18f\n", result);
  printf ("estimated error = % .18f\n", abserror);

  gsl_integration_workspace_free (w);

  return result;
}

And the following error comes up:

"cannot convert 'Rcpp::Function* {aka Rcpp::Function_Impl<Rcpp::PreserveStorage>*}' to 'double (*)(double, void*)' "

The problem is in the line where I declare what the function to integrate is (i.e. "F.function = & f1;").

I looked up similar problems but couldn't find anything listed... Any hint will be greatly appreciated!

Many thanks

Answer 1

I created a working example (and timed it against C) in which you can pass an arbitrary user-defined R function to the GSL function QAWF. You should be able to generalize it to other gsl functions as well.

See the example here: https://sites.google.com/site/andrassali/computing/user-supplied-functions-in-rcppgsl

As noted above, the competing C implementation is much-much faster.

Answer 2

Quick comments:

1. You write: My C code invokes a user defined R function (SomeRFunction in my code below) saved into the global environment. So your C code will still slow down a lot at each function evaluation due to the cost of calling R, and the slower speed of R.

2. My RcppDE package (which is also for optmisation) as as example of using external pointers (Rcpp::XPtr) to pass a user-defined C function down to the optimiser. Same flexibility. better speed.

3. Your compile error is exactly at that intersection---you can't "just pass" an Rcpp object to a void pointer. Rcpp::XPtr helps you there, but you too need to know what you are doing so a the working example may be helpful.

4. This is really an Rcpp question but you didn't add the tag, so I just did.

Answer 3

F.function expects a function of this signature double (*) (double x, void * params) so a function taking a double then a void*. You need to help Rcpp::Function if you want this to fly.

This is typical C apis stripping away types. So what you need to give as function is something that understand what you have as params and is able to call the R function, with the right parameters and convert to a double. The simple way to do this is to consider the R function as data, so augment the struct to make it something like this:

struct my_f_params { double Rsq; int pA; Function& rfun };

This way, the function you pass as function can be something like this:

double proxy_fun( double x, void* params ){
  my_f_params* data = reinterpret_cast<my_f_params*>(params) ;
  return as<double>( data->rfun(x, data->Rsq, data->pA ) ) ;
} ;

So you can build your gsl_function something like this:

my_f_params params = { Rsq, pA, f1 };
gsl_function F ; 
F.function = &proxy_fun ;
F.params = &params;

With a bit of C++11 variadic templates and tuples, you could generalize this to anything instead of the (double, int) pair but it is more work.