Weight computation for scaled interpolation

[DylanMMarques]

Hello,

I am working with Interpolations.jl and I am using the weight computation feature. However, the results are not what I was expecting. This code works fine if the interpolation is not scaled but it is not correct due to the scale step. Is it possible to calculate the weight for a scaled interpolation? Or does the weight calculation interface works differently than this with a scaled interpolation?

using Interpolations
x_X = range(0, 100., length = 10)
y = 1:length(x_X)
itp = interpolate(y, BSpline(Linear()))
itp_scaled = scale(itp, (x_X))

wis = Interpolations.weightedindexes((Interpolations.value_weights,), Interpolations.itpinfo(itp_scale)..., (2.5,))
y[wis[1]] # 2.5

itp_scaled(2.5) # 1.225

Thanks

[mkitti]

All scale interpolation does is scale the input:

Interpolations.jl/src/scaling/scaling.jl

Lines 102 to 108 in af35e77

	coordlookup(r::AbstractUnitRange, x) = x - first(r) + oneunit(eltype(r))
	# coordlookup(i::Bool, r::AbstractRange, x) = i ? coordlookup(r, x) : convert(typeof(coordlookup(r,x)), x)
	coordlookup(r::StepRange, x) = (x - r.start) / r.step + oneunit(eltype(r))
	coordlookup(r::AbstractRange, x) = (x - first(r)) / step(r) + oneunit(eltype(r))
	boundstep(r::AbstractRange) = 0.5*step(r)
	rescale_gradient(r::AbstractRange, g) = g / step(r)

It does not change the weights.

[DylanMMarques]

The scale function stretch the x-axis of the interpolation. Therefore, the weights are dependent on the scale because the weights are dependent on the x-axis of the interpolated function.
The output of the function:

wis = Interpolations.weightedindexes((Interpolations.value_weights,), Interpolations.itpinfo(itp_scale)..., (2.5,))
# wis = (Interpolations.WeightedAdjIndex{2, Float64}(2, (0.5, 0.5)),)

The itp_scale is:
f(0) = y[1] and f(x_X[1]) = y[2], so the weights of x = 2.5 should be between y[1] and y[2], not y[2] and y[3] as the value of wis is.

[mkitti]

Are you asking how it currently works or are you suggesting an improvement? Interpolations.weightedindexes is not meant as a public interface.

[DylanMMarques]

I am tying to confirm that I understand how it works because, to me, it doesn't seem to work properly with an itp scaled by a StepRangeLength (the function is probably not meant for that). I analysed the code with more detail and this makes sense as Interpolation.itpinfo removes any reference about the scaling part.

If you agree with this, I would ask if there is any way of computing the weighted indexes of a scaled interpolation. If not, I am suggesting an improvement.

[mkitti]

Correct. The underlying Interpolation has no idea it is being scaled. What you would want to do is implement additional methods on ScaledInterpolation directly.

[DylanMMarques]

I am interested on the weights as I need to calculate an interpolation based on a (very sparse) matrix. This matrix is constructed based on the weights. For my application, it would be good to have a public interface of the weights that would work for any interpolation. Do you feel that such interface could be useful to more people? Or do you feel that the application of such interface is too small and might not be worth?

[mkitti]

It is possible. What would be the most consistent way to calculate the scaled weights from non-scaled weights? Would you be interested in submitting a pull request to implement Interpolations.weightedindexes on ScaledInterpolation?

[DylanMMarques]

Sure, I think that something like this should work nicely:

using Interpolations
x1 = range(0, 20, length = 11)
x2 = range(0, 20, length = 11)
y = x1 .+ x2'
itp = interpolate(y, BSpline(Linear()))
itp = scale(itp, x1, x2)

function weightedindexes(arg, scaleditp::ScaledInterpolation, x::NTuple{N}) where {N}
    ind = ntuple(i-> (x[i] - first(scaleditp.ranges[i])) / step(scaleditp.ranges[i]) + 1, N)
    return Interpolations.weightedindexes(arg, Interpolations.itpinfo(scaleditp)..., ind)
end

x = (5.4, 5.4)
wis = weightedindexes((Interpolations.value_weights,), itp, x)
y[wis...] == itp(x...) 

What do you think?

[mkitti]

This roughly looks fine to me. Could you turn this into a pull request?

[DylanMMarques]

I ran a few tests for the algorithm and it does not work for Interpolations.gradient_weights. For the gradients, the weight must be divided by step(scaleditp.ranges). I tried to do this but I am not very familiar with Interpolations.jl internal interface unfortunately, so it is quite challenging to me doing that. Do you think that you could update the code to return the correct values for the Interpolations.gradient_weights?

[koehlerson]

I would be interested in this as well. My use case is that I want to find the contributing points and their weights, so I don't need gradient_weights. Maybe I find time to make a PR