---
title: "Using splines2 with Rcpp"
author: Wenjie Wang
date: "`r Sys.Date()`"
bibliography:
  - ../inst/bib/splines2.bib
vignette: >
  %\VignetteIndexEntry{Using splines2 with Rcpp}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
output:
  rmarkdown::html_vignette:
    number_sections: yes
    toc: yes
---

# Introduction

In this vignette, we introduce how to use the C++ header-only library that
**splines2** contains with the **Rcpp** package [@eddelbuettel2013springer] 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** [@sanderson2016armadillo] library with
the help of **RcppArmadillo** [@eddelbuettel2014csda] 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.

```cpp
#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:

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

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

```cpp
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

```cpp
// 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.

```cpp
// 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.

```cpp
// 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).

```cpp
// 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:

```cpp
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.

```cpp
// 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`.

```cpp
// 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.

```cpp
// 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:

```cpp
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.

```cpp
// 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,

```cpp
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:

```cpp
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:

```cpp
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:

```cpp
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.

```cpp
// 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 {-}