Interpolations.jl
Interpolations.jl copied to clipboard
Canonical way of doing concrete type annotations
I'm trying to use Interpolations as a replacement for vectors of the equivalent y-interpolation points in a struct. That is, for a function y(x)
discretized at points xr
, yr
, instead of storing xr::Vector{Float64}, yr::Vector{Float64}
and accessing them separately, I'd like to be able to interpolate (xr, yr) at the struct instantiation, store the resulting Extrapolation, and later ask for y(x)
.
However, to do this efficiently, I would want a concrete type annotation, i.e., y::Extrapolation{T,N,ITPT,IT,ET}
, and in practice this annotation can get impractically long; for the case of two Float64 arrays, it's
julia> LinearInterpolation(xr, yr, extrapolation_bc=Flat()) |> typeof
Interpolations.Extrapolation{Float64, 1, Interpolations.GriddedInterpolation{Float64, 1, Float64, Gridded{Linear{Throw{OnGrid}}}, Tuple{Vector{Float64}}}, Gridded{Linear{Throw{OnGrid}}}, Flat{Nothing}}
and when I add in Unitful.jl values, it gets even longer. It isn't practical to annotate each interpolation with one of these (and it would break abstraction barriers), so is there a canonical/recommended way to make this concrete?
Alternatively, is there a more Julian way of using an interpolation as a drop-in substitution for a vector without remaking the interpolation many times or using an abstract field? Adding these long annotations makes me feel like there's an easier implementation, but I can't find it or figure it out.
Thank you!
You could just leave the type of the field as a parameter:
julia> struct Container{T}
a::T
b::T
end
julia> using Interpolations
julia> A = LinearInterpolation(xr, yr, extrapolation_bc=Flat());
julia> B = LinearInterpolation(sort(yr), xr, extrapolation_bc=Flat());
julia> C = Container(A,B)
Container{Interpolations.Extrapolation{Float64, 1, Interpolations.GriddedInterpolation{Float64, 1, Float64, Gridded{Linear{Throw{OnGrid}}}, Tuple{Vector{Float64}}}, Gridded{Linear{Throw{OnGrid}}}, Flat{Nothing}}}(10-element extrapolate(interpolate((::Vector{Float64},), ::Vector{Float64}, Gridded(Linear())), Flat()) with element type Float64:
0.03892484282890618
0.4326413224354688
0.7086107470140933
0.9468036185143996
0.6129393881428883
0.006543941392038288
0.732982887955143
0.6066581318875
0.42783287190296404
0.0663913828245517, 10-element extrapolate(interpolate((::Vector{Float64},), ::Vector{Float64}, Gridded(Linear())), Flat()) with element type Float64:
0.028064095354493168
0.028917369318009478
0.045579425288025166
0.05588372727912427
0.16536653973633275
0.2467204261726983
0.3768350265988216
0.710573746547038
0.7667599456032246
0.9829913001975853)
julia> typeof(C)
Container{Interpolations.Extrapolation{Float64, 1, Interpolations.GriddedInterpolation{Float64, 1, Float64, Gridded{Linear{Throw{OnGrid}}}, Tuple{Vector{Float64}}}, Gridded{Linear{Throw{OnGrid}}}, Flat{Nothing}}}
That works when all the data are the same type, but I've got several interpolations all of different types because they're carrying Unitful.jl units (I'm dealing with an atmosphere so I'm interpolating pressure, temperature, density etc. against length). I could make that work with parametric types T1...Tn, but that'd make the container type signature awkward to interpret. Another option is to do this for Length{Float64} vs Float64 and re-add the dimensions after calling the interpolating function, which is workable but I'd like to see if there's something more elegant (stripping/re-adding units repeatedly feels against multiple dispatch).
Do you actually need the extrapolation part?
Yes, these are quantities that are usually specified through predefined data arrays at discrete sample points, but I want to treat them as continuous functions without rebuilding the function at each call.
Is the domain that you are sampling not fixed or bounded?
Users might request values beyond the range, in which case I'd like the Flat boundary condition to be used without having to specify, e.g. y(max(x, x1))
On Wed, Apr 20, 2022, 6:07 PM Mark Kittisopikul @.***> wrote:
Is the domain that you are sampling not fixed or bounded?
— Reply to this email directly, view it on GitHub https://github.com/JuliaMath/Interpolations.jl/issues/488#issuecomment-1104193250, or unsubscribe https://github.com/notifications/unsubscribe-auth/AIK476UQAJM2553HJIFQ7CTVGA2UXANCNFSM5TZ2PIUA . You are receiving this because you authored the thread.Message ID: @.***>
that'd make the container type signature awkward to interpret
Who or what is doing the interpreting? You can modify how the type is displayed for human consumption if that is what you mean.
The way to do this is to use parameters. I'm not sure if I can provide more advice than that.
Perhaps you should consider relaxing your concreteness requirement. If you are serious about the abstraction barriers, perhaps the most important aspect is what type is returned when requesting a value at some coordinate. In that case AbstractExtrapolation{T}
or AbstractInterpolation{T}
may be a more than sufficient description where T
is the expected return type.
I think my use case was a better fit to DataInterpolations.jl.