StackOverflowError in gradient computation

[gianmariomanca]

using Interpolations
xs = 0.5:0.1:5
A = [log(x) for x in xs]
interp_cubic = CubicSplineInterpolation(xs, A)
g = zeros(1)
Interpolations.gradient!(g,interp_cubic,1.0)

you get

julia> Interpolations.gradient!(g,interp_cubic,1.0)
1-element Vector{Float64}:
 1.0001435205929399

julia> g
1-element Vector{Float64}:
 1.0001435205929399

so far so good.

Then I tried

g = zeros(2)
Interpolations.gradient!(g,interp_cubic,[1.0,1.05])

and I get

julia> Interpolations.gradient!(g,interp_cubic,[1.0,1.05])
ERROR: StackOverflowError:
Stacktrace:
 [1] gradient!(dest::Vector{Float64}, itp::Interpolations.Extrapolation{Float64, 1, ScaledInterpolation{Float64, 1, Interpolations.BSplineInterpolation{Float64, 1, OffsetArrays.OffsetVector{Float64, Vector{Float64}}, BSpline{Cubic{Line{OnGrid}}}, Tuple{Base.OneTo{Int64}}}, BSpline{Cubic{Line{OnGrid}}}, Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}}, BSpline{Cubic{Line{OnGrid}}}, Throw{Nothing}}, x::Vector{Float64}) (repeats 79984 times)
   @ Interpolations ~/.julia/packages/Interpolations/QLbbu/src/Interpolations.jl:423

any chance this is a bug?

[mkitti]

Yes, it looks like the fallback version of this method ended falling back to itself, repeatedly:

Interpolations.jl/src/Interpolations.jl

Lines 422 to 424 in 89ae91e

	function gradient!(dest, itp::AbstractInterpolation, x::Vararg{UnexpandedIndexTypes})
	gradient!(dest, itp, to_indices(itp, x)...)
	end

How about this for a workaround?

julia> g = zeros(2)
2-element Vector{Float64}:
 0.0
 0.0

julia> x = [1.0, 1.05]
2-element Vector{Float64}:
 1.0
 1.05

julia> for i in eachindex(g)
           Interpolations.gradient!(@view(g[i:i]), interp_cubic, x[i])
       end

julia> g
2-element Vector{Float64}:
 1.0001435205929399
 0.9523565763011438

[mkitti]

Here's a more compact version:

julia> g = zeros(2)
2-element Vector{Float64}:
 0.0
 0.0

julia> Interpolations.gradient!.(Iterators.partition(g, 1), Ref(interp_cubic), [1.0, 1.05])
2-element Vector{SubArray{Float64, 1, Vector{Float64}, Tuple{UnitRange{Int64}}, true}}:
 [1.0001435205929399]
 [0.9523565763011438]

julia> g
2-element Vector{Float64}:
 1.0001435205929399
 0.9523565763011438

[gianmariomanca]

They run but they are slower for me than the example mentioned on discourse:

https://discourse.julialang.org/t/multi-therading-and-dierckx-jl-interpolations-jl-and-gradients/67452/2

But for the gradient I have to do something like:

const jac_model_y(x) = [-Interpolations.gradient(interp_model_y,xx)[1] for xx in (xg .- x)]
(where xg is a constant vector and x is a scalar shift)
because Interpolations.gradient returns a 1-element SVector.

[timholy]

One could prevent the StackOverflowError with

y = to_indices(itp, x)
y isa Tuple{Vararg{UnexpandedIndices}} && error(...)
gradient!(dest, itp, y...) 

StackOverflowErrors are nasty: sometimes they take seconds or minutes to resolve, and depending on platform can corrupt your running session. I often quit immediately after getting a StackOverflowError and start a fresh session. Not true for other error types, of course.