Broadcasting throws StackOverflowError

[miguelbiron]

First of all, thank you for such a great package! Now, I wanted to report an issue that arose when updating to the latest release of Interpolations. Basically, dot-broadcasting is broken -- at least for SteffenMonotonicInterpolation. Simple example:

using Interpolations, Plots

percentile_values = [0.0, 0.01, 0.1, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 91.0, 92.0, 93.0, 94.0, 95.0, 96.0, 97.0, 98.0, 99.0, 99.9, 99.99, 100.0]
y = sort(randn(length(percentile_values)))
itp_cdf = extrapolate(interpolate(y, percentile_values, SteffenMonotonicInterpolation()), Flat());
plot(itp_cdf,-3,3)                 # StackOverflowError
t = -3.0:0.01:3.0 
itp_cdf.(t)                        # StackOverflowError
interpolated_cdf = map(itp_cdf, t) # works

The full trace of the error I get in the two cases is

ERROR: StackOverflowError:
Stacktrace:
 [1] Base.Broadcast.BroadcastStyle(#unused#::Type{Interpolations.MonotonicInterpolation{Float64, Float64, Float64, Float64, Float64, SteffenMonotonicInterpolation, Vector{Float64}, Vector{Float64}}}) (repeats 79984 times)
   @ Interpolations ~/.julia/packages/Interpolations/UAW5r/src/gpu_support.jl:68

Version information

julia> versioninfo()
Julia Version 1.7.3
Commit 742b9abb4d (2022-05-06 12:58 UTC)
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: Intel(R) Core(TM) i7-6820HQ CPU @ 2.70GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-12.0.1 (ORCJIT, skylake)
Environment:
  JULIA_EDITOR = code
  JULIA_NUM_THREADS = 4

[mkitti]

What version of Interpolations.jl were you von and what version are you on now?

[mkitti]

@N5N3 this results in an infinite recursion:

https://github.com/JuliaMath/Interpolations.jl/blame/master/src/gpu_support.jl#L68

[mkitti]

This might fix the problem. I will test it thoroughly later this week.

import Interpolations: MonotonicInterpolation, root_storage_type
Interpolations.root_storage_type(::Type{T}) where {T<:MonotonicInterpolation} = root_storage_type(fieldtype(T, 2))

julia> itp_cdf.(t)
601-element Vector{Float64}:
   0.0
   0.0
   0.0
   0.0
   0.0
   0.0
   0.0
   ⋮
 100.0
 100.0
 100.0
 100.0
 100.0
 100.0

[N5N3]

Opps. root_storage_type was added only because BroadcastStyle for OffsetArray always returns a DefaultArrayStyle. (And of cource we should not do type piratcy here) Edit: we should just add a default fallback for AbstractInterpolation.

[mkitti]

@N5N3 let me know if you can send a pull request. Otherwise, I will get to this towards the end of the week.