Using splines2 with Rcpp
Wenjie Wang
20230821
Source:vignettes/splines2wircpp.Rmd
splines2wircpp.Rmd
Introduction
In this vignette, we introduce how to use the C++ headeronly 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 procedurebased functions in the R interface, the
C++ interface follows the commonlyused objectoriented 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 usingdirectives:
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 Bsplines; 
MSpline
for Msplines; 
ISpline
for Isplines; 
CSpline
for Csplines; 
NaturalSpline
andNaturalSplineK
for natural cubic splines; 
PeriodicMSpline
for periodic Msplines; 
PeriodicBSpline
for periodic Bsplines;
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 nondefault constructor is invoked when internal knots are
explicitly specified as the second argument. Taking Bsplines as an
example, the first nondefault constructor of a BSpline
object is
// 1. specify x, internal knots, degree, and boundary knots
(const vec& x,
BSplineconst vec& internal_knots,
const unsigned int degree = 3,
const vec& boundary_knots = vec());
The second nondefault 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
(const vec& x,
BSplineconst unsigned int spline_df,
const unsigned int degree = 3,
const vec& boundary_knots = vec());
The third nondefault 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
(const vec& x,
BSplineconst unsigned int degree,
const vec& knot_sequence);
The fourth nondefault 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
(const SplineBase* pSplineBase); BSpline
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:
{ &msp_obj }; BSpline bsp_obj
Constructors of PeriodicMSpline
and
PeriodicBSpline
The PeriodicMSpline
and PeriodicBSpline
classes are intended for constructing the periodic Msplines and
periodic Bsplines, respectively, which provide the same set of
nondefault constructors with BSpline
. The only difference
is that the knot sequence specified for the third nondefault
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 nondefault
constructor is called when internal knots are explicitly specified.
// 1. specify x, internal knots, and boundary knots
(const vec& x,
NaturalSplineconst vec& internal_knots,
const vec& boundary_knots = vec());
The second nondefault 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
(const vec& x,
NaturalSplineconst unsigned int spline_df,
const vec& boundary_knots = vec());
The third nondefault 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
(const SplineBase* pSplineBase); NaturalSpline
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 theCSpline
class)
The specific function signatures are as follows:
(const bool complete_basis = true);
mat basis(const unsigned int derivs = 1,
mat derivativeconst bool complete_basis = true);
(const bool complete_basis = true); mat integral
We can set and get the spline specifications through the following setter and getter functions, respectively.
// setter functions
* 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);
SplineBase
// getter functions
();
vec get_x();
vec get_internal_knots();
vec get_boundary_knots();
vec get_knot_sequenceunsigned 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,
{ arma::regspace(0, 0.1, 1) }; // 0, 0.1, ..., 1
vec x { x, 5 }; // df = 5 (and degree = 3, by default)
BSpline obj // change degree to 2 and get basis
{ obj.set_degree(2)>basis() }; mat basis_mat
The corresponding first derivatives and integrals of the basis functions can be obtained as follows:
{ bs.derivative() };
mat derivative_mat { bs.integral() }; mat integral_mat
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 nondefault constructor is as follows:
(const vec& x,
BernsteinPolyconst 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:
(const bool complete_basis = true);
mat basis(const unsigned int derivs = 1,
mat derivativeconst bool complete_basis = true);
(const bool complete_basis = true); mat integral
In addition, we may set and get the specifications through the following setter and getter functions, respectively.
// setter functions
* 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&);
BernsteinPoly
// getter functions
();
vec get_xunsigned int get_degree();
unsigned int get_order();
(); vec get_boundary_knots
The setter function returns a pointer to the current object.