# Using splines2 with Rcpp

#### Wenjie Wang

#### 2023-08-19

# 1 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.

# 2 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

# 3 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;

## 3.1 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 };

## 3.2 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.

## 3.3 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);

## 3.4 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`.

# 4 Generalized Bernstein Polynomials

The `BernsteinPoly` class is provided for the generalized Bernstein polynomials.

## 4.1 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.

## 4.2 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.