# 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

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:

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.