Create periodic B-Splines using Spline2 R-Package

periodic_bsplines.cpp

#include <RcppArmadillo.h>
// include header file from splines2 package
#include <splines2Armadillo.h>

// [[Rcpp::export]]
arma::mat splines2_periodic(const arma::vec &x,
                            const arma::vec &knots,
                            const unsigned int deg,
                            const bool &intercept = true)
{

    unsigned int order = deg + 1;
    arma::uvec inner_idx = arma::regspace<arma::uvec>(order,
                                                      knots.n_elem - order - 1);
    arma::uvec bound_idx = {deg, knots.n_elem - order};

    // Create periodic mspline object
    splines2::PeriodicMSpline ps(x, knots(inner_idx), deg, knots(bound_idx));
    arma::mat ps_mat = ps.basis(true);

    if (!intercept)
        ps_mat.shed_col(0);

    // We will use https://doi.org/10.6339/21-JDS1020 eq. (1) to convert
    // the mspline basis to a bspline basis

    // We need this sequence to calculate the weights
    arma::vec knots_ext = knots.subvec(bound_idx(0), bound_idx(1));
    knots_ext = arma::join_cols(knots_ext,
                                knots(inner_idx.head(deg)) + knots(bound_idx(1)) - knots(bound_idx(0)));

    for (unsigned int i = 0; i < ps_mat.n_cols; i++)
    {
        double w = knots_ext(1 - intercept + i + order) -
                   knots_ext(1 - intercept + i);
        ps_mat.col(i) *= w / order;
    }

    return ps_mat;
}

/*** R

x <- 0:100 / 100
deg <- 2

knots <- c(-0.2, -0.1, 0.0, 0.1, 0.2, 0.5, 1.0, 1.5, 2.1)

res <- splines2_periodic(x, knots, deg, TRUE)
ts.plot(res)
rowSums(res)
*/

script.R

Rcpp::sourceCpp("periodic_bsplines.cpp")