Using cinterpolate

This package provides a minimal set of interpolation methods (piecewise constant, linear and spline) designed to be compatible with R’s approx and spline functions but callable from C. It will only be of interest to people writing C or C++ packages.

The package is designed to be used with R’s LinkingTo: support and is header only. This is a somewhat awkward situation for C (rather than C++). The approach is the same as taken by ring.

Package preparation

  • In your DESCRIPTION, add a line LinkingTo: cinterpolate. You also need to add cinterpolate to Imports and ensure that the package is loaded before use, as it uses R’s R_GetCCallable interface to call the actual interpolation functions.

  • In your src/ directory, add a file cinterpolate.c containing just the line #include <cinterpolate/cinterpolate.c>.

  • Anywhere in your C (or C++) code you want to use the interpolation code buffer, include the line #include <cinterpolate/cinterpolate.h> to include the prototypes and use the interface as described below.

(I am not sure what the best practice way of doing this with a standalone shared library compiled with R CMD SHLIB is though; probably best to make a package.)

The API

There are only three functions in the cinterpolate API; one to build the object (cinterpolate_alloc), one for carrying out interpolation (cinterpolate_eval) and one for freeing the object after calculations have been run (cinterpolate_free). If you allocate a cinterpolate object then you are responsible for freeing it (even on error elsewhere in code). Not doing this will cause leaks.

#ifndef CINTERPOLTE_CINTERPOLATE_H_
#define CINTERPOLTE_CINTERPOLATE_H_
#include <stddef.h> // size_t
#include <stdbool.h> // bool

// Allow use from C++
#ifdef __cplusplus
extern "C" {
#endif

// There are only three functions in the interface; allocation,
// evaluation and freeing.

// Allocate an interpolation object.
//
//   type: The mode of interpolation. Must be one of "constant",
//       "linear" or "spline" (an R error is thrown if a different
//       value is given).
//
//   n: The number of `x` points to interpolate over
//
//   ny: the number of `y` points per `x` point.  This is 1 in the
//       case of zimple interpolation as used by Rs `interpolate()`
//
//   x: an array of `x` values of length `n`
//
//   y: an array of `ny` sets of `y` values.  This is in R's matrix
//       order (i.e., the first `n` values are the first series to
//       interpolate over).
//
//   fail_on_extrapolate: if true, when an extrapolation occurs throw
//       an error; if false return NA_REAL
//
//   auto_free: automatically clean up the interpolation object on
//       return to R. This uses `R_alloc` for allocations rather than
//       `Calloc` so freeing will always happen (even on error
//       elsewhere in the code). However, this prevents returning back
//       a pointer to R that will last longer than the call into C
//       code.
//
// The return value is an opaque pointer that can be passed through to
// `cinterpolate_eval` and `cinterpolate_free`
void *cinterpolate_alloc(const char *type, size_t n, size_t ny,
                         double *x, double *y, bool fail_on_extrapolate,
                         bool auto_free);

// Evaluate the interpolated function at a new `x` point.
//
//   x: A new, single, `x` point to interpolate `y` values to
//
//   obj: The interpolation object, as returned by `cinterpolate_alloc`
//
//   y: An array of length `ny` to store the interpolated values
//
// The return value is 0 if the interpolation is successful (with x
// lying within the range of values that the interpolation function
// supports), -1 otherwise
int cinterpolate_eval(double x, void *obj, double *y);

// Clean up all allocated memory
//
//   obj: The interpolation object, as returned by `cinterpolate_alloc`
void cinterpolate_free(void *obj);

#ifdef __cplusplus
}
#endif

#endif

A complete example of use is included in the package as system.file("example", package = "cinterpolate").

The DESCRIPTION looks like

Package: example
Title: Testing 'cinterpolate' Package Use
Version: 0.0.1
Description: Testing that using 'cinterpolate' from another package works.
License: CC0
Authors@R: person("Rich", "FitzJohn", role = c("aut", "cre"),
    email = "[email protected]")
LinkingTo: cinterpolate
Imports: cinterpolate
Suggests: testthat

Note the use of LinkingTo: and Imports: here.

The NAMESPACE file ensures that the package’s shared library is loaded (useDynLib(example)) and that cinterpolate’s functions will be available by importing the package import(cinterpolate) (importFrom(cinterpolate, interpolate_function) would also be fine).

useDynLib(example)
import(cinterpolate)

The actual usage from C looks like:

#include <cinterpolate/cinterpolate.c>
#include <stdbool.h>
#include <R.h>
#include <Rinternals.h>
#include <R_ext/Utils.h>

SEXP test(SEXP r_x, SEXP r_y, SEXP r_xout, SEXP r_type) {
  const char * type = CHAR(STRING_ELT(r_type, 0));

  bool is_matrix = isMatrix(r_y);

  size_t n = length(r_x);
  size_t ny = is_matrix ? ncols(r_y) : 1;
  void *obj = cinterpolate_alloc(type, n, ny, REAL(r_x), REAL(r_y),
                                 false, true);

  size_t nxout = length(r_xout);
  SEXP r_yout = PROTECT(is_matrix ?
                        allocMatrix(REALSXP, ny, nxout) :
                        allocVector(REALSXP, nxout));
  double *xout = REAL(r_xout), *yout = REAL(r_yout);
  for (size_t i = 0; i < nxout; ++i) {
    cinterpolate_eval(xout[i], obj, yout);
    yout += ny;
  }

  UNPROTECT(1);
  return r_yout;
}

  • Because this is the only use of cinterpolate in the package, we can directly include cinterpolate/cinterpolate.c

  • The interpolate_alloc function is used with auto_clean = true so there is no use of interpolate_free - because this interpolator only needs to survive for this single C function this method of cleanup is probably better

  • There is a single allocation of an interpolation object but several calls to interpolate new values.