1-point ConstantInterpolation gets BoundsError with cache_parameters

[visr]

Describe the bug 🐞

A ConstantInterpolation works fine for single data points as long as extrapolation is enabled. However it cannot be constructed when using cache_parameters.

Expected behavior

I expect cache_parameters to work just like with the setting turned off.

Minimal Reproducible Example 👇

using DataInterpolations
ConstantInterpolation([2], [0.0]; extrapolate = true, cache_parameters = true)

Error & Stacktrace ⚠️

ERROR: BoundsError: attempt to access 0-element Vector{Float64} at index [1]
Stacktrace:
 [1] getindex
   @ .\essentials.jl:13 [inlined]
 [2] first
   @ .\abstractarray.jl:452 [inlined]
 [3] cumulative_integral
   @ a:\.julia\packages\DataInterpolations\C8cNn\src\interpolation_utils.jl:196 [inlined]
 [4] ConstantInterpolation(u::Vector{…}, t::Vector{…}; dir::Symbol, extrapolate::Bool, cache_parameters::Bool, assume_linear_t::Float64)
   @ DataInterpolations a:\.julia\packages\DataInterpolations\C8cNn\src\interpolation_caches.jl:284

Environment (please complete the following information):

• Output of using Pkg; Pkg.status()

DataInterpolations v6.6.0

• Output of versioninfo()

Julia Version 1.10.7
Commit 4976d05258 (2024-11-26 15:57 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: 24 × 12th Gen Intel(R) Core(TM) i7-12800HX
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-15.0.7 (ORCJIT, alderlake)
Threads: 12 default, 0 interactive, 6 GC (on 24 virtual cores)
Environment:
  JULIA_DEPOT_PATH = a:\.julia;
  JULIA_EDITOR = code
  JULIA_NUM_THREADS = 12
  JULIA_PKG_PRESERVE_TIERED_INSTALLED = true
  JULIA_PKG_USE_CLI_GIT = true
  JULIA_PROJECT = @.

[sathvikbhagavan]

I think it is fixed on master.

julia> ConstantInterpolation([2], [0.0]; extrapolation=DataInterpolations.ExtrapolationType.Constant, cache_parameters = true)
ConstantInterpolation with 1 points, in left direction
┌──────┬─────┐
│ time │   u │
├──────┼─────┤
│  0.0 │ 2.0 │
└──────┴─────┘

We just need to do a release. @ChrisRackauckas, can you release it? This would be a breaking release as extrapolate keyword was removed and extrapolation was added for various extrapolation options in #356

[visr]

Thanks, I should've checked master, just released as v7. Indeed the construction now works regardless of cache_parameters. However, it seems now extrapolation on either side is broken:

itp = ConstantInterpolation([2], [0.0]; extrapolation=ExtrapolationType.Constant)
itp(1.0)

ERROR: InexactError: Int64(NaN)
Stacktrace:
 [1] Int64(x::Float64)
   @ Base .\float.jl:912
 [2] _derivative
   @ a:\.julia\dev\DataInterpolations\src\derivatives.jl:162 [inlined]
 [3] derivative(A::ConstantInterpolation{Vector{…}, Vector{…}, Vector{…}, Int64, (1,)}, t::Float64, order::Int64)
   @ DataInterpolations a:\.julia\dev\DataInterpolations\src\derivatives.jl:9
 [4] derivative
   @ a:\.julia\dev\DataInterpolations\src\derivatives.jl:2 [inlined]
 [5] _extrapolate_right(A::ConstantInterpolation{Vector{Int64}, Vector{Float64}, Vector{Float64}, Int64, (1,)}, t::Float64)
   @ DataInterpolations a:\.julia\dev\DataInterpolations\src\interpolation_methods.jl:38
 [6] _interpolate(A::ConstantInterpolation{Vector{Int64}, Vector{Float64}, Vector{Float64}, Int64, (1,)}, t::Float64)
   @ DataInterpolations a:\.julia\dev\DataInterpolations\src\interpolation_methods.jl:5
 [7] (::ConstantInterpolation{Vector{Int64}, Vector{Float64}, Vector{Float64}, Int64, (1,)})(t::Float64)
   @ DataInterpolations a:\.julia\dev\DataInterpolations\src\DataInterpolations.jl:27

This worked on v6.6. I get the same error for two datapoints with constant u, itp = ConstantInterpolation([2,2], [0.0,0.1]; extrapolation=ExtrapolationType.Constant). And if u is not constant it extrapolates to NaN: itp = ConstantInterpolation([2,2.1], [0.0,0.1]; extrapolation=ExtrapolationType.Constant).

[ChrisRackauckas]

@SouthEndMusic can you look into this case?

[SouthEndMusic]

I looked into it, and the problem comes from derivative calculation here:

DataInterpolations.jl/src/derivatives.jl

Lines 161 to 163 in f37ff2a

	function _derivative(A::ConstantInterpolation{<:AbstractVector}, t::Number, iguess)
	return isempty(searchsorted(A.t, t)) ? zero(A.u[1]) : eltype(A.u)(NaN)
	end

This code makes sure that evaluating the derivative of a constant interpolation in a datapoint yields a NaN of the appropriate type. This is triggered because extrapolation specifically evaluates the derivative at the boundary data points. However, your eltype of u is Int64, which cannot be converted to NaN.

I'm unfamiliar with the history how this came to be, but I suggest just returning a derivative of 0 in the boundary points when eltype(A.u) <: Integer.

The sad thing is that we don't even need derivative calculations for constant extrapolation numerically, it's only added in the extrapolation code for type stability reasons.

P.s.: @visr I was going to say that I was happy that caching was picked up by more users before I saw your name, you don't count 🙃