Implement extrapolation

Author: timholy

src/Interpolations.jl

@@ -204,9 +204,9 @@ end
 # getindex is supported only for Integer indices (deliberately)
 import Base: getindex
 @propagate_inbounds getindex(itp::AbstractInterpolation{T,N}, i::Vararg{Integer,N}) where {T,N} = itp(i...)
-@inline function getindex(itp::AbstractInterpolation{T,1}, i::Integer, j::Integer) where T
+@propagate_inbounds function getindex(itp::AbstractInterpolation{T,1}, i::Integer, j::Integer) where T
     @boundscheck (j == 1 || Base.throw_boundserror(itp, (i, j)))
-    @inbounds itp(i)
+    itp(i)
 end
 
 # deprecate getindex for other numeric indices
@@ -215,7 +215,7 @@ end
 include("nointerp/nointerp.jl")
 include("b-splines/b-splines.jl")
 # include("gridded/gridded.jl")
-# include("extrapolation/extrapolation.jl")
+include("extrapolation/extrapolation.jl")
 # include("scaling/scaling.jl")
 include("utils.jl")
 include("io.jl")

src/b-splines/b-splines.jl

@@ -53,8 +53,10 @@ axes(itp::BSplineInterpolation) = itp.parentaxes
 
 lbounds(itp::BSplineInterpolation) = _lbounds(itp.parentaxes, itpflag(itp), gridflag(itp))
 ubounds(itp::BSplineInterpolation) = _ubounds(itp.parentaxes, itpflag(itp), gridflag(itp))
-_lbounds(axs::NTuple{N,Any}, itp, gt) where N = ntuple(d->lbound(axs[d], iextract(itp,d), iextract(gt,d)), Val(N))
-_ubounds(axs::NTuple{N,Any}, itp, gt) where N = ntuple(d->ubound(axs[d], iextract(itp,d), iextract(gt,d)), Val(N))
+_lbounds(axs, itp, gt) = (lbound(axs[1], getfirst(itp), getfirst(gt)), _lbounds(Base.tail(axs), getrest(itp), getrest(gt))...)
+_ubounds(axs, itp, gt) = (ubound(axs[1], getfirst(itp), getfirst(gt)), _ubounds(Base.tail(axs), getrest(itp), getrest(gt))...)
+_lbounds(::Tuple{}, itp, gt) = ()
+_ubounds(::Tuple{}, itp, gt) = ()
 
 # The unpadded defaults
 lbound(ax, ::BSpline, ::OnCell) = first(ax) - 0.5

src/b-splines/cubic.jl

@@ -26,7 +26,7 @@ Cubic
 
 function base_rem(::Cubic, bounds, x)
     xf = floorbounds(x, bounds)
-    xf -= ifelse(xf > bounds[2]-1, oneunit(xf), zero(xf))
+    xf -= ifelse(xf > last(bounds)-1, oneunit(xf), zero(xf))
     δx = x - xf
     fast_trunc(Int, xf), δx
 end
@@ -71,9 +71,9 @@ hessian_weights(::Cubic, δx) = (1-δx, 3*δx-2, 3*(1-δx)-2, δx)
 padded_axis(ax::AbstractUnitRange, ::BSpline{<:Cubic}) = first(ax)-1:last(ax)+1
 padded_axis(ax::AbstractUnitRange, ::BSpline{Cubic{Periodic}}) = ax
 
-# Due to padding we can extend the bounds
-lbound(ax, ::BSpline{Cubic{BC}}, ::OnGrid) where BC = first(ax) - 0.5
-ubound(ax, ::BSpline{Cubic{BC}}, ::OnGrid) where BC = last(ax) + 0.5
+# # Due to padding we can extend the bounds
+# lbound(ax, ::BSpline{Cubic{BC}}, ::OnGrid) where BC = first(ax) - 0.5
+# ubound(ax, ::BSpline{Cubic{BC}}, ::OnGrid) where BC = last(ax) + 0.5
 
 """
 `Cubic`: continuity in function value, first and second derivatives yields

src/b-splines/indexing.jl

@@ -4,10 +4,12 @@
     @boundscheck (checkbounds(Bool, itp, x...) || Base.throw_boundserror(itp, x))
     expand_value(itp, x)
 end
-# @inline function (itp::BSplineInterpolation{T,1})(x::Integer, y::Integer) where T
-#     @boundscheck (y == 1 || Base.throw_boundserror(itp, (x,y)))
-#     expand_value(itp, (x,))
-# end
+@propagate_inbounds function (itp::BSplineInterpolation{T,N})(x::Vararg{Number,M}) where {T,M,N}
+    inds, trailing = split_trailing(itp, x)
+    @boundscheck (check1(trailing) || Base.throw_boundserror(itp, x))
+    @assert length(inds) == N
+    itp(inds...)
+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)
@@ -46,7 +48,7 @@ Interpolate `itp` at `x`.
 function expand_value(itp::AbstractInterpolation, x::Tuple)
     coefs = coefficients(itp)
     degree = interpdegree(itp)
-    ixs, rxs = expand_indices_resid(degree, bounds(itp), x)
+    ixs, rxs = splitgrouped(expand_indices_resid(degree, axes(itp), x))
     cxs = expand_weights(value_weights, degree, rxs)
     expand(coefs, cxs, ixs)
 end
@@ -59,7 +61,7 @@ Calculate the interpolated gradient of `itp` at `x`.
 function expand_gradient(itp::AbstractInterpolation, x::Tuple)
     coefs = coefficients(itp)
     degree = interpdegree(itp)
-    ixs, rxs = expand_indices_resid(degree, bounds(itp), x)
+    ixs, rxs = splitgrouped(expand_indices_resid(degree, axes(itp), x))
     cxs = expand_weights(value_weights, degree, rxs)
     gxs = expand_weights(gradient_weights, degree, rxs)
     expand(coefs, (cxs, gxs), ixs)
@@ -68,7 +70,7 @@ end
 function expand_gradient!(dest, itp::AbstractInterpolation, x::Tuple)
     coefs = coefficients(itp)
     degree = interpdegree(itp)
-    ixs, rxs = expand_indices_resid(degree, bounds(itp), x)
+    ixs, rxs = splitgrouped(expand_indices_resid(degree, axes(itp), x))
     cxs = expand_weights(value_weights, degree, rxs)
     gxs = expand_weights(gradient_weights, degree, rxs)
     expand!(dest, coefs, (cxs, gxs), ixs)
@@ -82,7 +84,7 @@ Calculate the interpolated hessian of `itp` at `x`.
 function expand_hessian(itp::AbstractInterpolation, x::Tuple)
     coefs = coefficients(itp)
     degree = interpdegree(itp)
-    ixs, rxs = expand_indices_resid(degree, bounds(itp), x)
+    ixs, rxs = splitgrouped(expand_indices_resid(degree, axes(itp), x))
     cxs = expand_weights(value_weights, degree, rxs)
     gxs = expand_weights(gradient_weights, degree, rxs)
     hxs = expand_weights(hessian_weights, degree, rxs)
@@ -271,17 +273,16 @@ function expand!(dest, coefs, (vweights, gweights, hweights)::NTuple{3,HasNoInte
     dest
 end
 
-function expand_indices_resid(degree, bounds, x)
-    ixbase, δxs = splitpaired(_base_rem(degree, bounds, x))
-    expand_indices(degree, ixbase, bounds, δxs), δxs
+function expand_indices_resid(degree, axs, x)
+    item = expand_index_resid(getfirst(degree), axs[1], x[1])
+    (item, expand_indices_resid(getrest(degree), Base.tail(axs), Base.tail(x))...)
 end
+expand_indices_resid(degree, ::Tuple{}, ::Tuple{}) = ()
 
-@inline _base_rem(degree::Union{Degree,NoInterp}, bounds, x) =
-    (base_rem(degree, bounds[1], x[1]), _base_rem(degree, Base.tail(bounds), Base.tail(x))...)
-@inline _base_rem(degree::Union{Degree,NoInterp}, ::Tuple{}, ::Tuple{}) = ()
-@inline _base_rem(degree::Tuple{Vararg{Union{Degree,NoInterp},N}}, bounds::Tuple{Vararg{Any,N}}, x::Tuple{Vararg{Number,N}}) where N =
-    (base_rem(degree[1], bounds[1], x[1]), _base_rem(Base.tail(degree), Base.tail(bounds), Base.tail(x))...)
-@inline _base_rem(::Tuple{}, ::Tuple{}, ::Tuple{}) = ()
+function expand_index_resid(degree, ax, x::Number)
+    ix, δx = base_rem(degree, ax, x)
+    expand_index(degree, ix, ax, δx), δx
+end
 
 expand_weights(f, degree::Union{Degree,NoInterp}, ixs) =
     (f(degree, ixs[1]), expand_weights(f, degree, Base.tail(ixs))...)
@@ -290,14 +291,14 @@ expand_weights(f, degree::Union{Degree,NoInterp}, ::Tuple{}) = ()
 expand_weights(f, degree::Tuple{Vararg{Union{Degree,NoInterp},N}}, ixs::NTuple{N,Number}) where N =
     f.(degree, ixs)
 
-expand_indices(degree::Union{Degree,NoInterp}, ixs, axs, δxs) =
-    (expand_index(degree, ixs[1], axs[1], δxs[1]), expand_indices(degree, Base.tail(ixs), Base.tail(axs), Base.tail(δxs))...)
-expand_indices(degree::Union{Degree,NoInterp}, ::Tuple{}, ::Tuple{}, ::Tuple{}) = ()
+# expand_indices(degree::Union{Degree,NoInterp}, ixs, axs, δxs) =
+#     (expand_index(degree, ixs[1], axs[1], δxs[1]), expand_indices(degree, Base.tail(ixs), Base.tail(axs), Base.tail(δxs))...)
+# expand_indices(degree::Union{Degree,NoInterp}, ::Tuple{}, ::Tuple{}, ::Tuple{}) = ()
 
-expand_indices(degree::Tuple{Vararg{Union{Degree,NoInterp},N}}, ixs::NTuple{N,Number}, axs::NTuple{N,Tuple{Real,Real}}, δxs::NTuple{N,Number}) where N =
-    expand_index.(degree, ixs, axs, δxs)
+# expand_indices(degree::Tuple{Vararg{Union{Degree,NoInterp},N}}, ixs::NTuple{N,Number}, axs::NTuple{N,Tuple{Real,Real}}, δxs::NTuple{N,Number}) where N =
+#     expand_index.(degree, ixs, axs, δxs)
 
-expand_index(degree, ixs, bounds::Tuple{Real,Real}, δxs) = expand_index(degree, ixs, axfrombounds(bounds), δxs)
+# expand_index(degree, ixs, bounds::Tuple{Real,Real}, δxs) = expand_index(degree, ixs, axfrombounds(bounds), δxs)
 
 checklubounds(ls, us, xs) = _checklubounds(true, ls, us, xs)
 _checklubounds(tf::Bool, ls, us, xs) = _checklubounds(tf & (ls[1] <= xs[1] <= us[1]),
@@ -318,14 +319,19 @@ function getindex_return_type(::Type{BSplineInterpolation{T,N,TCoefs,IT,GT,Pad}}
 end
 
 # This handles round-towards-the-middle for points on half-integer edges
-roundbounds(x, bounds::Tuple{Integer,Integer}) = round(x)
-roundbounds(x, (l, u)) = ifelse(x == l, ceil(l), ifelse(x == u, floor(u), round(x)))
+roundbounds(x::Integer, bounds) = x
+function roundbounds(x, bounds)
+    l, u = first(bounds), last(bounds)
+    h = half(x)
+    xh = x+h
+    ifelse(x < u+half(u), floor(xh), ceil(xh)-1)
+end
 
-floorbounds(x, bounds::Tuple{Integer,Integer}) = floor(x)
-function floorbounds(x, (l, u)::Tuple{Real,Real})
-    ceill = ceil(l)
-    ifelse(l <= x <= ceill, ceill, floor(x))
+floorbounds(x::Integer, ax) = x
+function floorbounds(x, ax)
+    l = first(ax)
+    h = half(x)
+    ifelse(x < l, floor(x+h), floor(x+zero(h)))
 end
 
-axfrombounds((l, u)::Tuple{Integer,Integer}) = UnitRange(l, u)
-axfrombounds((l, u)) = UnitRange(ceil(Int, l), floor(Int, u))
+half(x) = oneunit(x)/2

src/b-splines/linear.jl

@@ -23,12 +23,12 @@ Linear
 
 function base_rem(::Linear, bounds, x)
     xf = floorbounds(x, bounds)
-    xf -= ifelse(xf >= floor(bounds[2]), oneunit(xf), zero(xf))
+    xf -= ifelse(xf >= last(bounds), oneunit(xf), zero(xf))
     δx = x - xf
     fast_trunc(Int, xf), δx
 end
 
-expand_index(::Linear, xi::Number, ax::AbstractUnitRange, δx) = (xi, xi+(δx>0))
+expand_index(::Linear, xi::Number, ax::AbstractUnitRange, δx) = (xi, xi + ((δx!=0) | (xi < last(ax))))
 
 value_weights(::Linear, δx) = (1-δx, δx)
 gradient_weights(::Linear, δx) = (-oneunit(δx), oneunit(δx))

src/b-splines/quadratic.jl

@@ -50,9 +50,9 @@ expand_index(::Quadratic{BC}, xi::Number, ax::AbstractUnitRange, δx) where BC<:
 padded_axis(ax::AbstractUnitRange, ::BSpline{<:Quadratic}) = first(ax)-1:last(ax)+1
 padded_axis(ax::AbstractUnitRange, ::BSpline{Quadratic{BC}}) where BC<:Union{Periodic,InPlace,InPlaceQ} = ax
 
-# Due to padding we can extend the bounds
-lbound(ax, ::BSpline{Quadratic{BC}}, ::OnGrid) where BC = first(ax) - 0.5
-ubound(ax, ::BSpline{Quadratic{BC}}, ::OnGrid) where BC = last(ax) + 0.5
+# # Due to padding we can extend the bounds
+# lbound(ax, ::BSpline{Quadratic{BC}}, ::OnGrid) where BC = first(ax) - 0.5
+# ubound(ax, ::BSpline{Quadratic{BC}}, ::OnGrid) where BC = last(ax) + 0.5
 
 function inner_system_diags(::Type{T}, n::Int, ::Quadratic) where {T}
     du = fill(convert(T, SimpleRatio(1,8)), n-1)