Using cinterpolate

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 = "rich.fitzjohn@gmail.com")
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.