Skip to contents

Using integratecpp

Source: vignettes/integratecpp.Rmd
integratecpp.Rmd

The following demonstrates on a simplified problem how using integratecpp simplifies numerical integration in Rcpp compared to using R’s C API. 1 Suppose, we want need to evaluate the upper incomplete Gamma function \[ \Gamma{(\alpha, x)} = \int_{x}^{\infty}{ \exp{\{ - u \}} u^{\alpha-1} {\mathrm{d}u} } , \quad x > 0 , \] for values \(\alpha \in \mathbb{R}\), \(x \geq 0\).23

An implementation with integratecpp

As integratecpp requires C++11, we use the Rcpp plugin cpp11 and link to the includes directory of integratecpp using the depends attribute integratecpp:

// C++11

// [[Rcpp::plugins(cpp11)]]
// [[Rcpp::depends(integratecpp)]]

Next, we include the STL header <cmath> for using std::pow and we include the header <integratecpp.h>:

#include <cmath>
#include <integratecpp.h>

Finally, we implement our function using the integratecpp interface to Rdqagi:

// [[Rcpp::export]]
double gamma_inc(const double alpha, const double x) {
    // NOTE: integrand implemented as capturing lambda
    const auto fn = [alpha](const double u) {
        return std::exp(-u) * std::pow(u, alpha - 1.);
    };

    // NOTE: throws if integration error occurs
    const auto result = integratecpp::integrate(
        fn, x, std::numeric_limits<double>::infinity());

    return result.value;
}
# R

## demo with sane input
gamma_inc(-0.5, 1)
#> [1] 0.1781477

## demo with not so sane input
tryCatch(gamma_inc(0, 0), error = function(cond) print(cond))
#> <integratecpp::bad_integrand_error in eval(expr, envir, enclos): extremely bad integrand behaviour>
tryCatch(gamma_inc(0, -1), error = function(cond) print(cond))
#> <integratecpp::integration_runtime_error in eval(expr, envir, enclos): non-finite function value>

Note, that in this basic form, using integratecpp’s wrapper requires only the following two things:

Optionally, integration parameters can be changed, exceptions can be captured and handled, and the integration error estimate can be further examined.

An implementation without integratecpp

The following implementation uses R’s C API instead of the integratecpp wrapper and hints how the wrapper works internally. First, we need to include the header for R’s C entry-points for the integration routines Rdqags and Rdqagi:

// C++11

#include <R_ext/Applic.h>

Then, we require some STL headers for handling data:

#include <algorithm> // std::transform
#include <utility>   // std::make_pair, std::pair
#include <vector>    // std::vector

Also, we need more STL headers and the Rcpp header for exception and error handling:

#include <exception>  // std::exception_ptr, std::make_exception_ptr, 
                      // std::rethrow_exception, std::current_exception
#include <stdexcept>  // std::exception, std::runtime_error
#include <Rcpp.h>     // Rcpp::stop

Finally, we can implement the upper incomplete Gamma function:

// [[Rcpp::export]]
double gamma_inc_direct(const double alpha, const double x) {
    // NOTE: integrand implemented as capturing lambda
    const auto fn = [alpha](const double u) {
        return std::exp(-u) * std::pow(u, alpha - 1.);
    };

    // create bounds parameters for `Rdqagi`
    double bound = x;
    int inf = 1;

    // initialize output variables for integration results
    double result;
    double abserr;
    int last;
    int neval;

    // initialize configuration parameters
    int limit = 100;
    double epsrel = 0.0001220703125;
    double epsabs = epsrel;
    int lenw = 4 * limit;

    // initialize working array
    auto iwork = std::vector<int>(limit);
    auto work = std::vector<double>(lenw);

    // initialize variable for error code
    int ier = 0;

    // NOTE: `Rdqagi` requires a function pointer with signature
    // `void(*)(double *, int, void *)` and a void pointer
    // `void *` passed to the callback as the last argument
    const auto fn_callback = [](double *x, int n, void *ex) {
        // cast `void *` to original type and create references
        // to the `Callable` and the exception pointer
        using ex_t = std::pair<decltype(fn), std::exception_ptr>;
        auto &fn_integrand = (*static_cast<ex_t *>(ex)).first;
        auto &e_ptr = (*static_cast<ex_t *>(ex)).second;
        
        // evaluate `Callable` at *in-out* pointer-to-array `x`
        // and catch possible exceptions in exception pointer
        try {
            std::transform(&x[0], &x[n], &x[0], fn_integrand);
        } catch (const std::exception& e) {
            e_ptr = std::current_exception();
        }
        
        // store a runtime exception in the exception pointer 
        // if any results are infinite 
        if (!std::all_of(&x[0], &x[n], 
                         [](const double x) { 
                             return std::isfinite(x); 
                         })) {
            e_ptr = std::make_exception_ptr(
                std::runtime_error("non-finite function value"));
          }
    };

    // create pair with `Callable` and exception pointer 
    // (`nullptr` on initialization) as external data for 
    // callback function
    auto ex = std::make_pair(fn, std::exception_ptr());
    
    // call C-method `Rdqagi`
    Rdqagi(fn_callback, &ex, &bound, &inf, &epsabs, &epsrel, &result,
           &abserr, &neval, &ier, &limit, &lenw, &last, iwork.data(),
           work.data());
    
    // rethrow possible exceptions during function evaluation
    if (ex.second) std::rethrow_exception(ex.second);
    
    // throw `Rcpp::exception` if integration error occurred
    if (ier > 0) Rcpp::stop("Integration error");

    return result;
}
# R

## demo with sane input
gamma_inc_direct(-0.5, 1)
#> [1] 0.1781477

## demo with not so sane input
tryCatch(gamma_inc_direct(0, 0), error = function(cond) print(cond))
#> <Rcpp::exception in eval(expr, envir, enclos): Integration error>
tryCatch(gamma_inc_direct(0, -1), error = function(cond) print(cond))
#> <std::runtime_error in eval(expr, envir, enclos): non-finite function value>