Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Error in set_knot_sequence for PeriodicMSpline #18

Closed
BerriJ opened this issue Mar 22, 2023 · 5 comments
Closed

Error in set_knot_sequence for PeriodicMSpline #18

BerriJ opened this issue Mar 22, 2023 · 5 comments

Comments

@BerriJ
Copy link

BerriJ commented Mar 22, 2023

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
Copy link
Owner

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
Copy link
Author

BerriJ commented Mar 22, 2023

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
Copy link
Owner

wenjie2wang commented Mar 23, 2023

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
Copy link
Owner

wenjie2wang commented Mar 25, 2023

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
Copy link
Owner

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants