Skip to contents

Using splines2 with Rcpp

Wenjie Wang

2024-09-02

Source: vignettes/splines2-wi-rcpp.Rmd
splines2-wi-rcpp.Rmd

Introduction

In this vignette, we introduce how to use the C++ header-only library that splines2 contains with the Rcpp package (Eddelbuettel 2013) for constructing spline basis functions directly in C++. The introduction is intended for package developers who would like to use splines2 package in C++ by adding splines2 to the LinkingTo field of the package DESCRIPTION file.

Header Files and Namespace

Different from the procedure-based functions in the R interface, the C++ interface follows the commonly-used object-oriented design in C++ for ease of usage and maintenance. The implementations use the Armadillo (Sanderson 2016) library with the help of RcppArmadillo (Eddelbuettel and Sanderson 2014) and require C++11. We assume that C++11 is enabled and the header file named splines2Armadillo.h is included for access to all the classes and implementations in the namespace splines2 henceforth.

#include <RcppArmadillo.h>
#include <splines2Armadillo.h>  // include header files from splines2
// [[Rcpp::plugins(cpp11)]]

To use Rcpp::sourceCpp(), one may need to add [[Rcpp::depends()]] as follows:

// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::depends(splines2)]]

For ease of demonstration, we assume the following using-directives:

using namespace arma
using namespace splines2

Classes for Spline Basis Functions

A virtual base class named SplineBase is implemented to support a variety of classes for spline basis functions including

  • BSpline for B-splines;
  • MSpline for M-splines;
  • ISpline for I-splines;
  • CSpline for C-splines;
  • NaturalSpline and NaturalSplineK for natural cubic splines;
  • PeriodicMSpline for periodic M-splines;
  • PeriodicBSpline for periodic B-splines;

Constructors of BSpline, MSpline, ISpline, and CSpline

The BSpline, MSpline, ISpline, and CSpline classes share the same constructors inherited from the SplineBase class. There are four constructors in addition to the default constructor.

The first non-default constructor is invoked when internal knots are explicitly specified as the second argument. Taking B-splines as an example, the first non-default constructor of a BSpline object is

// 1. specify x, internal knots, degree, and boundary knots
BSpline(const vec& x,
        const vec& internal_knots,
        const unsigned int degree = 3,
        const vec& boundary_knots = vec());

The second non-default constructor is called when an unsigned integer is specified as the second argument, which represents the degree of freedom (DF) of the complete spline basis functions (different from the df argument in the R interface) is specified. Then the number of internal knots is computed as spline_df - degree - 1 and the placement of internal knots uses quantiles of specified x within the boundary.

// 2. specify x, spline DF, degree, and boundary knots
BSpline(const vec& x,
        const unsigned int spline_df,
        const unsigned int degree = 3,
        const vec& boundary_knots = vec());

The third non-default constructor is intended for the basis functions with an extended knot sequence, where the multiplicities of the knots can be more than one.

// 3. specify x, degree, and (extended) knot sequence
BSpline(const vec& x,
        const unsigned int degree,
        const vec& knot_sequence);

The fourth non-default constructor is explicit and takes a pointer to a base class object, which can be useful when we want to create a new object using the same specification (x, degree, internal_knots, and boundary_knots) of an existing object (not necessarily a BSpline object).

// 4. create a new object from a base class pointer
BSpline(const SplineBase* pSplineBase);

This constructor also allows us to easily switch between different types of splines. For example, we can create a BSpline object named bsp_obj from an existing MSpline object named msp_obj with the same specification as follows:

BSpline bsp_obj { &msp_obj };

Constructors of PeriodicMSpline and PeriodicBSpline

The PeriodicMSpline and PeriodicBSpline classes are intended for constructing the periodic M-splines and periodic B-splines, respectively, which provide the same set of non-default constructors with BSpline. The only difference is that the knot sequence specified for the third non-default constructor must be a simple knot sequence.

Constructors of NaturalSpline and NaturalSplineK

The classes NaturalSpline and NaturalSplineK are intended for natural cubic splines. The former corresponds to the function splines2::naturalSpline() (or splines2::nsp()) in R, while the latter is the engine of the function splines2::nsk(). They have the same constructors that do not allow the specification of the degree. Taking NaturalSpline as an example, the first non-default constructor is called when internal knots are explicitly specified.

// 1. specify x, internal knots, and boundary knots
NaturalSpline(const vec& x,
              const vec& internal_knots,
              const vec& boundary_knots = vec());

The second non-default constructor is called when an unsigned integer representing the degree of freedom of the complete spline basis functions (different from the df argument in the R interface) is specified. Then the number of internal knots is computed as spline_df - 2 and the placement of internal knots uses quantiles of specified x.

// 2. specify x, spline DF, and boundary knots
NaturalSpline(const vec& x,
              const unsigned int spline_df,
              const vec& boundary_knots = vec());

The third non-default constructor is explicit and takes a pointer to a base class object. It can be useful when we want to create a new object using the same specification (x, internal_knots, boundary_knots, etc.) of an existing object.

// 3. create a new object from a base class pointer
NaturalSpline(const SplineBase* pSplineBase);

Function Members

The main methods are

  • basis() for spline basis matrix
  • derivative() for derivatives of spline basis
  • integral() for integrals of spline basis (except for the CSpline class)

The specific function signatures are as follows:

mat basis(const bool complete_basis = true);
mat derivative(const unsigned int derivs = 1,
               const bool complete_basis = true);
mat integral(const bool complete_basis = true);

We can set and get the spline specifications through the following setter and getter functions, respectively.

// setter functions
SplineBase* set_x(const vec&);
SplineBase* set_x(const double);
SplineBase* set_internal_knots(const vec&);
SplineBase* set_boundary_knots(const vec&);
SplineBase* set_knot_sequence(const vec&);
SplineBase* set_degree(const unsigned int);
SplineBase* set_order(const unsigned int);

// getter functions
vec get_x();
vec get_internal_knots();
vec get_boundary_knots();
vec get_knot_sequence();
unsigned int get_degree();
unsigned int get_order();
unsigned int get_spline_df();

The setter function returns a pointer to the current object so that the specification can be chained for convenience. For example,

vec x { arma::regspace(0, 0.1, 1) }; // 0, 0.1, ..., 1
BSpline obj { x, 5 };                // df = 5 (and degree = 3, by default)
// change degree to 2 and get basis
mat basis_mat { obj.set_degree(2)->basis() };

The corresponding first derivatives and integrals of the basis functions can be obtained as follows:

mat derivative_mat { bs.derivative() };
mat integral_mat { bs.integral() };

Notice that there is no available integral() method for CSpline and no meaningful degree related methods for NaturalSpline.

Generalized Bernstein Polynomials

The BernsteinPoly class is provided for the generalized Bernstein polynomials.

Constructors

The main non-default constructor is as follows:

BernsteinPoly(const vec& x,
              const unsigned int degree,
              const vec& boundary_knots = vec());

In addition, two explicit constructors are provided for BernsteinPoly* and SplineBase*, which set x, degree, and boundary_knots from the objects that the pointers point to.

Function Members

The main methods are

  • basis() for the basis functions
  • derivative() for the derivatives of basis functions
  • integral() for the integrals of basis functions

The specific function signatures are as follows:

mat basis(const bool complete_basis = true);
mat derivative(const unsigned int derivs = 1,
               const bool complete_basis = true);
mat integral(const bool complete_basis = true);

In addition, we may set and get the specifications through the following setter and getter functions, respectively.

// setter functions
BernsteinPoly* set_x(const vec&);
BernsteinPoly* set_x(const double);
BernsteinPoly* set_degree(const unsigned int);
BernsteinPoly* set_order(const unsigned int);
BernsteinPoly* set_internal_knots(const vec&); // placeholder, does nothing
BernsteinPoly* set_boundary_knots(const vec&);

// getter functions
vec get_x();
unsigned int get_degree();
unsigned int get_order();
vec get_boundary_knots();

The setter function returns a pointer to the current object.

Reference

Eddelbuettel, Dirk. 2013. Seamless R and C++ Integration with Rcpp. Springer.
Eddelbuettel, Dirk, and Conrad Sanderson. 2014. “RcppArmadillo: Accelerating R with High-Performance C++ Linear Algebra.” Computational Statistics and Data Analysis 71: 1054–63.
Sanderson, Conrad. 2016. “Armadillo: An Open Source C++ Linear Algebra Library for Fast Prototyping and Computationally Intensive Experiments.” Journal of Open Source Software 1: 26.