Fix the evaluation of bsplines infinitesimally close to the boundary

[EmilyBourne]

Fixes #599

[tpadioleau]

Waiting for the main branch to publish the docker images.

[tpadioleau]

@EmilyBourne We must use ncells instead of nbasis right ?

[EmilyBourne]

@EmilyBourne We must use ncells instead of nbasis right ?

No. nbasis is correct. The returned index should be the index of one of the basis elements. There are nbasis possibilities. ncells is much smaller in the non-periodic case. The problem is that rounding errors can go both ways. The second to last basis reaches 0 at the boundary so it is also affected by rounding errors. So the correct index is either the last or the second to last one. This is probably why I was just checking that we were returning a basis element (previously we returned Idx(nbasis) which is out of bounds).

[EmilyBourne]

@tpadioleau Do you prefer that I revert to checking that the index is in the acceptable domain or would you rather I check if it is one of the two plausible values?

[tpadioleau]

@EmilyBourne We must use ncells instead of nbasis right ?

No. nbasis is correct. The returned index should be the index of one of the basis elements. There are nbasis possibilities. ncells is much smaller in the non-periodic case. The problem is that rounding errors can go both ways. The second to last basis reaches 0 at the boundary so it is also affected by rounding errors. So the correct index is either the last or the second to last one. This is probably why I was just checking that we were returning a basis element (previously we returned Idx(nbasis) which is out of bounds).

I am surprised, shouldn't there be a difference of at most one between back_idx and bspl_basis_domain.back() ?

[tpadioleau]

@tpadioleau Do you prefer that I revert to checking that the index is in the acceptable domain or would you rather I check if it is one of the two plausible values?

Well in my opinion the more precise the test the better

[EmilyBourne]

@EmilyBourne We must use ncells instead of nbasis right ?

No. nbasis is correct. The returned index should be the index of one of the basis elements. There are nbasis possibilities. ncells is much smaller in the non-periodic case. The problem is that rounding errors can go both ways. The second to last basis reaches 0 at the boundary so it is also affected by rounding errors. So the correct index is either the last or the second to last one. This is probably why I was just checking that we were returning a basis element (previously we returned Idx(nbasis) which is out of bounds).

I am surprised, shouldn't there be a difference of at most one between back_idx and bspl_basis_domain.back() ?

Not sure what your question is. "second to last" means bspl_basis_domain.back()-1

[tpadioleau]

@EmilyBourne We must use ncells instead of nbasis right ?

No. nbasis is correct. The returned index should be the index of one of the basis elements. There are nbasis possibilities. ncells is much smaller in the non-periodic case. The problem is that rounding errors can go both ways. The second to last basis reaches 0 at the boundary so it is also affected by rounding errors. So the correct index is either the last or the second to last one. This is probably why I was just checking that we were returning a basis element (previously we returned Idx(nbasis) which is out of bounds).

I am surprised, shouldn't there be a difference of at most one between back_idx and bspl_basis_domain.back() ?

Not sure what your question is. "second to last" means bspl_basis_domain.back()-1

The CI shows there is a difference that corresponds to the degree of the splines.

337/762 Test #338: BSplinesFixture.Rounding_NonUniform<std::tuple<std::integral_constant<unsigned long, 5ul>, std::integral_constant<unsigned long, 100ul>, std::integral_constant<bool, false> >> ...........***Failed    0.02 sec
Note: Google Test filter = BSplinesFixture/34.Rounding_NonUniform
[==========] Running 1 test from 1 test suite.
[----------] Global test environment set-up.
[----------] 1 test from BSplinesFixture/34, where TypeParam = std::tuple<std::integral_constant<unsigned long, 5ul>, std::integral_constant<unsigned long, 100ul>, std::integral_constant<bool, false> >
[ RUN      ] BSplinesFixture/34.Rounding_NonUniform
/src/tests/splines/bsplines.cpp:161: Failure
Expected equality of these values:
  back_idx
    Which is: (99)
  bspl_basis_domain.back()
    Which is: (104)

[EmilyBourne]

@EmilyBourne We must use ncells instead of nbasis right ?

No. nbasis is correct. The returned index should be the index of one of the basis elements. There are nbasis possibilities. ncells is much smaller in the non-periodic case. The problem is that rounding errors can go both ways. The second to last basis reaches 0 at the boundary so it is also affected by rounding errors. So the correct index is either the last or the second to last one. This is probably why I was just checking that we were returning a basis element (previously we returned Idx(nbasis) which is out of bounds).

I am surprised, shouldn't there be a difference of at most one between back_idx and bspl_basis_domain.back() ?

Not sure what your question is. "second to last" means bspl_basis_domain.back()-1

The CI shows there is a difference that corresponds to the degree of the splines.

337/762 Test #338: BSplinesFixture.Rounding_NonUniform<std::tuple<std::integral_constant<unsigned long, 5ul>, std::integral_constant<unsigned long, 100ul>, std::integral_constant<bool, false> >> ...........***Failed    0.02 sec
Note: Google Test filter = BSplinesFixture/34.Rounding_NonUniform
[==========] Running 1 test from 1 test suite.
[----------] Global test environment set-up.
[----------] 1 test from BSplinesFixture/34, where TypeParam = std::tuple<std::integral_constant<unsigned long, 5ul>, std::integral_constant<unsigned long, 100ul>, std::integral_constant<bool, false> >
[ RUN      ] BSplinesFixture/34.Rounding_NonUniform
/src/tests/splines/bsplines.cpp:161: Failure
Expected equality of these values:
  back_idx
    Which is: (99)
  bspl_basis_domain.back()
    Which is: (104)

I have had a closer look. What I said earlier was incorrect. There is no rounding error because the point that I am testing here is out of bounds. The value that is returned should be a value idx such that the d+1 bspline elements [idx, idx+degree] can all be accessed. For a periodic domain this means that idx == bspl_basis_domain.back() as there are degree ghost elements. For a non-periodic domain idx == bspl_basis_domain.back() - degree. I have updated the test. It is now passing again

[tpadioleau]

I have had a closer look. What I said earlier was incorrect. There is no rounding error because the point that I am testing here is out of bounds. The value that is returned should be a value idx such that the d+1 bspline elements [idx, idx+degree] can all be accessed. For a periodic domain this means that idx == bspl_basis_domain.back() as there are degree ghost elements.

But we also do bspl_basis_domain.back() - BSplinesX::degree() even for the periodic case ?

[EmilyBourne]

I have had a closer look. What I said earlier was incorrect. There is no rounding error because the point that I am testing here is out of bounds. The value that is returned should be a value idx such that the d+1 bspline elements [idx, idx+degree] can all be accessed. For a periodic domain this means that idx == bspl_basis_domain.back() as there are degree ghost elements.

But we also do bspl_basis_domain.back() - BSplinesX::degree() even for the periodic case ?

Maybe bspl_basis_domain is badly named now. It is no longer restricted to the nbasis elements. It is the full domain

[tpadioleau]

For a periodic domain this means that idx == bspl_basis_domain.back() as there are degree ghost elements.

Did you mean idx == bspl_basis_domain.back() - BSplinesX::degree() ?