using C function from other package in Rcpp

Question

I'm trying to call a C routine from the cubature package in a c++ function to perform multidimensional integration.

The basic R example I'm trying to reproduce is

library(cubature)
integrand <- function(x) sin(x)
adaptIntegrate(integrand, 0, pi)

I could just call this R function from Rcpp following this recipe from the gallery, but there would be some performance penalty in switching back and forth from c/c++ to R. It seems more sensible to directly call the C function from C++.

The C routine adapt_integrate is exported from cubature with

 // R_RegisterCCallable("cubature", "adapt_integrate", (DL_FUNC) adapt_integrate);

I don't understand how to call it from c++, however. Here's my lame attempt,

sourceCpp(code = '
#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
double integrand(double x){
 return(sin(x));
}

// [[Rcpp::depends(cubature)]]
// [[Rcpp::export]]
Rcpp::List integratecpp(double llim, double ulim)
{
  Rcpp::Function p_cubature = R_GetCCallable("cubature", "adapt_integrate");

  Rcpp::List result = p_cubature(integrand, llim, ulim);
  return(result);
}
'
)

integratecpp(0, pi)

This fails to compile; clearly I'm doing something very silly and missing some important steps to convert the output of R_GetCCallable into an Rcpp::Function (or call it directly?). I've read several related posts dealing with function pointers, but haven't seen an example using an external C function.

score 11 · Accepted Answer · answered Dec 09 '13 at 16:51

11

Unfortunately cubature does not ship the headers in inst/include, so you have to borrow that from them and do something like this in your code:

typedef void (*integrand) (unsigned ndim, const double *x, void *,
           unsigned fdim, double *fval);

int adapt_integrate(
    unsigned fdim, integrand f, void *fdata,
    unsigned dim, const double *xmin, const double *xmax, 
    unsigned maxEval, double reqAbsError, double reqRelError, 
    double *val, double *err)
{
    typedef int (*Fun)(unsigned,integrand,void*,unsigned,
        const double*,const double*, unsigned, double, double, double*, double*) ;
    Fun fun = (Fun) R_GetCCallable( "cubature", "adapt_integrate" ) ;           
    return fun(fdim,f,fdata,dim,xmin,xmax,maxEval,reqAbsError, reqRelError,val,err); 
}

It might be a good idea to negociate with the maintainer of cubature that he ships declarations in inst/include so that you'd only have to use LinkingTo.

answered Dec 09 '13 at 16:51

Romain Francois

17,432
3
51
77

Many thanks for assembling these missing pieces. I will have to rethink the problem, unfortunately, because from what I see `adapt_integrate` will not easily accept an integrand that I defined using armadillo's data structures. For completeness, would you be able to add a minimal example of use? – baptiste Dec 09 '13 at 17:16
What it gives you is access to the function pointer that `cubature` registered. I don't know what you are then supposed to do with the C function... – Romain Francois Dec 09 '13 at 17:20
indeed, [considering an example of use](http://ab-initio.mit.edu/wiki/index.php/Cubature#Example) I see trouble ahead: `adapt_integrate_v` expects pointers to objects such as `*fdata`, the integrand expects pointers such as `*fval`, while the arguments I really want to pass are e.g. `arma::colvec` objects. I don't think I'll be able to make a bridge between the two. I may have to stick with R-level interface, or implement my own 2D quadrature in c++. – baptiste Dec 09 '13 at 17:37
Sure. I was only answering how to call this function. If you want to call it with some other data type, you have to either do some data conversion work or negociate another interface with the maintainer of `cubature`. – Romain Francois Dec 09 '13 at 18:14

Dirk Eddelbuettel · Answer 2 · 2013-12-09T19:48:51.440

Didn't see this question earlier, and it looks like @Romain addressed it.

For completeness, a working example of how to do this when all parties play along is provided by the xts and RcppXts packages. In xts, we do this (for about ten functions) in the (source) file inst/include/xtsAPI.h:

SEXP attribute_hidden xtsLag(SEXP x, SEXP k, SEXP pad) {     
    static SEXP(*fun)(SEXP,SEXP,SEXP) = NULL;         
    if (fun == NULL)                                  
        fun = (SEXP(*)(SEXP,SEXP,SEXP)) R_GetCCallable("xts","lagXts");   
    return fun(x, k, pad);                               
}

along with the usual business of R_registerRoutines and R_RegisterCCallable.

In RcppXts this is picked up (in an Rcpp Module) as

function("xtsLag", 
         &xtsLag,    
         List::create(Named("x"), Named("k"), Named("pad")),   
         "Extract the coredata from xts object");

which works pretty well. Someone reprimanded me to write the xts side more compactly (as the if NULL is spurious) which I will get to ... eventually.

thanks for the pointers, unfortunately this is a case where working out the proper glue between the two pieces of code would take me way more time and suffering than actually doing the calculations with the slower R implementation. I would otherwise consider linking directly to the original cubature code rather than through the R package, which is using an older version. — baptiste, Dec 10 '13 at 10:56

score 2 · Answer 3 · answered Dec 11 '16 at 19:54

2

This question is three years old now but I want to point out that multidimensional integration with Rcpp may be easier now that the RcppNumerical library is available: https://github.com/yixuan/RcppNumerical

The routines for computing integrals are based on Thomas Hahn's Cuba package and are also available in the R2Cuba library on CRAN, so if you can accept using the Cuba routines over the function adaptIntegrate() from Cubature, this package may be of interest.

answered Dec 11 '16 at 19:54

user31313

131
5

thanks for the pointer. However it's probably best to keep it as a comment, since the question wasn't specifically about integration, but rather about accessing a C function from another package. I do appreciate the link though, I ended up using my own copy of cubature, but this looks like a worthy alternative (particularly if one is already using Eigen, but I'm used to Armadillo). One advantage of cubature though is the ability to deal with vector-valued integrands. – baptiste Dec 11 '16 at 21:22

using C function from other package in Rcpp

3 Answers3

Linked

Related