-
-
Notifications
You must be signed in to change notification settings - Fork 46
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
Unstable quadratic spline #353
Comments
@SouthEndMusic you were the last one to look at the indexing of quadratic, seems like it could be off by one or something in one of the dispatches? |
The quadratic spline interpolation does what it is supposed to: it's You could represent this data better with a |
Hello, thank you for your answer!
|
Hi Marco, When it comes to detecting an instability: it depends on how you do define that. In general the best strategy is to choose the type of interpolation that has the properties that you want. As to why |
One of the weird things this interpolation type does is assume that the derivative in the first point is 0, throwing away a perfectly good degree of freedom. |
@SouthEndMusic , what do you mean with "this interpolation type"? The quadratic spline in general, the algorithm implemented here...? |
@marcobonici the algorithm for I hardcoded a degree 2 B-spline fit to your data, and it works well: using DataInterpolations
using GLMakie
# Don't use x[end-1] in knot vector
# This is a choice, but something has to be done to make
# the number of degrees of freedom match the amount of data
k = zeros(length(x) + 3)
k[1:3] .= x[1]
k[(end-2):end] .= x[end]
k[4:(end-3)] .= x[2:(end-2)]
function B(t, i)
return if t < k[i]
zero(t)
elseif t <= k[i+1]
enum = (t - k[i])^2
denom = (k[i+2] - k[i]) * (k[i+1] - k[i])
if (enum == 0 && denom == 0)
i == 1 ? one(t) : zero(t)
else
enum/denom
end
elseif t <= k[i+2]
(t - k[i]) * (k[i+2] - t) / ((k[i+2] - k[i]) * (k[i+2] - k[i+1])) +
(k[i+3] - t) * (t - k[i+1]) / ((k[i+3] - k[i+1]) * (k[i+2] - k[i+1]))
elseif t <= k[i+3]
((k[i+3] - t)^2) / ((k[i+3] - k[i+1]) * (k[i+3] - k[i+2]))
else
zero(t)
end
end
# A * c = y
n = length(x)
A = zeros(n, n)
for j in 1:n
for i in 1:n
A[j,i] = B(x[j], i)
end
end
c = A \ y
function itp(t)
out = zero(t)
for i in 1:n
out += c[i] * B(t, i)
end
out
end
new_x = LinRange(x[1], x[end], 10000)
fig = Figure()
ax_bf = Axis(fig[1,1], title = "Basis functions")
scatter!(ax_bf, x, zero(x))
for i in 1:n
lines!(ax_bf, new_x, B.(new_x, i))
end
ax_fit = Axis(fig[2,1], title = "Interpolation")
lines!(ax_fit, new_x, itp.(new_x))
scatter!(ax_fit, x,y, label = "Data")
fig We should probably incorporate something like this in I don't know exactly what the current |
I have a POC of an improved lightweight implementation, I will open a PR when I find the time. |
Thank you very much! :) |
Describe the bug 🐞
When using the quadratic spline on some smooth data, I get some oddly unstable behaviour which I do not get when using Linear, Cubic, or Akima interpolation.
scipy
quadratic spline doesn't give the problems I am getting here.Expected behavior
I would expect the interpolated function to be smooth.
Minimal Reproducible Example 👇
Here is the MWE, with the input data and the few lines of code used to interpolate.
Error & Stacktrace⚠️
There is no stack trace, the code is just giving a weird output.
Environment (please complete the following information):
using Pkg; Pkg.status()
using Pkg; Pkg.status(; mode = PKGMODE_MANIFEST)
versioninfo()
The text was updated successfully, but these errors were encountered: