Vector-valued Untiful Support

[agerlach]

Interpolation of unitful types works for ℝ -> ℝ but fails for ℝ -> ℝⁿ

using Interpolations, Unitful

x = collect(1:10)
y = [[i, 2i] for i in x]
yunit = [[i*u"m/s", 2i*u"m/s^2"]  for i in x]

interp = LinearInterpolation(x*u"s",yunit)         #ℝ -> ℝ² w/ units, Error
interp = LinearInterpolation(x*u"s",x*u"m/s")      #ℝ -> ℝ  w/ units, OK
interp = LinearInterpolation(x*u"s",y)             #ℝ -> ℝ² w/o units, OK

ERROR: LoadError: MethodError: no method matching Interpolations.GriddedInterpolation(::Type{Quantity{Int64, D, U} where {D, U}}, ::Tuple{Vector{Quantity{Int64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}}}, ::Vector{Vector{Quantity{Int64, D, U} where {D, U}}}, ::Gridded{Linear{Throw{OnGrid}}})
Closest candidates are:
  Interpolations.GriddedInterpolation(::Type{TWeights}, ::Tuple{Vararg{Union{AbstractVector{T}, Tuple} where T, N}}, ::AbstractArray{TCoefs, N}, ::IT) where {N, TCoefs, TWeights<:Real, IT<:Union{NoInterp, Tuple{Vararg{Union{NoInterp, Gridded}, N} where N}, Gridded}, pad} at /Users/me/.julia/dev/Interpolations/src/gridded/gridded.jl:37
  Interpolations.GriddedInterpolation(::Type{TWeights}, ::Tuple{Vararg{AbstractUnitRange, N}}, ::AbstractArray{TCoefs, N}, ::IT) where {N, TCoefs, TWeights<:Real, IT<:Union{NoInterp, Tuple{Vararg{Union{NoInterp, Gridded}, N} where N}, Gridded}, pad} at /Users/me/.julia/dev/Interpolations/src/gridded/gridded.jl:60
Stacktrace:
 [1] interpolate(#unused#::Type{Quantity{Int64, D, U} where {D, U}}, #unused#::Type{Vector{Quantity{Int64, D, U} where {D, U}}}, knots::Tuple{Vector{Quantity{Int64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}}}, A::Vector{Vector{Quantity{Int64, D, U} where {D, U}}}, it::Gridded{Linear{Throw{OnGrid}}})
   @ Interpolations ~/.julia/dev/Interpolations/src/gridded/gridded.jl:135
 [2] interpolate(knots::Tuple{Vector{Quantity{Int64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}}}, A::Vector{Vector{Quantity{Int64, D, U} where {D, U}}}, it::Gridded{Linear{Throw{OnGrid}}})
   @ Interpolations ~/.julia/dev/Interpolations/src/gridded/gridded.jl:152
 [3] LinearInterpolation(range::Vector{Quantity{Int64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}}, vs::Vector{Vector{Quantity{Int64, D, U} where {D, U}}}; extrapolation_bc::Throw{Nothing})
   @ Interpolations ~/.julia/dev/Interpolations/src/convenience-constructors.jl:9
 [4] LinearInterpolation(range::Vector{Quantity{Int64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}}, vs::Vector{Vector{Quantity{Int64, D, U} where {D, U}}})
   @ Interpolations ~/.julia/dev/Interpolations/src/convenience-constructors.jl:9
 [5] top-level scope
   @ ~/Library/Application Support/Code/User/globalStorage/buenon.scratchpads/scratchpads/34cdd249d862f9176704eaed2c12df8d/scratch13.jl:54
in expression starting at /Users/me/Library/Application Support/Code/User/globalStorage/buenon.scratchpads/scratchpads/34cdd249d862f9176704eaed2c12df8d/scratch13.jl:54

GriddedInterpolation requires Type{TWeights}<:Real but in this case it is an Unitful.Quanity which is a Number but not a Real. For the 1D case it looks like there is a special disbatch https://github.com/JuliaMath/Interpolations.jl/blob/master/src/requires/unitful.jl. However, there is none for the AbstractArray{AbstractArray{}} case.

I was going to put together a PR, but I see a couple paths forward and may have stumbled on some inconsistencies w tweights that I wanted to discuss first.

Possible Solutions I see:

1. add an addition disbatch. The current Quantity disbatch reinterprets the vector via ustrip. It is not clear to me how this could be handled cleanly in the array of arrays case w/o looping through each element and creating unnecessary overhead.
2. Loosening the type limit in GriddedInterpolation from <:Real to <:Number
3. Add a non-Qauntity-specific fallback for tweight, e.g. this currently stands as

tweight(A::AbstractArray) = Float64
tweight(A::AbstractArray{T}) where T<:AbstractFloat = T
tweight(A::AbstractArray{<:AbstractVector{T}}) where {T} = T
tweight(A::AbstractArray{Rational{Int}}) = Rational{Int}
tweight(A::AbstractArray{T}) where {T<:Integer} = typeof(float(zero(T)))

The first line provides a fallback for non-AbstractFloat eltypes. Allowing a similar fallback for Array{Array{}} solves this issue and passes all tests

tweight(A::AbstractArray) = Float64
tweight(A::AbstractArray{<:AbstractVector}) = Float64
tweight(A::AbstractArray{T}) where T<:AbstractFloat = T
tweight(A::AbstractArray{<:AbstractVector{T}}) where {T<:AbstractFloat} = T
tweight(A::AbstractArray{Rational{Int}}) = Rational{Int}
tweight(A::AbstractArray{T}) where {T<:Integer} = typeof(float(zero(T)))

However, it is unclear to me on what the expected behavior for tweight should be. E.g. using the current release you get

x = collect(1:10).  # Vector{Int64}
Interpolations.tweight(x) #Float64
Interpolations.tweight([[i,i] for i in x]). #Int64

w/ option 3 above the last line would now return Float64. Return Float64 seems more consistent w/ the other code to me.

[mkitti]

Another direction I am contemplating is whether we need a new package dedicated to handling the composition of Interpolations and Unitful. Perhaps it could be called UnitfulInterpolations.jl.

The reason for this is that it would be convenient to write code for Interpolations that is aware of Unitful rather than non-specific changes.

[agerlach]

What advantages would a package that is explicitly aware of Unitful have?

[mileslucas]

@mkitti I think if we need to get to the point where an extra package, UnitfulInterpolations.jl, has to exist to enable the functionality we want, then somewhere along the way we've lost the composability of Julia. This code-base is quite spread out with a ton of tiny one-liners, which makes it really hard to follow the type-info personally, but perhaps there are tricks that can be played downstream.

For example, in https://github.com/francescoalemanno/Trapz.jl/pull/11 I avoided changing any of the dispatch and did not have to depend on Unitful to enable full support of AbstractQuantitys through using oneunit.

Sorry this isn't concrete code suggestions, I just think there's already a pretty high barrier to entry for Interpolations.jl (I personally can't use it without pulling up docs everytime, and I've contributed here before 😅 ) and forcing another package to be downloaded and used is just a no-go, especially if I'm trying to use Interpolations.jl generically in my own package code.

[mkitti]

The question is if this can be implemented via Requires or if there needs to be a deeper integration than is possible via https://github.com/JuliaMath/Interpolations.jl/blob/master/src/requires/unitful.jl

One thing I am concerned about is if the needs of autodiff and unitful are fully compatible here. I've been spot patching return type issues, but a comprehensive review is needed to ensure consistency. This may be the source of the inconsistencies noted above is.

Anyways, this issue has been hanging around for a while, and I'm still not sure if I good answer for you in the abstract. @agerlach , why don't put this in the form a PR which would make it easier for me to see how this change composes with other packages.