PCHIPInterpolation endpoints differences #415
Description
Describe the bug 🐞
Found new differences in results between the pchip implementation of DataInterpolations.jl and SciPy/Matlab. This is with DataInterpolations v8.0.0, after #400. This time the differences seem to occur only at the endpoints and when the "slopes" of the last two "intervals" are different (see MRE). This gets really bad when extrapolating.
Expected behavior
Results should be consistent between libraries.
Minimal Reproducible Example 👇
Case 1: Different slopes at the endpoints
Code for the figure
using DataInterpolations
using CondaPkg
using PythonCall
using MATLAB
## Interpolate
u = [0.9, 0.8, 0.86, 0.65, 0.44, 0.76, 0.73, 0.8]
t = [-2, -1, 0.0022, 0.68, 1.41, 2.22, 2.46, 2.76]
t_eval = t[1]:0.01:t[end]
# DataInterpolations.jl
A_datainterpolations = DataInterpolations.PCHIPInterpolation(u, t; extrapolation = ExtrapolationType.Extension)(t_eval)
# SciPy
# CondaPkg.add("scipy") # to install scipy in your local environment the first time you run the code
pyinterpolate = pyimport("scipy.interpolate")
A_python = pyconvert.(Array,pyinterpolate.PchipInterpolator(t, u)(t_eval))
# Matlab
using MATLAB
@mput t
@mput u
@mput t_eval
mat"A_matlab = interp1(t, u, t_eval,'pchip')"
@mget A_matlab
# Plot
using GLMakie
fig = Figure(size = (1600, 1000))
ax = Axis(fig[1:2, 1]; title = "PCHIP")
lines!(ax, t_eval, A_datainterpolations; linestyle = :dash, linewidth = 3, label = "DataInterpolations.jl")
lines!(ax, t_eval, A_python; linestyle = :solid, linewidth = 3, label = "SciPy")
lines!(ax, t_eval, A_matlab; linestyle = :dashdot, linewidth = 3, label = "MATLAB")
scatter!(ax, t, u, marker = :utriangle, markersize = 15, color = :black, label = "Data")
vlines!(ax, [t[1], t[end]], color = :black, linewidth = 2, linestyle = :dash)
axislegend()
ax = Axis(fig[1, 2]; title = "Difference")
scatterlines!(ax, t_eval, A_python .- A_datainterpolations; linewidth = 3, label = "Scipy - DataInterpolations.jl")
vlines!(ax, [t[1], t[end]], color = :black, linewidth = 2, linestyle = :dash)
axislegend()
ax = Axis(fig[2, 2]; title = "Difference", yticks = [-2e-16, 0, 2e-16])
scatterlines!(ax, t_eval, A_python .- A_matlab; linewidth = 3, label = "Scipy - MATLAB")
vlines!(ax, [t[1], t[end]], color = :black, linewidth = 2, linestyle = :dash)
axislegend()
display(fig)
Case 2: Different slopes at the endpoints + extrapolation
Code for the figure
using DataInterpolations
using CondaPkg
using PythonCall
using MATLAB
## Interpolate
u = [0.9, 0.8, 0.86, 0.65, 0.44, 0.76, 0.73, 0.8]
t = [-2, -1, 0.0022, 0.68, 1.41, 2.22, 2.46, 2.76]
t_eval = -4:0.01:3.5
# DataInterpolations.jl
A_datainterpolations = DataInterpolations.PCHIPInterpolation(u, t; extrapolation = ExtrapolationType.Extension)(t_eval)
# SciPy
# CondaPkg.add("scipy") # to install scipy in your local environment the first time you run the code
pyinterpolate = pyimport("scipy.interpolate")
A_python = pyconvert.(Array,pyinterpolate.PchipInterpolator(t, u)(t_eval))
# Matlab
using MATLAB
@mput t
@mput u
@mput t_eval
mat"A_matlab = interp1(t, u, t_eval,'pchip')"
@mget A_matlab
# Plot
using GLMakie
fig = Figure(size = (1600, 1000))
ax = Axis(fig[1:2, 1]; title = "PCHIP")
lines!(ax, t_eval, A_datainterpolations; linestyle = :dash, linewidth = 3, label = "DataInterpolations.jl")
lines!(ax, t_eval, A_python; linestyle = :solid, linewidth = 3, label = "SciPy")
lines!(ax, t_eval, A_matlab; linestyle = :dashdot, linewidth = 3, label = "MATLAB")
scatter!(ax, t, u, marker = :utriangle, markersize = 15, color = :black, label = "Data")
vlines!(ax, [t[1], t[end]], color = :black, linewidth = 2, linestyle = :dash)
axislegend()
ax = Axis(fig[1, 2]; title = "Difference")
scatterlines!(ax, t_eval, A_python .- A_datainterpolations; linewidth = 3, label = "Scipy - DataInterpolations.jl")
vlines!(ax, [t[1], t[end]], color = :black, linewidth = 2, linestyle = :dash)
axislegend()
ax = Axis(fig[2, 2]; title = "Difference", yticks = [-2e-16, 0, 2e-16])
scatterlines!(ax, t_eval, A_python .- A_matlab; linewidth = 3, label = "Scipy - MATLAB")
vlines!(ax, [t[1], t[end]], color = :black, linewidth = 2, linestyle = :dash)
axislegend()
display(fig)
Case 3: Same slopes at one end + extrapolation
Code for the figure
using DataInterpolations
using CondaPkg
using PythonCall
using MATLAB
## Interpolate
u = [0.9, 0.8, 0.76, 0.65, 0.44, 0.76, 0.73, 0.8]
t = [-2, -1, 0.0022, 0.68, 1.41, 2.22, 2.46, 2.76]
t_eval = -4:0.01:3.5
# DataInterpolations.jl
A_datainterpolations = DataInterpolations.PCHIPInterpolation(u, t; extrapolation = ExtrapolationType.Extension)(t_eval)
# SciPy
# CondaPkg.add("scipy") # to install scipy in your local environment the first time you run the code
pyinterpolate = pyimport("scipy.interpolate")
A_python = pyconvert.(Array,pyinterpolate.PchipInterpolator(t, u)(t_eval))
# Matlab
using MATLAB
@mput t
@mput u
@mput t_eval
mat"A_matlab = interp1(t, u, t_eval,'pchip')"
@mget A_matlab
# Plot
using GLMakie
fig = Figure(size = (1600, 1000))
ax = Axis(fig[1:2, 1]; title = "PCHIP")
lines!(ax, t_eval, A_datainterpolations; linestyle = :dash, linewidth = 3, label = "DataInterpolations.jl")
lines!(ax, t_eval, A_python; linestyle = :solid, linewidth = 3, label = "SciPy")
lines!(ax, t_eval, A_matlab; linestyle = :dashdot, linewidth = 3, label = "MATLAB")
scatter!(ax, t, u, marker = :utriangle, markersize = 15, color = :black, label = "Data")
vlines!(ax, [t[1], t[end]], color = :black, linewidth = 2, linestyle = :dash)
axislegend()
ax = Axis(fig[1, 2]; title = "Difference")
scatterlines!(ax, t_eval, A_python .- A_datainterpolations; linewidth = 3, label = "Scipy - DataInterpolations.jl")
vlines!(ax, [t[1], t[end]], color = :black, linewidth = 2, linestyle = :dash)
axislegend()
ax = Axis(fig[2, 2]; title = "Difference", yticks = [-2e-16, 0, 2e-16])
scatterlines!(ax, t_eval, A_python .- A_matlab; linewidth = 3, label = "Scipy - MATLAB")
vlines!(ax, [t[1], t[end]], color = :black, linewidth = 2, linestyle = :dash)
axislegend()
display(fig)
