Safer knot deduplication

[mateuszbaran]

This small PR ensures that more pathological cases of duplicated knots (for example [1.0, 1.0, 1.0, nextfloat(1.0), nextfloat(1.0)]) are correctly handled.

[codecov]

Codecov Report

Merging #468 (45ebf78) into master (eefe448) will decrease coverage by 8.69%.
The diff coverage is 100.00%.

❗ Current head 45ebf78 differs from pull request most recent head 6349308. Consider uploading reports for the commit 6349308 to get more accurate results

@@            Coverage Diff             @@
##           master     #468      +/-   ##
==========================================
- Coverage   85.11%   76.42%   -8.70%     
==========================================
  Files          26       26              
  Lines        1747     1421     -326     
==========================================
- Hits         1487     1086     -401     
- Misses        260      335      +75     

Impacted Files	Coverage Δ
src/gridded/gridded.jl	86.00% <100.00%> (-5.53%)	⬇️
src/requires/unitful.jl	0.00% <0.00%> (-100.00%)	⬇️
src/nointerp/nointerp.jl	44.44% <0.00%> (-35.56%)	⬇️
src/utils.jl	41.66% <0.00%> (-33.02%)	⬇️
src/b-splines/linear.jl	73.33% <0.00%> (-26.67%)	⬇️
src/b-splines/constant.jl	76.00% <0.00%> (-24.00%)	⬇️
src/b-splines/b-splines.jl	68.13% <0.00%> (-21.87%)	⬇️
src/b-splines/indexing.jl	50.96% <0.00%> (-15.21%)	⬇️
src/scaling/scaling.jl	78.37% <0.00%> (-12.21%)	⬇️
src/extrapolation/extrapolation.jl	61.76% <0.00%> (-9.03%)	⬇️
... and 16 more

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update eefe448...6349308. Read the comment docs.

[mateuszbaran]

Hm... Julia 1.7 has a regression in type inference related to ndims. I think I can only fix it by forcing more method specialization. What do you think?

[mkitti]

I am still not convinced that we are honoring the user's intentions here. In your example we have two unique knots: 1.0 and nextfloat(1.0). The interpolated values at each of the unique knots should be one of the values given for those knots. Let's say the values are [1,2,3,4,5].

julia> A = [1.0, 1.0, 1.0, nextfloat(1.0), nextfloat(1.0)]
5-element Vector{Float64}:
 1.0
 1.0
 1.0
 1.0000000000000002
 1.0000000000000002

julia> B = deduplicate_knots!(copy(A))
5-element Vector{Float64}:
 1.0
 1.0000000000000002
 1.0000000000000004
 1.0000000000000007
 1.0000000000000009

julia> V = [1,2,3,4,5]
5-element Vector{Int64}:
 1
 2
 3
 4
 5

itp = interpolate((B,), V, Gridded(Linear()))

julia> collect(zip(A, itp.(A)))
5-element Vector{Tuple{Float64, Float64}}:
 (1.0, 1.0)
 (1.0, 1.0)
 (1.0, 1.0)
 (1.0000000000000002, 2.0)
 (1.0000000000000002, 2.0)

Above 1.0 is mapped to 1 which is fine. That's the first value I gave to the knot at 1.0.
nextfloat(1.0) however, is mapped to 2.0. However, the values I assigned to nextfloat(1.0) are either 3, 4, or 5.

My expectation is that itp(nextfloat(1.0)) in V[3:5] be true, but it is not.

We may thus need a two argument version of deduplicate_knots! where the values are also given and some of them are ommitted.

The behavior of itp below is closer to what I might expect from deduplicating this.

julia> V = [1,4,5]
3-element Vector{Int64}:
 1
 4
 5

julia> B = [1.0, nextfloat(1.0), nextfloat(nextfloat(1.0))]
3-element Vector{Float64}:
 1.0
 1.0000000000000002
 1.0000000000000004

julia> itp = interpolate((B,), V, Gridded(Linear()))
3-element interpolate((::Vector{Float64},), ::Vector{Int64}, Gridded(Linear())) with element type Float64:
 1.0
 4.0
 5.0

julia> itp.(A)
3-element Vector{Float64}:
 1.0
 4.0
 5.0

julia> collect(zip(A, itp.(A)))
3-element Vector{Tuple{Float64, Float64}}:
 (1.0, 1.0)
 (1.0000000000000002, 4.0)
 (1.0000000000000004, 5.0)

nextfloat(1.0) now maps to one of the values originally associated with that knot, 4.0. If there were further values, then values greater than nextfloat(1.0) might map to 5.0. Here we had to not only change the knots, but we also had to change the values and drop some of the knot-value pairs.

[mateuszbaran]

Right, my solution is not perfect (though it's fine for my use case). Actually I'm not sure what the best solution would be. Consider knots [prevfloat(1.0), 1.0, 1.0, nextfloat(1.0)] and values [0.0, -1.0, 1.0, 0.0]. It's impossible to correct it because there is no "right" value to pick for knot 1.0. The sloppy solution is doing what I've done in this PR. If we can remove knots such that the range of values is retained (for example we could drop the second knot if values in my example were [0.0, 0.5, 1.0, 0.0], but not the third one) I think we can reasonably do so. Otherwise the user has to fix the data themselves.

[mkitti]

In general, I think the user should probably fix the data themselves. I consider deduplicate_knots! to be merely a convenience function. In the latest example, the value at the knot 1.0 should be either -1.0 or 1.0. The circumstances force us to drop a knot to be the most consistent.

[N5N3]

As for ndims, eltype and so on, I think we'd better replace the following:

Interpolations.jl/src/Interpolations.jl

Lines 150 to 161 in eefe448

	itptype(::Type{AbstractInterpolation{T,N,IT}}) where {T,N,IT<:DimSpec{InterpolationType}} = IT
	itptype(::Type{ITP}) where {ITP<:AbstractInterpolation} = itptype(supertype(ITP))
	itptype(itp::AbstractInterpolation) = itptype(typeof(itp))
	gridtype(::Type{AbstractInterpolation{T,N,IT}}) where {T,N,IT<:DimSpec{InterpolationType}} = GT
	gridtype(::Type{ITP}) where {ITP<:AbstractInterpolation} = gridtype(supertype(ITP))
	gridtype(itp::AbstractInterpolation) = gridtype(typeof(itp))
	ndims(::Type{AbstractInterpolation{T,N,IT}}) where {T,N,IT<:DimSpec{InterpolationType}} = N
	ndims(::Type{ITP}) where {ITP<:AbstractInterpolation} = ndims(supertype(ITP))
	ndims(itp::AbstractInterpolation) = ndims(typeof(itp))
	eltype(::Type{AbstractInterpolation{T,N,IT}}) where {T,N,IT<:DimSpec{InterpolationType}} = T
	eltype(::Type{ITP}) where {ITP<:AbstractInterpolation} = eltype(supertype(ITP))
	eltype(itp::AbstractInterpolation) = eltype(typeof(itp))

with:

itptype(::Type{<:AbstractInterpolation{T,N,IT}}) where {T,N,IT<:DimSpec{InterpolationType}} = IT
itptype(itp::AbstractInterpolation) = itptype(typeof(itp))
ndims(::Type{<:AbstractInterpolation{T,N,<:DimSpec{InterpolationType}}}) where {T,N} = N
ndims(itp::AbstractInterpolation) = ndims(typeof(itp))
eltype(::Type{<:AbstractInterpolation{T,N,<:DimSpec{InterpolationType}}}) where {T,N} = T
eltype(itp::AbstractInterpolation) = eltype(typeof(itp))

I delete gridtype as it's broken (GT is undefined)

[mkitti]

@N5N3, I'm unclear how this is related. Could you open a new issue or PR?

[mkitti]

To move this PR forward, could we implement this behavior as an keyword option? I think this may be useful to others, but I do not think it should be the default behavior.

[mateuszbaran]

This PR is a bit stuck on unrelated Julia 1.7 failures so I was planning to go back to it once that is solved in some way. I'm also not quite sure at the moment how exactly this should work.

[mkitti]

Let me worry about 1.7.

[mkitti]

I ended up going with a warning:

julia> Interpolations.deduplicate_knots!([1.0, 1.0, 1.0, nextfloat(1.0), nextfloat(1.0)])
┌ Warning: Successive repeated knots detected. Consider using `move_knots` keyword to Interpolations.deduplicate_knots!
│   last_knot = 1.0
│   knots[i - 1] = 1.0000000000000004
│   knots[i] = 1.0000000000000002
└ @ Interpolations C:\Users\kittisopikulm\.julia\dev\Interpolations\src\gridded\gridded.jl:120
5-element Vector{Float64}:
 1.0
 1.0000000000000002
 1.0000000000000004
 1.0000000000000007
 1.0000000000000009

[mkitti]

Fix #467

[mateuszbaran]

Thanks for finishing it! I was still considering adding the node removal feature here but I was busy with other work.

[jaemolihm]

@mkitti Does this PR fix #107 ?