Using cinterpolate

Rich FitzJohn

2023-11-29

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

(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;
}