wenjie2wang / splines2

Regression Spline Functions and Classes
https://wwenjie.org/splines2
GNU General Public License v3.0
42 stars 3 forks source link

Error in set_knot_sequence for PeriodicMSpline #18

Closed BerriJ closed 1 year ago

BerriJ commented 1 year ago

Hi,

First of all, thanks for this fantastic package.

At the moment, I'm working with splines2::PeriodicMSpline. Following the docs there is a method set_knot_sequence and I assumed that this can be used to set the entire knot sequence (so inner knots + boundary knots). However, this fails with an error. Consider the following small example:

#include <RcppArmadillo.h>
#include <splines2Armadillo.h>

// [[Rcpp::export]]
arma::mat splines_periodic(const bool &knot_sequence)
{

    int deg = 3;
    arma::vec x = arma::linspace(0, 1, 1000);

    // Inner + boundary knots
    arma::vec knots = {0, 0.2, 0.4, 0.6, 0.8, 1};

    // Create periodic spline object
    splines2::PeriodicMSpline ps_obj;
    ps_obj.set_x(x);
    ps_obj.set_degree(deg);

    if (knot_sequence)
    {
        // This fails with: "The length of specified knot sequence is too small."
        ps_obj.set_knot_sequence(knots);
    }
    else
    {
        // This works:
        ps_obj.set_internal_knots(knots.subvec(1, 4));
        ps_obj.set_boundary_knots({knots(0), knots(5)});
    }
    return ps_obj.basis(true);
}

Am I interpreting the functionality of set_knot_sequence wrong?

Thanks so much.

BerriJ

wenjie2wang commented 1 year ago

Thanks for the question!

I did not try the example code. But could you try the following for using the set_knot_sequence?

arma::vec knots = {0, 0, 0, 0, 0.2, 0.4, 0.6, 0.8, 1, 1, 1, 1};

I believe what you want to specify is a simple knot sequence (also called clamped knot vector), where each boundary knot is repeated degree + 1 times.

BerriJ commented 1 year ago

Thank you for your quick answer. However, I can not get this to work.

Using your suggestion leads to: Error in splines_periodic(TRUE) : Mat::tail_cols(): size out of bounds. I removed the first and last element in your suggested knot sequence (so that the boundary knots are repeated only degree times, but that leads to Expected a simple knot sequence..

wenjie2wang commented 1 year ago

Thanks for the quick feedback! I had a closer look and I think it is indeed unexpected. I will fix it once I get chance.

wenjie2wang commented 1 year ago

This issue should be resolved by the referenced commit. Consider the following example:

// issue-18.cpp

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

#include <RcppArmadillo.h>
#include <splines2Armadillo.h>

// [[Rcpp::export]]
arma::mat splines_periodic()
{
    unsigned int deg = 3;
    arma::vec x = arma::linspace(0, 1, 1000);

    // a simple knot sequence = inner + repeated boundary knots
    arma::vec knots = {0, 0, 0, 0, 0.2, 0.4, 0.6, 0.8, 1, 1, 1, 1};

    // Create periodic spline object
    splines2::PeriodicMSpline ps_obj;
    ps_obj.set_x(x);
    ps_obj.set_degree(deg);

    // it must be a simple knot sequence
    ps_obj.set_knot_sequence(knots);

    return ps_obj.basis(true);
}

// [[Rcpp::export]]
arma::mat splines_periodic2()
{
    unsigned int deg = 3;
    arma::vec x = arma::linspace(0, 1, 1000);

    // a simple knot sequence = inner + repeated boundary knots
    arma::vec knots = {0, 0, 0, 0, 0.2, 0.4, 0.6, 0.8, 1, 1, 1, 1};

    // Create periodic spline object
    splines2::PeriodicMSpline ps_obj { x, deg, knots };

    return ps_obj.basis(true);
}
## issue-18.R

Rcpp::sourceCpp("issue-18.cpp")

res <- splines_periodic()
matplot(res, type = "l")

res2 <- splines_periodic2()
all.equal(res, res2)

For PeriodicMSpline, the specified knot sequence must be a simple knot sequence. For consistency, I added a constructor to create a PeriodicMSpline object with a simple knot sequence directly.

wenjie2wang commented 1 year ago

The latest version of {splines2} (v0.4.8) is available on CRAN.