jipolanco / BSplineKit.jl

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

Extrapolations? #19

Closed alecloudenback closed 1 year ago

alecloudenback commented 3 years ago

Is Extrapolations within the scope of this package's design goals? As an end user (for Yields.jl), I'm looking for a fast, pure-julia way to do 1-D interpolation and extrapolation (ideally with different conditions like Flat, Line, etc).

This seems almost as fast as Interpolations.jl but more flexible, as non-uniform cubic spline isn't supported AFAICT.

I'm thinking about simply wrapping a BSplineKit Interpolation and checking the bounds and performing my own extrapolation, including using the derivative for the line at the boundary condition.

I was wondering if there is a recommended way to: 1) recover the boundary of the spline from Interpolation? 2) recommended/prebuild way to find get the slope at the boundary?

jipolanco commented 3 years ago

Hi, thanks for your interest in BSplineKit!

Extrapolations were not one of the original goals of this package, but that doesn't mean they can't be supported. In fact I think it would be easy and useful to have them.

I will look into that when I can. As you suggest, a good solution would be to support the same boundary conditions defined in Interpolations.jl (Flat, Line, ...).

In the meantime, the following should allow you to get the value and slope at the (left) boundary. For now it's a bit verbose, but I could add some helper functions to make this easier (and also more efficient...):

using BsplineKit
using Random

xs = -1:0.1:1
ys = randn(length(xs))
itp = interpolate(xs, ys, BSplineOrder(4))

S = spline(itp)  # spline passing through data points
B = basis(S)     # B-spline basis

a, b = boundaries(B)  # left and right boundaries
@assert a == xs[begin] && b == xs[end]

# For now, we construct the full spline S′(x).
# There are faster ways of doing this that should be implemented...
S′ = diff(S, Derivative(1))

# Evaluate S and S′ at x = a
y_left = S(a)
yprime_left = S′(a)
alecloudenback commented 2 years ago

Thanks for the suggested workaround - I finally implemented this in https://github.com/JuliaActuary/Yields.jl/pull/85

The core code is very similar to what you proposed here: https://github.com/JuliaActuary/Yields.jl/blob/e36a755faf2eb1d725691d87fdcdd2802ded0820/src/utils.jl#L86-L133

dmbates commented 1 year ago

Let me second the request for evaluation of extrapolations, especially for cases like natural cubic spline interpolations. I occasionally have cases where a calculation produces a point just beyond the range of the interpolation x values and suddenly the spline value snaps to zero.

By the way, @jipolanco, thanks again for providing the periodic end conditions. If you are interested in why I wanted those, an example usage is at https://dmbates.quarto.pub/plotting-profiles

jipolanco commented 1 year ago

Interesting (and nice figures!), thanks for the link. Good to know periodic splines work well for you.

And thanks for the reminder, it should be easy to implement extrapolations. I'll try to add them within the next few days.

jipolanco commented 1 year ago

So I just implemented a basic extrapolation interface, available in v0.14. The usage is very close to that in Interpolations.jl:

itp = interpolate(-1:0.1:1, ...)
ext = extrapolate(itp, Flat())
ext(1.2)  # evaluate extrapolation

The documentation was also updated with some examples.

For now there are only two extrapolation strategies: Flat (same as in Interpolations.jl) and Smooth. Others can be easily added in the future if there's interest.