Cubic interpolation of a vector of matrices with grid/range doesn't work

[andreasvarga]

I just started using Interpolations, being interested to interpolate time values stored in a vector of array. While linear interpolation works

julia> itp = LinearInterpolation([1,2,3],[rand(2,2),rand(2,2),rand(2,2)])
3-element extrapolate(interpolate((::Vector{Int64},), ::Vector{Matrix{Float64}}, Gridded(Linear())), Throw()) with element type Matrix{Float64}:
 [0.45075441351126755 0.6948520432892077; 0.9632771824389577 0.13939170840849813]
 [0.7724086132599711 0.4960608819423973; 0.8859328713683932 0.44873824726319933]
 [0.11167858142523124 0.7060973358621785; 0.4175384295030793 0.9972636473469578]

I encountered the following error when trying to use cubic splines with the gridding option

julia> itp = interpolate([rand(2,2),rand(2,2),rand(2,2)], BSpline(Cubic(Line(OnGrid()))))
ERROR: MethodError: no method matching zero(::Type{Matrix{Float64}})
Closest candidates are:
  zero(::Union{Type{P}, P}) where P<:Dates.Period at C:\Users\Andreas\AppData\Local\Programs\Julia-1.7.1\share\julia\stdlib\v1.7\Dates\src\periods.jl:53
  zero(::SA) where SA<:StaticArrays.StaticArray at C:\Users\Andreas\.julia\packages\StaticArrays\A0XaR\src\linalg.jl:92
  zero(::OffsetArrays.OffsetArray) at C:\Users\Andreas\.julia\packages\OffsetArrays\I5Pfg\src\OffsetArrays.jl:395  ...
Stacktrace:
 [1] copy_with_padding(#unused#::Type{Matrix{Float64}}, A::Vector{Matrix{Float64}}, it::BSpline{Cubic{Line{OnGrid}}})
   @ Interpolations C:\Users\Andreas\.julia\packages\Interpolations\Glp9h\src\b-splines\prefiltering.jl:27       
 [2] prefilter(#unused#::Type{Float64}, #unused#::Type{Matrix{Float64}}, A::Vector{Matrix{Float64}}, it::BSpline{Cubic{Line{OnGrid}}})
   @ Interpolations C:\Users\Andreas\.julia\packages\Interpolations\Glp9h\src\b-splines\prefiltering.jl:39       
 [3] interpolate(#unused#::Type{Float64}, #unused#::Type{Matrix{Float64}}, A::Vector{Matrix{Float64}}, it::BSpline{Cubic{Line{OnGrid}}})
   @ Interpolations C:\Users\Andreas\.julia\packages\Interpolations\Glp9h\src\b-splines\b-splines.jl:160
 [4] interpolate(A::Vector{Matrix{Float64}}, it::BSpline{Cubic{Line{OnGrid}}})
   @ Interpolations C:\Users\Andreas\.julia\packages\Interpolations\Glp9h\src\b-splines\b-splines.jl:184
 [5] top-level scope
   @ REPL[26]:1

I hope it is my fault (using the wrong options) and I would be very grateful for any hint to overcome this error. I am using Julia 1.7.1 and Interpolations v0.13.5.

[dkarrasch]

Apparently, you can interpolate data "cubicly" only for data types that have zero defined. Probably, you can fix this by making your array SMatrix-valued (from StaticArrays.jl).

[andreasvarga]

Is this possibly related to

julia> zero([rand(3,3),rand(3,3)])
ERROR: MethodError: no method matching zero(::Type{Matrix{Float64}})
Closest candidates are:
  zero(::Union{Type{P}, P}) where P<:Dates.Period at C:\Users\Andreas\AppData\Local\Programs\Julia-1.7.1\share\julia\stdlib\v1.7\Dates\src\periods.jl:53
  zero(::SA) where SA<:StaticArrays.StaticArray at C:\Users\Andreas\.julia\packages\StaticArrays\A0XaR\src\linalg.jl:92
  zero(::RationalTransferFunction) at C:\Users\Andreas\Documents\software\Julia\DescriptorSystems.jl\src\types\RationalTransferFunction.jl:434
  ...
Stacktrace:
 [1] zero(x::Vector{Matrix{Float64}})
   @ Base .\abstractarray.jl:1133
 [2] top-level scope
   @ REPL[8]:1

which fails too ? Actually, I wonder if

julia> zero.([rand(3,3),rand(3,3)])
2-element Vector{Matrix{Float64}}:
 [0.0 0.0 0.0; 0.0 0.0 0.0; 0.0 0.0 0.0]
 [0.0 0.0 0.0; 0.0 0.0 0.0; 0.0 0.0 0.0]

produces what has to be expected.

[andreasvarga]

I managed using componentwise computation for the [i,j] entries through all matrices. However, I look forward to the next versions supporting this feature for cubic splines too (similarly as for linear interpolation).

[mkitti]

I'm starting to catch up here. I think you want to use Gridded Cubic if I understand the problem correctly.

[andreasvarga]

I would like to interpolate time data stored in a vector of matrices say A, where A[i] contains the value of a time-dependent matrix A(t) at a time moment t_i. The time moments form a uniform grid and can be alternatively specified by a range of values.

In my problem, A(t) is periodic with period T and the time range is (0:N-1)dt and Ndt = T.

The built interpolation structure itp is used to define a function t -> itp(mod(t,T)) to evaluate A(t) for arbitrary t.

The linear interpolation can be used to work directly with an extended range (0:N)*dt and an extended data set [A;A[1]]. However, this doesn't work for cubic interpolation.

Mark Kittisopikul @.***> schrieb am Do., 24. Feb. 2022, 22:50:

I'm starting to catch up here. I think you want to use Gridded Cubic if I understand the problem correctly.

— Reply to this email directly, view it on GitHub https://github.com/JuliaMath/Interpolations.jl/issues/481#issuecomment-1050299449, or unsubscribe https://github.com/notifications/unsubscribe-auth/ALJDHEBJXBNDXFBELJH7ECDU42RYVANCNFSM5PEDRMOA . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

You are receiving this because you authored the thread.Message ID: @.***>

[mkitti]

Just to be clear on the issue here, currently we have Linear() |> Gridded working.

julia> interpolate(([1,2,3],), [rand(2,2), rand(2,2), rand(2,2)], Linear() |> Gridded)
3-element interpolate((::Vector{Int64},), ::Vector{Matrix{Float64}}, Gridded(Linear())) with element type Matrix{Float64}:
 [0.9931055694014964 0.6928134801613715; 0.9929877647442734 0.053553541846786845]
 [0.2853354863085713 0.28932791645525635; 0.790812658397169 0.9976569037234465]
 [0.7379706351081073 0.3894714019525468; 0.6146898249659012 0.8106959142441317]

What we need is OnGrid() |> Line |> Cubic |> Gridded.

julia> interpolate(([1,2,3],), [rand(2,2), rand(2,2), rand(2,2)], OnGrid() |> Line |> Cubic |> Gridded)
ERROR: only Linear, Constant, and NoInterp supported, got Gridded(Cubic(Line(OnGrid())))
Stacktrace:
 [1] error(s::String)
   @ Base .\error.jl:33
 [2] check_gridded
   @ C:\Users\kittisopikulm\.julia\packages\Interpolations\Glp9h\src\gridded\gridded.jl:73 [inlined]
 [3] Interpolations.GriddedInterpolation(#unused#::Type{Float64}, knots::Tuple{Vector{Int64}}, A::Vector{Matrix{Float64}}, it::Gridded{Cubic{Line{OnGrid}}})
   @ Interpolations C:\Users\kittisopikulm\.julia\packages\Interpolations\Glp9h\src\gridded\gridded.jl:40
 [4] interpolate(#unused#::Type{Float64}, #unused#::Type{Matrix{Float64}}, knots::Tuple{Vector{Int64}}, A::Vector{Matrix{Float64}}, it::Gridded{Cubic{Line{OnGrid}}})
   @ Interpolations C:\Users\kittisopikulm\.julia\packages\Interpolations\Glp9h\src\gridded\gridded.jl:149
 [5] interpolate(knots::Tuple{Vector{Int64}}, A::Vector{Matrix{Float64}}, it::Gridded{Cubic{Line{OnGrid}}})
   @ Interpolations C:\Users\kittisopikulm\.julia\packages\Interpolations\Glp9h\src\gridded\gridded.jl:166
 [6] top-level scope
   @ REPL[34]:1

[linwaytin]

I'm wondering how difficult to finish the cubic spline for matrix data. From the mathematical point of view, it seems to be the same as the scalar case. What is the difficulty?

[timholy]

Someone who wants it just needs to invest the effort. The people who have contributed thus far to the package presumably didn't need it.

[linwaytin]

I would like to fix this issue, but I want to understand more about it first. Does anyone have some idea?