Reuse weights with extrapolation

[fredrikpaues]

I'm finding the documentation on how to reuse weights quite difficult to comprehend. But I would like to compute weights for an extrapolant/interpolant. My best attempt at replicating the code in the documentation but with my own interpolation object (below) results in a MethodError.

julia> using Interpolations
julia> x = collect(range(0.1, 0.9, 100));
julia> v = 2x;
julia> itp = LinearInterpolation(x, v; extrapolation_bc=((Flat(), Line()),));
julia> itp(0.5)  # Works just fine :)
1.0
julia> wis = Interpolations.weightedindexes((Interpolations.value_weights,), Interpolations.itpinfo(itp)..., 0.5);
ERROR: MethodError: no method matching weightedindexes(::Tuple{typeof(Interpolations.value_weights)}, ::Tuple{Gridded{Linear{Throw{OnGrid}}}}, ::Tuple{Base.OneTo{Int64}}, ::Float64)
Closest candidates are:
  weightedindexes(::F, ::Tuple{Vararg{Interpolations.Flag, N}}, ::Tuple{Vararg{AbstractVector, N}}, ::Tuple{Vararg{Number, N}}) where {F, N} at C:\...\.julia\packages\Interpolations\Glp9h\src\b-splines\indexing.jl:64

[fredrikpaues]

One step closer, the last line could be replaced with

julia> wis = Interpolations.weightedindexes((Interpolations.value_weights,), Interpolations.itpinfo(itp)..., (0.5,));

whereby

julia> v[wis...]
0.19191919191919196

which confuses me.

[mkitti]

Try

julia> wis = Interpolations.weightedindexes((Interpolations.value_weights,), Interpolations.itpinfo(itp.itp)..., (0.5,));

julia> itp.itp.coefs[wis...]
1.0

[fredrikpaues]

Thanks! But it seems that the extrapolation scheme is not respected.

julia> using Interpolations

julia> x = collect(range(0.1, 0.9, 9)); v = 2x;

julia> itp = LinearInterpolation(x, v; extrapolation_bc=((Flat(), Line()),));

julia> itp(1.0)  # Linear extrapolation
2.0

julia> wis = Interpolations.weightedindexes((Interpolations.value_weights,), Interpolations.itpinfo(itp.itp)..., (1.0,))
(Interpolations.WeightedAdjIndex{2, Float64}(8, (-1.0, 2.0)),)

julia> v[wis...]  # Linear extrapolation is respected
2.0

julia> itp.itp.coefs[wis2...]  # Linear extrapolation is respected
2.0

julia> itp(0.0)  # The flat extrapolation kicks in
0.2

julia> wis = Interpolations.weightedindexes((Interpolations.value_weights,), Interpolations.itpinfo(itp.itp)..., (0.0,))
(Interpolations.WeightedAdjIndex{2, Float64}(1, (2.0, -1.0)),)

julia> v[wis...]  # The flat extrapolation is ignored
0.0

julia> itp.itp.coefs[wis...]  # The flat extrapolation is ignored
0.0

[mkitti]

itp is the extrapolation object. itp.itp is the underlying interpolation instance. If you want the branching behavior between the two, I suggest using the higher level interface.

[fredrikpaues]

I don't quite understand what you mean by "the higher level interface". Do you mean using itp rather than itp.itp in the call to Interpolations.itpinfo?

But that doesn't do the desired interpolations either:

julia> wis = Interpolations.weightedindexes((Interpolations.value_weights,), Interpolations.itpinfo(itp)..., (1.0,))
(Interpolations.WeightedAdjIndex{2, Float64}(1, (1.0, 0.0)),)

julia> itp(1.0)
2.0

julia> v[wis...]
0.2

julia> itp.itp.coefs[wis...]
0.2

Or do you mean that I shouldn't be using Interpolations.weightedindexes at all? And instead just use the interpolant itp?

My best guess for what v[wis...] and itp.itp.coefs[wis...] output is by how much the endpoints needs to be adjusted.

[mkitti]

wis knows nothing about extrapolation. It is purely used for interpolation.

[fredrikpaues]

I see. Then I guess it is not the tool that I need. Thank you for your time! Perhaps you could that this as a feature suggestion 🙂

[eirikbrandsaas]

Cross-referencing https://github.com/JuliaMath/Interpolations.jl/issues/363.

[eirikbrandsaas]

Given that Fredrik also found the documentation on how to reuse weights difficult to comprehend and I want to understand this, I'd be happy to write something for the documentation. However, to do so there are two things I don't understand:

1. I'm not able to reuse weights for scaled interpolation, only gridded
2. Using the weights is actually slower than doing the "whole" interpolation:

using Interpolations
using BenchmarkTools

## Setup grids
xs = 0.2:0.1:0.8
A = 2xs

## Scaled 1D interpolation (re-using weights gives wrong results?)
itp_1d = LinearInterpolation((xs, ), A) 
@btime itp_1d(0.5)
wis = Interpolations.weightedindexes((Interpolations.value_weights,), Interpolations.itpinfo(itp_1d.itp)...,(0.5,))
A[wis...] # Wrong answer
itp_1d.itp.itp.coefs[wis...] # Different wrong answer (and requires another .itp)

## Gridded 1D interpolation (re-use weights works)
itp_1d = LinearInterpolation(collect(xs), A) 
@btime itp_1d(0.5)
wis = Interpolations.weightedindexes((Interpolations.value_weights,), Interpolations.itpinfo(itp_1d.itp)...,(0.75,))
@btime A[wis...] # Right answer, but slow
@btime A[wis[1]] # Right answer, but slow
@btime itp_1d.itp.coefs[wis...] # Right answer, but slow
##

[eirikbrandsaas]

There have been some changes to the package, so that A[wis...] no longer works, but one have to use Interpolations.InterpGetindex(A)[wis...] instead.

@mkitti : After rereading this issue, and other issues, (e.g., https://github.com/JuliaMath/Interpolations.jl/issues/479), is this correct: weightedindexes should only be used for gridded interpolation? For example, it should not be used for scaled interpolation nor extrapolation?

[mkitti]

I need to reorient myself here again. I'll try to respond later this week.