JuliaMath / JuliaMath/Interpolations.jl
Vector-valued Untiful Support
- Dominant language
- Julia
- Stars
- 575
- Forks
- 117
- PR merge metrics
- No merged PRs in 30d
Description
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
```
```julia
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
```julia
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
```julia
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
```julia
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.
Contributor guide
No contributing guide indexed for this repository
Assessment
This issue has not been assessed yet.