Interpolating into a series of matrices

I’m writing a package that needs to generically interpolate a matrix from a time-series of matrices.

For example, maybe it’s a series of spatial matrices of monthly mean values, that I would like to interpolate to a single matrix of values for a particular day somewhere in the series.

Is there an options in an interpolations package that would handle this kind of thing, giving the options of various splines or linear interpolation?

I also want to run it on a GPU (greedy…), so broadcast operations are good here, including ways I can roll my own. But hopefully without having to define all the spline etc. methods myself.

2 Likes

AFAIK Interpolations.jl should “just work”:

julia> data = [[1 2; 3 4], [5 6; 7 8]]
2-element Array{Array{Int64,2},1}:
 [1 2; 3 4]
 [5 6; 7 8]

julia> itp = interpolate(data, BSpline(Linear()))
2-element interpolate(::Array{Array{Int64,2},1}, BSpline(Linear())) with element type Array{Float64,2}:
 [1 2; 3 4]
 [5 6; 7 8]

julia> itp[1.0]
2×2 Array{Float64,2}:
 1.0  2.0
 3.0  4.0

julia> itp[1.1]
2×2 Array{Float64,2}:
 1.4  2.4
 3.4  4.4
2 Likes

Although apparently not for splines:

julia> itp = interpolate(data, BSpline(Cubic(Line(OnGrid()))))
ERROR: MethodError: no method matching zero(::Type{Array{Int64,2}})
Closest candidates are:
  zero(::Type{Pkg.Resolve.FieldValue}) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.5/Pkg/src/Resolve/fieldvalues.jl:38
  zero(::Type{DateTime}) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.5/Dates/src/types.jl:404
  zero(::Type{Measures.Length{:mm,Float64}}) at /home/raf/.julia/packages/Plots/M1wcx/src/layouts.jl:12
  ...
Stacktrace:
 [1] copy_with_padding(::Type{Array{Int64,2}}, ::Array{Array{Int64,2},1}, ::BSpline{Cubic{Line{OnGrid}}}) at /home/raf/.julia/packages/Interpolations/TBuvH/src/b-splines/prefiltering.jl:27
 [2] prefilter(::Type{Float64}, ::Type{Array{Int64,2}}, ::Array{Array{Int64,2},1}, ::BSpline{Cubic{Line{OnGrid}}}) at /home/raf/.julia/packages/Interpolations/TBuvH/src/b-splines/prefiltering.jl:39
 [3] interpolate(::Type{Float64}, ::Type{Array{Int64,2}}, ::Array{Array{Int64,
2},1}, ::BSpline{Cubic{Line{OnGrid}}}) at /home/raf/.julia/packages/Interpolations/TBuvH/src/b-splines/b-splines.jl:160
 [4] interpolate(::Array{Array{Int64,2},1}, ::BSpline{Cubic{Line{OnGrid}}}) at /home/raf/.julia/packages/Interpolations/TBuvH/src/b-splines/b-splines.jl:179
 [5] top-level scope at REPL[40]:1

Ah, interesting. Depending on the size of your array, you could use SMatrix from StaticArrays, which does have a zero defined:

julia> data = [@SMatrix[1.0 2; 3 4], @SMatrix[3.0 4; 5 6]]
2-element Array{SArray{Tuple{2,2},Float64,2,4},1}:
 [1.0 2.0; 3.0 4.0]
 [3.0 4.0; 5.0 6.0]

julia> itp = interpolate(data, BSpline(Cubic(Line(OnGrid()))))
2-element interpolate(OffsetArray(::Array{SArray{Tuple{2,2},Float64,2,4},1}, 0:3), BSpline(Cubic(Line(OnGrid())))) with element type SArray{Tuple{2,2},Float64,2,4}:
 [0.9999999999999999 2.0; 3.0 4.0]
 [3.0 4.0; 5.0 6.0]

julia> itp[1.1]
2×2 SArray{Tuple{2,2},Float64,2,4} with indices SOneTo(2)×SOneTo(2):
 1.2  2.2
 3.2  4.2

Alternatively, if you know the particular sizes of your data but don’t want to use StaticArrays, you could create a wrapper struct around a Matrix{T} and define the relevant methods (probably just zero and basic arithmetic). I did that in the past when I needed to interpolate Vector{Float64}s and it worked fine.

1 Like

2000 * 4000 :sweat_smile:

Strange that Array doesn’t have zero defined.

They will always be wrapped in GeoArray from GeoData.jl. I can just define zero on that. Thanks

Array doesn’t have zero because there isn’t a way to know how big the zeroed array should be (since zero gets passed the type, not the value). I think SizedArray from StaticArrays should work though, and without the compile time issues for large matrices.

edit: to be more clear, there’s a method for zero that takes a type, and a method that takes a value, so zero(rand(2,2)) works, but zero(typeof(rand(2,2))) does not, and it seems like Interpolations uses the latter one.

3 Likes

That makes more sense. So really this is probably a bug/oversight in Interpolations.jl using zero on the type instead of the array.

Namely here:

https://github.com/JuliaMath/Interpolations.jl/blob/eb5690069a8c5cbe97b6ac2e44f58ea1567bee8d/src/b-splines/prefiltering.jl#L19-L31

It’s using TC instead of A on line 27, but probably because of things happening below in prefilter

Wish it had a maintainer to clarify that.

It’s also inefficient in that you can’t pass in the padded_similar array so it reuses the same one instead of allocating a new one every time.

Perhaps open an issue. (How do you expect the maintainers to know about it otherwise?)

1 Like

It’s looking for a maintainer currently - but looking more closely they are trying to help with PRs to encourage someone to take it over.

In which case opening an issue is even more important, since scattered questions like this would just get lost in the transition.