jipolanco / BSplineKit.jl

A collection of B-spline tools in Julia
https://jipolanco.github.io/BSplineKit.jl/dev/
MIT License
50 stars 9 forks source link

Added Linear() extrapolation method #93

Closed fintzij closed 3 months ago

fintzij commented 3 months ago

Added linear extrapolation method. The slope at the boundaries is calculated via finite differences.

(Hope the pull request is OK and this is helpful. Thanks for creating this package!)

codecov-commenter commented 3 months ago

Codecov Report

Attention: Patch coverage is 0% with 9 lines in your changes missing coverage. Please review.

Project coverage is 94.91%. Comparing base (a33f3c7) to head (ed6b81c). Report is 9 commits behind head on master.

Files Patch % Lines
src/SplineExtrapolations/SplineExtrapolations.jl 0.00% 9 Missing :warning:
Additional details and impacted files ```diff @@ Coverage Diff @@ ## master #93 +/- ## ========================================== - Coverage 95.61% 94.91% -0.71% ========================================== Files 35 35 Lines 2237 2224 -13 ========================================== - Hits 2139 2111 -28 - Misses 98 113 +15 ```

:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Have feedback on the report? Share it here.

jipolanco commented 3 months ago

Thanks, this can be useful indeed!

I'd just be a bit careful about using finite differences with eps() increments for at least two reasons.

First, because calling eps() without argument is the same as calling eps(1.0), which only makes sense when the boundary location a or b is approximately 1 (and is a Float64). For example, note that 10.0 + eps() == 10.0. The solution is to use eps(a) instead (or nextfloat/prevfloat).

The second problem is that, if one takes two "consecutive" numbers x and y = x + eps(x) and evaluate a function f(x) at these two points, it is possible for the two results to be the same. For example:

julia> x = 2.0
2.0

julia> y = x + eps(x)      # or y = nextfloat(x)
2.0000000000000004

julia> sqrt(x) == sqrt(y)  # true!
true

A solution would be to use the actual derivatives of the spline at the boundary locations to perform the extrapolations. Right now BSplineKit.jl doesn't allow to get a derivative at a single point x only, but only the full derivative of the spline which must then be evaluated at x.

This will be clearer with an example:

using BSplineKit
xs = 0.2:0.2:10.2
ys = sqrt.(xs)
S = interpolate(xs, ys, BSplineOrder(4))
a, b = boundaries(basis(S))

To obtain the derivative of S at a, we could do the following:

S′ = Derivative() * S
S′(a)  # 1.0738781882215265

But this is actually overkill as it requires first computing the spline derivative S′, which can be relatively expensive if the spline has many points. A more efficient solution would be to obtain S′(a) using automatic differentiation:

using ForwardDiff
ForwardDiff.derivative(S, a)  # 1.0738781882215263

In my opinion we can follow that last strategy. This will require adding ForwardDiff.jl as a dependency, but that's fine for me.

fintzij commented 3 months ago

Thanks, I agree that using ForwardDiff is a much better solution. I've pushed an update Linear() function but for some reason can't get the dependencies to work properly when I do a test run in the interpolations.jl file in the examples. Probably something to do with building the package locally. Does it work properly if you pull my changes in?

jipolanco commented 3 months ago

Not sure, but it may be related to the other changes made in Project.toml (besides adding ForwardDiff.jl).

fintzij commented 3 months ago

Yup, that'll do it. Everything is working on my end now, I think it's ready for review.

jipolanco commented 3 months ago

Looks good to me, thanks!

fintzij commented 3 months ago

Excellent, thank you too!