Improved syntax for gridded/scaled interpolation

[marius311]

If I've got some gridded x and y data already sitting around, its a bit of a pain to interpolate. AFAICT the best way to do it is something like this:

x = [2,4,6,8]
y = [1,4,3,6]

itp = scale(interpolate(y,BSpline(Linear())),x[1]:(x[2]-x[1]):x[end])

Is there some way this could be shortened to something like

itp = interpolate(x,y,BSpline(Linear()))

instead?

[floswald]

well you can always do this:

julia> i=interpolate((x,),y,Gridded(Linear()))
4-element Interpolations.GriddedInterpolation{Float64,1,Float64,Interpolations.Gridded{Interpolations.Linear},Tuple{Array{Float64,1}},0}:
 -0.5
  1.0
  2.5
  4.0

however I'm with you. If the data are reguarly spaced, the first call is faster I think, so it would be good to be able to easily use it in that case together with scale. I haven't been able to use that function, in fact your example gives me:

julia> x = [2,4,6,8]
4-element Array{Int64,1}:
 2
 4
 6
 8

julia> y = [1,4,3,6]
4-element Array{Int64,1}:
 1
 4
 3
 6

julia> itp = scale(interpolate(y,BSpline(Linear())),x[1]:(x[2]-x[1]):x[end])
ERROR: MethodError: no method matching interpolate(::Array{Int64,1}, ::Interpolations.BSpline{Interpolations.Linear})
Closest candidates are:
  interpolate{IT<:Union{Interpolations.BSpline{D<:Interpolations.Degree},Interpolations.NoInterp,Tuple{Vararg{Union{Interpolations.BSpline,Interpolations.NoInterp},N}}},GT<:Union{Interpolations.GridType,Interpolations.NoInterp,Tuple{Vararg{Union{Interpolations.GridType,Interpolations.NoInterp},N}}}}(::AbstractArray{T,N}, ::IT<:Union{Interpolations.BSpline{D<:Interpolations.Degree},Interpolations.NoInterp,Tuple{Vararg{Union{Interpolations.BSpline,Interpolations.NoInterp},N}}}, ::GT<:Union{Interpolations.GridType,Interpolations.NoInterp,Tuple{Vararg{Union{Interpolations.GridType,Interpolations.NoInterp},N}}}) at /Users/florian.oswald/.julia/v0.5/Interpolations/src/b-splines/b-splines.jl:57

julia> y = [1,4,3,6.0]
4-element Array{Float64,1}:
 1.0
 4.0
 3.0
 6.0

julia> x = [2,4,6,8.0]
4-element Array{Float64,1}:
 2.0
 4.0
 6.0
 8.0

julia> itp = scale(interpolate(y,BSpline(Linear())),x[1]:(x[2]-x[1]):x[end])
ERROR: MethodError: no method matching interpolate(::Array{Float64,1}, ::Interpolations.BSpline{Interpolations.Linear})

[sglyon]

Hey @floswald when using the linear splines you are forgetting to specify OnGrid or OnCell.

This should work

julia> using Interpolations

julia> x = [2,4,6,8]
4-element Array{Int64,1}:
 2
 4
 6
 8

julia> y = [1,4,3,6]
4-element Array{Int64,1}:
 1
 4
 3
 6

julia> itp = scale(interpolate(y,BSpline(Linear()), OnGrid()),x[1]:(x[2]-x[1]):x[end])
4-element Interpolations.ScaledInterpolation{Float64,1,Interpolations.BSplineInterpolation{Float64,1,Array{Float64,1},Interpolations.BSpline{Interpolations.Linear},Interpolations.OnGrid,0},Interpolations.BSpline{Interpolations.Linear},Interpolations.OnGrid,Tuple{StepRange{Int64,Int64}}}:
 -0.5
  1.0
  2.5
  4.0

[floswald]

Oh I see! I just copied the OP's code to be honest. Thanks though!

On Saturday, 15 October 2016, Spencer Lyon notifications@github.com wrote:

Hey @floswald https://github.com/floswald when using the linear splines
you are forgetting to specify OnGrid or OnCell.

This should work

julia> using Interpolations

julia> x = [2,4,6,8]
4-element Array{Int64,1}:
2
4
6
8

julia> y = [1,4,3,6]
4-element Array{Int64,1}:
1
4
3
6

julia> itp = scale(interpolate(y,BSpline(Linear()), OnGrid()),x[1]:(x[2]-x[1]):x[end])
4-element Interpolations.ScaledInterpolation{Float64,1,Interpolations.BSplineInterpolation{Float64,1,Array{Float64,1},Interpolations.BSpline{Interpolations.Linear},Interpolations.OnGrid,0},Interpolations.BSpline{Interpolations.Linear},Interpolations.OnGrid,Tuple{StepRange{Int64,Int64}}}:
-0.5
1.0
2.5
4.0

—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub
#124 (comment),
or mute the thread
https://github.com/notifications/unsubscribe-auth/AA-WdjWhL11KBNcG37rMMFkz8nGIAcZnks5q0UqCgaJpZM4KOLD-
.

[floswald]

hey @spencerlyon2 while you are here. can you sort me out with this thing? I just can't understand the result I'm getting after rescaling. Or I don't understand how to use it.

using Interpolations
v = [x^2 for x in 0:0.1:1]
itp=interpolate(v,BSpline(Linear()),OnGrid())
itp[1]
# 0.0
itp[11]
# 1.0
scale(itp,0:0.1:1)
itp[0]
# -0.010000000000000002
# why is this not equal to 0.0, i.e.  the value at the lowest index?

[sglyon]

I believe the issue is that the scale function doesn't alter the original itp object.

Try this:

julia> typeof(itp)
Interpolations.BSplineInterpolation{Float64,1,Array{Float64,1},Interpolations.BSpline{Interpolations.Linear},Interpolations.OnGrid,0}

julia> sitp = scale(itp,0:0.1:1);

julia> typeof(sitp)
Interpolations.ScaledInterpolation{Float64,1,Interpolations.BSplineInterpolation{Float64,1,Array{Float64,1},Interpolations.BSpline{Interpolations.Linear},Interpolations.OnGrid,0},Interpolations.BSpline{Interpolations.Linear},Interpolations.OnGrid,Tuple{FloatRange{Float64}}}

julia> itp[1]
0.0

julia> sitp[0.0]
0.0

[floswald]

Dang that was stupid! I don't know why but it felt so natural that it would
scale the original object. But has no ! there!

Thanks!

On Sunday, 16 October 2016, Spencer Lyon notifications@github.com wrote:

I believe the issue is that the scale function doesn't alter the original
itp object.

Try this:

julia> typeof(itp)
Interpolations.BSplineInterpolation{Float64,1,Array{Float64,1},Interpolations.BSpline{Interpolations.Linear},Interpolations.OnGrid,0}

julia> sitp = scale(itp,0:0.1:1);

julia> typeof(sitp)
Interpolations.ScaledInterpolation{Float64,1,Interpolations.BSplineInterpolation{Float64,1,Array{Float64,1},Interpolations.BSpline{Interpolations.Linear},Interpolations.OnGrid,0},Interpolations.BSpline{Interpolations.Linear},Interpolations.OnGrid,Tuple{FloatRange{Float64}}}

julia> itp[1]
0.0

julia> sitp[0.0]
0.0

—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub
#124 (comment),
or mute the thread
https://github.com/notifications/unsubscribe-auth/AA-WdoAFwNkgQF2J0lOgPwFIU2rB81okks5q0WrQgaJpZM4KOLD-
.

[fredRos]

Beyond supporting the request for simpler syntax, I wish it were more clearly stated in the readme how to handle interpolations without unit spacing between elements. For a function y(x), scale(interpolate(y, BSpline(Cubic(Line())), OnGrid()), x) only works for x::Range but not for x::Vector{Float64} which is a nuisance.

[sglyon]

The BSpline variants in this package all make the assumption that the knot vector (x points in y(x)) is equal to 1:length(x) in each dimension. This is built into the package as an assumption driving the derivation of the underlying mathematics.

scale allows you to shift the x points to be any evenly spaced domain. Calling scale on an interpolation object adds one extra operation before the actual interpolation happens: it remaps from the chosen x domain into 1:length(x) and hands off to the underlying routines I mentioned above. The reason we need it to be evenly spaced is that we need this mapping to be an obvious bijection so it is clear where on the original 1:length(x) domain we fall.

In order to support arbitrary arrays for x we would need to re-derive the mathematics and implement them. That is what the Gridded (not BSpline) family of interpolation types is for. So far we have only done this for Constant, NoInterp, and Linear. I have some math worked out for doing this for Cubic(Line()), but there are a couple conceptual issues holding me back from implementing it. I'm not sure when I'll have the time to do it.

For you your example the best we can currently do is

interpolate((x,), y, Gridded(Linear()))

[fredRos]

Thanks for the quick reply

On 25.04.2017 16:09, Spencer Lyon wrote: The |BSpline| variants in this package all make the assumption that the knot vector (|x| points in |y(x)|) is equal to |1:length(x)| in each dimension. This is built into the package as an assumption driving the derivation of the underlying mathematics.

Yes, I understand but I missed this clear statement in the readme. Perhaps I overlooked it. If not, please consider adding it.

|scale| allows you to shift the |x| points to be any /evenly spaced/ domain. Calling |scale| on an interpolation object adds one extra operation before the actual interpolation happens: it remaps from the chosen |x| domain into |1:length(x)| and hands off to the underlying routines I mentioned above. The reason we need it to be evenly spaced is that we need this mapping to be an obvious bijection so it is clear where on the original |1:length(x)| domain we fall.

Again, the fact that you reimplemented scale within Interpolations is not mentioned in the documentation. I just found it by chance at https://lectures.quantecon.org/jl/julia_libraries.html by googling around I understand you need the bijection. Why does it not work for an array of strictly ascending values? Or that's just what's in `Gridded`?

In order to support arbitrary arrays for |x| we would need to re-derive the mathematics and implement them. That is what the |Gridded| (not |BSpline|) family of interpolation types is for. So far we have only done this for |Constant|, |NoInterp|, and |Linear|. I have some math worked out for doing this for |Cubic(Line())|, but there are a couple conceptual issues holding me back from implementing it. I'm not sure when I'll have the time to do it. For you your example the best we can currently do is interpolate((x,), y, Gridded(Linear()))

In fact in my example the array consists of regularly spaced points, so I can stick with the `scale(...)` approach. But thanks for the clarification.

…

— You are receiving this because you commented. Reply to this email directly, view it on GitHub <#124 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/ABv63VM03O_xhRn-Y8eM5miNvm3oH5vtks5rzf6ogaJpZM4KOLD->.

[sglyon]

Hey @fredRos we should add a note to the readme about scale. I'll open an issue as a reminder.

The reason scale only works for evenly spaced range-like objects is that we leverage the uniform step size when we do the mapping back into 1:length(x).

The Gridded family of methods is meant to handle exactly the case you mentioned: arbitrary weakly increasing arrays for x