Re-enable commented-out tests and fix issues

Author: timholy

src/Interpolations.jl

@@ -20,6 +20,7 @@ export
     Natural,
     InPlace,
     InPlaceQ,
+    Throw,
 
     LinearInterpolation,
     CubicSplineInterpolation
@@ -199,7 +200,50 @@ weights(wi::WeightedIndex) = wi.weights
 indexes(wi::WeightedAdjIndex) = wi.istart
 indexes(wi::WeightedArbIndex) = wi.indexes
 
+# Make them iterable just like numbers are
+Base.iterate(x::WeightedIndex) = (x, nothing)
+Base.iterate(x::WeightedIndex, ::Any) = nothing
+Base.isempty(x::WeightedIndex) = false
+Base.length(x::WeightedIndex) = 1
+
+### Indexing with WeightedIndex
+
+# We inject indexing with `WeightedIndex` at a non-exported point in the dispatch heirarchy.
+# This is to avoid ambiguities with methods that specialize on the array type rather than
+# the index type.
+Base.to_indices(A, I::Tuple{Vararg{Union{Int,WeightedIndex}}}) = I
+@propagate_inbounds Base._getindex(::IndexLinear, A::AbstractVector, i::Int) = getindex(A, i)  # ambiguity resolution
+@inline function Base._getindex(::IndexStyle, A::AbstractArray{T,N}, I::Vararg{Union{Int,WeightedIndex},N}) where {T,N}
+    interp_getindex(A, I)
+end
 
+# This follows a "move processed indexes to the back" strategy, so J contains the yet-to-be-processed
+# indexes and I all the processed indexes.
+interp_getindex(A::AbstractArray{T,N}, J::Tuple{Int,Vararg{Any,L}}, I::Vararg{Int,M}) where {T,N,L,M} =
+    interp_getindex(A, Base.tail(J), I..., J[1])
+function interp_getindex(A::AbstractArray{T,N}, J::Tuple{WeightedIndex,Vararg{Any,L}}, I::Vararg{Int,M}) where {T,N,L,M}
+    wi = J[1]
+    _interp_getindex(A, indexes(wi), weights(wi), Base.tail(J), I...)
+end
+interp_getindex(A::AbstractArray{T,N}, ::Tuple{}, I::Vararg{Int,N}) where {T,N} =   # termination
+    @inbounds A[I...]  # all bounds-checks have already happened
+
+## Handle expansion of a single dimension
+# version for WeightedAdjIndex
+@inline _interp_getindex(A, i::Int, weights::NTuple{K,Number}, rest, I::Vararg{Int,M}) where {M,K} =
+    weights[1] * interp_getindex(A, rest, I..., i) + _interp_getindex(A, i+1, Base.tail(weights), rest, I...)
+@inline _interp_getindex(A, i::Int, weights::Tuple{Number}, rest, I::Vararg{Int,M}) where M =
+    weights[1] * interp_getindex(A, rest, I..., i)
+_interp_getindex(A, i::Int, weights::Tuple{}, rest, I::Vararg{Int,M}) where M =
+    error("exhausted weights, this should never happen")  # helps inference
+
+# version for WeightedArbIndex
+@inline _interp_getindex(A, indexes::NTuple{K,Int}, weights::NTuple{K,Number}, rest, I::Vararg{Int,M}) where {M,K} =
+    weights[1] * interp_getindex(A, rest, I..., indexes[1]) + _interp_getindex(A, Base.tail(indexes), Base.tail(weights), rest, I...)
+@inline _interp_getindex(A, indexes::Tuple{Int}, weights::Tuple{Number}, rest, I::Vararg{Int,M}) where M =
+    weights[1] * interp_getindex(A, rest, I..., indexes[1])
+_interp_getindex(A, indexes::Tuple{}, weights::Tuple{}, rest, I::Vararg{Int,M}) where M =
+    error("exhausted weights and indexes, this should never happen")
 
 """
     w = value_weights(degree, δx)
@@ -268,7 +312,9 @@ Base.to_index(::AbstractInterpolation, x::Number) = x
 #     itp(to_indices(itp, x)...)
 # end
 function gradient(itp::AbstractInterpolation, x::Vararg{UnexpandedIndexTypes})
-    gradient(itp, to_indices(itp, x)...)
+    xi = to_indices(itp, x)
+    xi == x && error("gradient of $itp not supported for position $x")
+    gradient(itp, xi...)
 end
 function gradient!(dest, itp::AbstractInterpolation, x::Vararg{UnexpandedIndexTypes})
     gradient!(dest, itp, to_indices(itp, x)...)
@@ -283,12 +329,15 @@ end
 # @inline function (itp::AbstractInterpolation)(x::Vararg{ExpandedIndexTypes})
 #     itp.(Iterators.product(x...))
 # end
-function gradient(itp::AbstractInterpolation, x::Vararg{ExpandedIndexTypes})
-    map(y->gradient(itp, y), Iterators.product(x...))
-end
-function hessian(itp::AbstractInterpolation, x::Vararg{ExpandedIndexTypes})
-    map(y->hessian(itp, y), Iterators.product(x...))
-end
+# function gradient(itp::AbstractInterpolation, x::Vararg{ExpandedIndexTypes})
+#     map(y->tgradient(itp, y), Iterators.product(x...))
+# end
+# function hessian(itp::AbstractInterpolation, x::Vararg{ExpandedIndexTypes})
+#     map(y->thessian(itp, y), Iterators.product(x...))
+# end
+#
+# tgradient(itp, y) = gradient(itp, y...)
+# thessian(itp, y) = hessian(itp, y...)
 
 # getindex is supported only for Integer indices (deliberately)
 import Base: getindex

src/b-splines/indexing.jl

@@ -1,42 +1,3 @@
-### Indexing with WeightedIndex
-
-# We inject indexing with `WeightedIndex` at a non-exported point in the dispatch heirarchy.
-# This is to avoid ambiguities with methods that specialize on the array type rather than
-# the index type.
-Base.to_indices(A, I::Tuple{Vararg{Union{Int,WeightedIndex}}}) = I
-@propagate_inbounds Base._getindex(::IndexLinear, A::AbstractVector, i::Int) = getindex(A, i)  # ambiguity resolution
-@inline function Base._getindex(::IndexStyle, A::AbstractArray{T,N}, I::Vararg{Union{Int,WeightedIndex},N}) where {T,N}
-    interp_getindex(A, I)
-end
-
-# This follows a "move processed indexes to the back" strategy, so J contains the yet-to-be-processed
-# indexes and I all the processed indexes.
-interp_getindex(A::AbstractArray{T,N}, J::Tuple{Int,Vararg{Any,L}}, I::Vararg{Int,M}) where {T,N,L,M} =
-    interp_getindex(A, Base.tail(J), I..., J[1])
-function interp_getindex(A::AbstractArray{T,N}, J::Tuple{WeightedIndex,Vararg{Any,L}}, I::Vararg{Int,M}) where {T,N,L,M}
-    wi = J[1]
-    _interp_getindex(A, indexes(wi), weights(wi), Base.tail(J), I...)
-end
-interp_getindex(A::AbstractArray{T,N}, ::Tuple{}, I::Vararg{Int,N}) where {T,N} =   # termination
-    @inbounds A[I...]  # all bounds-checks have already happened
-
-# version for WeightedAdjIndex
-_interp_getindex(A, i::Int, weights::NTuple{K,Number}, rest, I::Vararg{Int,M}) where {M,K} =
-    weights[1] * interp_getindex(A, rest, I..., i) + _interp_getindex(A, i+1, Base.tail(weights), rest, I...)
-_interp_getindex(A, i::Int, weights::Tuple{Number}, rest, I::Vararg{Int,M}) where M =
-    weights[1] * interp_getindex(A, rest, I..., i)
-_interp_getindex(A, i::Int, weights::Tuple{}, rest, I::Vararg{Int,M}) where M =
-    error("exhausted weights, this should never happen")  # helps inference
-
-# version for WeightedArbIndex
-_interp_getindex(A, indexes::NTuple{K,Int}, weights::NTuple{K,Number}, rest, I::Vararg{Int,M}) where {M,K} =
-    weights[1] * interp_getindex(A, rest, I..., indexes[1]) + _interp_getindex(A, Base.tail(indexes), Base.tail(weights), rest, I...)
-_interp_getindex(A, indexes::Tuple{Int}, weights::Tuple{Number}, rest, I::Vararg{Int,M}) where M =
-    weights[1] * interp_getindex(A, rest, I..., indexes[1])
-_interp_getindex(A, indexes::Tuple{}, weights::Tuple{}, rest, I::Vararg{Int,M}) where M =
-    error("exhausted weights and indexes, this should never happen")
-
-
 ### Primary evaluation entry points (itp(x...), gradient(itp, x...), and hessian(itp, x...))
 
 itpinfo(itp) = (tcollect(itpflag, itp), axes(itp))
@@ -52,6 +13,14 @@ end
     @assert length(inds) == N
     itp(inds...)
 end
+@inline function (itp::BSplineInterpolation{T,N})(x::Vararg{Union{Number,AbstractVector},N}) where {T,N}
+    @boundscheck (checkbounds(Bool, itp, x...) || Base.throw_boundserror(itp, x))
+    itps = tcollect(itpflag, itp)
+    wis = dimension_wis(value_weights, itps, axes(itp), x)
+    coefs = coefficients(itp)
+    ret = [coefs[i...] for i in Iterators.product(wis...)]
+    reshape(ret, shape(wis...))
+end
 
 @inline function gradient(itp::BSplineInterpolation{T,N}, x::Vararg{Number,N}) where {T,N}
     @boundscheck checkbounds(Bool, itp, x...) || Base.throw_boundserror(itp, x)
@@ -71,12 +40,14 @@ end
     dest .= hessian(itp, x...)
 end
 
-checkbounds(::Type{Bool}, itp::AbstractInterpolation, x::Vararg{Number,N}) where N =
+checkbounds(::Type{Bool}, itp::AbstractInterpolation, x::Vararg{ExpandedIndexTypes,N}) where N =
     checklubounds(lbounds(itp), ubounds(itp), x)
 
 checklubounds(ls, us, xs) = _checklubounds(true, ls, us, xs)
-_checklubounds(tf::Bool, ls, us, xs) = _checklubounds(tf & (ls[1] <= xs[1] <= us[1]),
-                                                        Base.tail(ls), Base.tail(us), Base.tail(xs))
+_checklubounds(tf::Bool, ls, us, xs::Tuple{Number, Vararg{Any}}) =
+    _checklubounds(tf & (ls[1] <= xs[1] <= us[1]), Base.tail(ls), Base.tail(us), Base.tail(xs))
+_checklubounds(tf::Bool, ls, us, xs::Tuple{AbstractVector, Vararg{Any}}) =
+    _checklubounds(tf & all(ls[1] .<= xs[1] .<= us[1]), Base.tail(ls), Base.tail(us), Base.tail(xs))
 _checklubounds(tf::Bool, ::Tuple{}, ::Tuple{}, ::Tuple{}) = tf
 
 # Leftovers from AbstractInterpolation

src/extrapolation/extrapolation.jl

@@ -39,7 +39,29 @@ count_interp_dims(::Type{<:Extrapolation{T,N,ITPT}}, n) where {T,N,ITPT} = count
     itp = parent(etp)
     eflag = etpflag(etp)
     xs = inbounds_position(eflag, bounds(itp), x, etp, x)
-    extrapolate_value(eflag, x, xs, Tuple(gradient(itp, xs...)), itp(xs...))
+    extrapolate_value(eflag, skip_flagged_nointerp(itp, x), skip_flagged_nointerp(itp, xs), Tuple(gradient(itp, xs...)), itp(xs...))
+end
+@inline function (etp::Extrapolation{T,N})(x::Vararg{Union{Number,AbstractVector},N}) where {T,N}
+    itp = parent(etp)
+    Tret = typeof(lispyprod(zero(T), x...))
+    ret = zeros(Tret, shape(x...))
+    for (i, y) in zip(eachindex(ret), Iterators.product(x...))
+        ret[i] = etp(y...)
+    end
+    return ret
+end
+
+@inline function gradient(etp::AbstractExtrapolation{T,N}, x::Vararg{Number,N}) where {T,N}
+    itp = parent(etp)
+    if checkbounds(Bool, itp, x...)
+        gradient(itp, x...)
+    else
+        eflag = tcollect(etpflag, etp)
+        xs = inbounds_position(eflag, bounds(itp), x, etp, x)
+        g = gradient(itp, xs...)
+        skipni = t->skip_flagged_nointerp(itp, t)
+        SVector(extrapolate_gradient.(skipni(eflag), skipni(x), skipni(xs), Tuple(g)))
+    end
 end
 
 checkbounds(::Bool, ::AbstractExtrapolation, I...) = true

src/extrapolation/filled.jl

@@ -33,6 +33,16 @@ end
     @assert length(inds) == N
     etp(inds...)
 end
+@inline function (etp::FilledExtrapolation{T,N})(x::Vararg{Union{Number,AbstractVector},N}) where {T,N}
+    itp = parent(etp)
+    Tret = typeof(lispyprod(zero(T), x...))
+    ret = fill(convert(Tret, etp.fillvalue), shape(x...))
+    axsib = inbounds(itp, x...)
+    any(isempty, axsib) && return ret
+    xib = getindex.(x, axsib)
+    getindex!(view(ret, keepvectors(axsib...)...), itp, xib...)
+    return ret
+end
 
 expand_index_resid_etp(deg, fillvalue, (l, u), x, etp::FilledExtrapolation, xN) =
     (l <= x <= u || Base.throw_boundserror(etp, xN))

src/gridded/gridded.jl

@@ -4,6 +4,12 @@ struct Gridded{D<:Degree} <: InterpolationType
     degree::D
 end
 
+function Base.show(io::IO, g::Gridded)
+    print(io, "Gridded(")
+    show(io, degree(g))
+    print(io, ')')
+end
+
 const GridIndex{T} = Union{AbstractVector{T}, Tuple}
 
 struct GriddedInterpolation{T,N,TCoefs,IT<:DimSpec{Gridded},K<:Tuple{Vararg{AbstractVector}}} <: AbstractInterpolation{T,N,IT}

src/gridded/indexing.jl

@@ -8,6 +8,15 @@ end
     itp(to_indices(itp, x)...)
 end
 
+@inline function gradient(itp::GriddedInterpolation{T,N}, x::Vararg{Number,N}) where {T,N}
+    @boundscheck (checkbounds(Bool, itp, x...) || Base.throw_boundserror(itp, x))
+    wis = weightedindexes((value_weights, gradient_weights), itpinfo(itp)..., x)
+    SVector(map(inds->coefficients(itp)[inds...], wis))
+end
+@propagate_inbounds function gradient!(dest, itp::GriddedInterpolation{T,N}, x::Vararg{Number,N}) where {T,N}
+    dest .= gradient(itp, x...)
+end
+
 itpinfo(itp::GriddedInterpolation) = (tcollect(itpflag, itp), itp.knots)
 
 weightedindex_parts(fs::F, itpflag::Gridded, ax, x) where F =
@@ -45,15 +54,23 @@ end
 function weightedindex(fs::F, deg::Degree, knotvec, x, iclamp) where F
     @inbounds l, u = knotvec[iclamp], knotvec[iclamp+1]
     δx = (x - l)/(u - l)
-    (position=iclamp, coefs=fmap(fs, deg, δx))
+    (position=iclamp, coefs=rescale_gridded(fs, fmap(fs, deg, δx), u-l))
 end
 
+rescale_gridded(fs::F, coefs, Δx) where F =
+    (rescale_gridded(fs[1], coefs[1], Δx), rescale_gridded(Base.tail(fs), Base.tail(coefs), Δx)...)
+rescale_gridded(::Tuple{}, ::Tuple{}, Δx) = ()
+rescale_gridded(::typeof(value_weights), coefs, Δx) = coefs
+rescale_gridded(::typeof(gradient_weights), coefs, Δx) = coefs./Δx
+rescale_gridded(::typeof(hessian_weights), coefs, Δx) = coefs./Δx.^2
+
 @inline function (itp::GriddedInterpolation{T,N})(x::Vararg{Union{Number,AbstractVector},N}) where {T,N}
     @boundscheck (checkbounds(Bool, itp, x...) || Base.throw_boundserror(itp, x))
     itps = tcollect(itpflag, itp)
     wis = dimension_wis(value_weights, itps, itp.knots, x)
     coefs = coefficients(itp)
-    [coefs[i...] for i in Iterators.product(wis...)]
+    ret = [coefs[i...] for i in Iterators.product(wis...)]
+    reshape(ret, shape(wis...))
 end
 
 function dimension_wis(f::F, itps, knots, xs) where F

src/io.jl

@@ -1,4 +1,4 @@
-function Base.showarg(io::IO, A::BSplineInterpolation{T,N,TW,ST,GT}, toplevel) where {T,N,TW,ST,GT}
+function Base.showarg(io::IO, A::BSplineInterpolation{T,N,TW,ST}, toplevel) where {T,N,TW,ST}
     print(io, "interpolate(")
     Base.showarg(io, A.coefs, false)
     print(io, ", ")
@@ -10,19 +10,19 @@ function Base.showarg(io::IO, A::BSplineInterpolation{T,N,TW,ST,GT}, toplevel) w
     end
 end
 
-# function Base.showarg(io::IO, A::GriddedInterpolation{T,N,TC,ST,K}, toplevel) where {T,N,TC,ST,K}
-#     print(io, "interpolate(")
-#     _showknots(io, A.knots)
-#     print(io, ", ")
-#     Base.showarg(io, A.coefs, false)
-#     print(io, ", ")
-#     _showtypeparam(io, ST)
-#     if toplevel
-#         print(io, ") with element type ",T)
-#     else
-#         print(io, ')')
-#     end
-# end
+function Base.showarg(io::IO, A::GriddedInterpolation{T,N,TC,ST,K}, toplevel) where {T,N,TC,ST,K}
+    print(io, "interpolate(")
+    _showknots(io, A.knots)
+    print(io, ", ")
+    Base.showarg(io, A.coefs, false)
+    print(io, ", ")
+    show(io, itpflag(A))
+    if toplevel
+        print(io, ") with element type ",T)
+    else
+        print(io, ')')
+    end
+end
 
 _showknots(io, A) = Base.showarg(io, A, false)
 function _showknots(io, tup::NTuple{N,Any}) where N
@@ -35,31 +35,31 @@ function _showknots(io, tup::NTuple{N,Any}) where N
     print(io, ')')
 end
 
-# function Base.showarg(io::IO, A::ScaledInterpolation{T}, toplevel) where {T}
-#     print(io, "scale(")
-#     Base.showarg(io, A.itp, false)
-#     print(io, ", ", A.ranges, ')')
-#     if toplevel
-#         print(io, " with element type ",T)
-#     end
-# end
+function Base.showarg(io::IO, A::ScaledInterpolation{T}, toplevel) where {T}
+    print(io, "scale(")
+    Base.showarg(io, A.itp, false)
+    print(io, ", ", A.ranges, ')')
+    if toplevel
+        print(io, " with element type ",T)
+    end
+end
 
-# function Base.showarg(io::IO, A::Extrapolation{T,N,TI,IT,GT,ET}, toplevel) where {T,N,TI,IT,GT,ET}
-#     print(io, "extrapolate(")
-#     Base.showarg(io, A.itp, false)
-#     print(io, ", ")
-#     _showtypeparam(io, ET)
-#     print(io, ')')
-#     if toplevel
-#         print(io, " with element type ",T)
-#     end
-# end
+function Base.showarg(io::IO, A::Extrapolation{T,N,TI,IT,ET}, toplevel) where {T,N,TI,IT,ET}
+    print(io, "extrapolate(")
+    Base.showarg(io, A.itp, false)
+    print(io, ", ")
+    show(io, etpflag(A))
+    print(io, ')')
+    if toplevel
+        print(io, " with element type ",T)
+    end
+end
 
-# function Base.showarg(io::IO, A::FilledExtrapolation{T,N,TI,IT,GT}, toplevel) where {T,N,TI,IT,GT}
-#     print(io, "extrapolate(")
-#     Base.showarg(io, A.itp, false)
-#     print(io, ", ", A.fillvalue, ')')
-#     if toplevel
-#         print(io, " with element type ",T)
-#     end
-# end
+function Base.showarg(io::IO, A::FilledExtrapolation{T,N,TI,IT}, toplevel) where {T,N,TI,IT}
+    print(io, "extrapolate(")
+    Base.showarg(io, A.itp, false)
+    print(io, ", ", A.fillvalue, ')')
+    if toplevel
+        print(io, " with element type ",T)
+    end
+end