CODEWITHSUNDEEP
c++crrcpp
baptiste
baptiste

Reputation: 77096

using C function from other package in Rcpp

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.

Upvotes: 19

Views: 3344

Answers (3)

user31313
user31313

Reputation: 131

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.

Upvotes: 2

Dirk is no longer here
Dirk is no longer here

Reputation: 368181

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.

Upvotes: 7

Romain Francois
Romain Francois

Reputation: 17632

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.

Upvotes: 11

Related Questions