Implement the new scheme for BSpline (and partial for NoInterp)

This makes several major changes:
- it eliminates the generated functions and switches emphasis to the value domain rather than the type domain
- it uses OffsetArrays to handle the padding. When using interpolate!, the
  axes may be narrowed to account for the boundary conditions.
- now, the axes of the output correspond to the region where we can guarantee
  that we reconstruct the original array values
- switches to the bounds-checking framework in base (may not be working yet)

Author: timholy

REQUIRE

@@ -4,3 +4,4 @@ WoodburyMatrices 0.1.5
 Ratios
 AxisAlgorithms 0.3.0
 OffsetArrays
+StaticArrays

src/Interpolations.jl

@@ -30,10 +30,10 @@ export
     # scaling/scaling.jl
 
 using LinearAlgebra, SparseArrays
-using WoodburyMatrices, Ratios, AxisAlgorithms, OffsetArrays
+using StaticArrays, WoodburyMatrices, Ratios, AxisAlgorithms, OffsetArrays
 
-import Base: convert, size, axes, getindex, promote_rule,
-             ndims, eltype, checkbounds
+using Base: @propagate_inbounds
+import Base: convert, size, axes, promote_rule, ndims, eltype, checkbounds
 
 abstract type Flag end
 abstract type InterpolationType <: Flag end
@@ -48,57 +48,183 @@ abstract type AbstractInterpolation{T,N,IT<:DimSpec{InterpolationType},GT<:DimSp
 abstract type AbstractInterpolationWrapper{T,N,ITPT,IT,GT} <: AbstractInterpolation{T,N,IT,GT} end
 abstract type AbstractExtrapolation{T,N,ITPT,IT,GT} <: AbstractInterpolationWrapper{T,N,ITPT,IT,GT} end
 
-struct Throw <: Flag end
-struct Flat <: Flag end
-struct Line <: Flag end
-struct Free <: Flag end
-struct Periodic <: Flag end
-struct Reflect <: Flag end
-struct InPlace <: Flag end
+"""
+    BoundaryCondition
+
+An abstract type with one of the following values (see the help for each for details):
+
+- `Throw()`
+- `Flat()`
+- `Line()`
+- `Free()`
+- `Periodic()`
+- `Reflect()`
+- `InPlace()`
+- `InPlaceQ()`
+"""
+abstract type BoundaryCondition <: Flag end
+struct Throw <: BoundaryCondition end
+struct Flat <: BoundaryCondition end
+struct Line <: BoundaryCondition end
+struct Free <: BoundaryCondition end
+struct Periodic <: BoundaryCondition end
+struct Reflect <: BoundaryCondition end
+struct InPlace <: BoundaryCondition end
 # InPlaceQ is exact for an underlying quadratic. This is nice for ground-truth testing
 # of in-place (unpadded) interpolation.
-struct InPlaceQ <: Flag end
+struct InPlaceQ <: BoundaryCondition end
 const Natural = Line
 
-@generated size(itp::AbstractInterpolation{T,N}) where {T, N} = Expr(:tuple, [:(size(itp, $i)) for i in 1:N]...)
-size(exp::AbstractExtrapolation, d) = size(exp.itp, d)
-bounds(itp::AbstractInterpolation{T,N}) where {T,N} = tuple(zip(lbounds(itp), ubounds(itp))...)
-bounds(itp::AbstractInterpolation{T,N}, d) where {T,N} = (lbound(itp,d),ubound(itp,d))
-@generated lbounds(itp::AbstractInterpolation{T,N}) where {T,N} = Expr(:tuple, [:(lbound(itp, $i)) for i in 1:N]...)
-@generated ubounds(itp::AbstractInterpolation{T,N}) where {T,N} = Expr(:tuple, [:(ubound(itp, $i)) for i in 1:N]...)
-lbound(itp::AbstractInterpolation{T,N}, d) where {T,N} = 1
-ubound(itp::AbstractInterpolation{T,N}, d) where {T,N} = size(itp, d)
+Base.IndexStyle(::Type{<:AbstractInterpolation}) = IndexCartesian()
+
+size(itp::AbstractInterpolation) = map(length, itp.parentaxes)
+size(exp::AbstractExtrapolation) = size(exp.itp)
+axes(itp::AbstractInterpolation) = itp.parentaxes
+axes(exp::AbstractExtrapolation) = axes(exp.itp)
+
+bounds(itp::AbstractInterpolation) = map(twotuple, lbounds(itp), ubounds(itp))
+twotuple(r::AbstractUnitRange) = (first(r), last(r))
+twotuple(x, y) = (x, y)
+bounds(itp::AbstractInterpolation, d) = bounds(itp)[d]
+lbounds(itp::AbstractInterpolation) = _lbounds(itp.parentaxes, itpflag(itp), gridflag(itp))
+ubounds(itp::AbstractInterpolation) = _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))
+
 itptype(::Type{AbstractInterpolation{T,N,IT,GT}}) where {T,N,IT<:DimSpec{InterpolationType},GT<:DimSpec{GridType}} = IT
 itptype(::Type{ITP}) where {ITP<:AbstractInterpolation} = itptype(supertype(ITP))
-itptype(itp::AbstractInterpolation ) = itptype(typeof(itp))
+itptype(itp::AbstractInterpolation) = itptype(typeof(itp))
+itpflag(itp::AbstractInterpolation) = itp.it
 gridtype(::Type{AbstractInterpolation{T,N,IT,GT}}) where {T,N,IT<:DimSpec{InterpolationType},GT<:DimSpec{GridType}} = GT
 gridtype(::Type{ITP}) where {ITP<:AbstractInterpolation} = gridtype(supertype(ITP))
 gridtype(itp::AbstractInterpolation) = gridtype(typeof(itp))
+gridflag(itp::AbstractInterpolation) = itp.gt
 ndims(::Type{AbstractInterpolation{T,N,IT,GT}}) where {T,N,IT<:DimSpec{InterpolationType},GT<:DimSpec{GridType}} = N
 ndims(::Type{ITP}) where {ITP<:AbstractInterpolation} = ndims(supertype(ITP))
 ndims(itp::AbstractInterpolation) = ndims(typeof(itp))
 eltype(::Type{AbstractInterpolation{T,N,IT,GT}}) where {T,N,IT<:DimSpec{InterpolationType},GT<:DimSpec{GridType}} = T
 eltype(::Type{ITP}) where {ITP<:AbstractInterpolation} = eltype(supertype(ITP))
 eltype(itp::AbstractInterpolation) = eltype(typeof(itp))
-count_interp_dims(::Type{T}, N) where {T<:AbstractInterpolation} = N
 
-# Generic indexing methods (fixes #183)
-Base.to_index(::AbstractInterpolation, i::Real) = i
+"""
+    n = count_interp_dims(ITP)
+
+Count the number of dimensions along which type `ITP` is interpolating.
+`NoInterp` dimensions do not contribute to the sum.
+"""
+count_interp_dims(::Type{ITP}) where {ITP<:AbstractInterpolation} = count_interp_dims(itptype(ITP), ndims(ITP))
+count_interp_dims(::Type{IT}, n) where {IT<:InterpolationType} = n * count_interp_dims(IT)
+count_interp_dims(it::Type{IT}, n) where IT<:Tuple{Vararg{InterpolationType,N}} where N =
+    _count_interp_dims(0, it...)
+@inline _count_interp_dims(c, ::IT1, args...) where IT1 =
+    _count_interp_dims(c + count_interp_dims(IT1), args...)
+_count_interp_dims(c) = c
+
+"""
+    w = value_weights(degree, δx)
+
+Compute the weights for interpolation of the value at an offset `δx` from the "base" position.
+`degree` describes the interpolation scheme.
+
+# Example
+
+```jldoctest
+julia> Interpolations.value_weights(Linear(), 0.2)
+(0.8, 0.2)
+```
+
+This corresponds to the fact that linear interpolation at `x + 0.2` is `0.8*y[x] + 0.2*y[x+1]`.
+"""
+function value_weights end
+
+"""
+    w = gradient_weights(degree, δx)
+
+Compute the weights for interpolation of the gradient at an offset `δx` from the "base" position.
+`degree` describes the interpolation scheme.
+
+# Example
+
+```jldoctest
+julia> Interpolations.gradient_weights(Linear(), 0.2)
+(-1.0, 1.0)
+```
+
+This defines the gradient of a linear interpolation at 3.2 as `y[4] - y[3]`.
+"""
+function gradient_weights end
+
+"""
+    w = hessian_weights(degree, δx)
+
+Compute the weights for interpolation of the hessian at an offset `δx` from the "base" position.
+`degree` describes the interpolation scheme.
 
-@inline function Base._getindex(::IndexCartesian, A::AbstractInterpolation{T,N}, I::Vararg{Int,N}) where {T,N} # ambiguity resolution
-    @inbounds r = getindex(A, I...)
-    r
+# Example
+
+```jldoctest
+julia> Interpolations.hessian_weights(Linear(), 0.2)
+(0.0, 0.0)
+```
+
+Linear interpolation uses straight line segments, so the second derivative is zero.
+"""
+function hessian_weights end
+
+
+gradient1(itp::AbstractInterpolation{T,1}, x) where {T} = gradient(itp, x)[1]
+hessian1(itp::AbstractInterpolation{T,1}, x) where {T} = hessian(itp, x)[1]
+
+### Supporting expansion of CartesianIndex
+
+const UnexpandedIndexTypes = Union{Number, AbstractVector, CartesianIndex}
+const ExpandedIndexTypes = Union{Number, AbstractVector}
+
+Base.to_index(::AbstractInterpolation, x::Number) = x
+
+# Commented out because you can't add methods to an abstract type.
+# @inline function (itp::AbstractInterpolation)(x::Vararg{UnexpandedIndexTypes})
+#     itp(to_indices(itp, x))
+# end
+function gradient(itp::AbstractInterpolation, x::Vararg{UnexpandedIndexTypes})
+    gradient(itp, to_indices(itp, x))
 end
-@inline function Base._getindex(::IndexCartesian, A::AbstractInterpolation{T,N}, I::Vararg{Real,N}) where {T,N}
-    @inbounds r = getindex(A, I...)
-    r
+function gradient!(dest, itp::AbstractInterpolation, x::Vararg{UnexpandedIndexTypes})
+    gradient!(dest, itp, to_indices(itp, x))
 end
+function hessian(itp::AbstractInterpolation, x::Vararg{UnexpandedIndexTypes})
+    hessian(itp, to_indices(itp, x))
+end
+function hessian!(dest, itp::AbstractInterpolation, x::Vararg{UnexpandedIndexTypes})
+    hessian!(dest, itp, to_indices(itp, x))
+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
+
+# 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...)
+@propagate_inbounds function getindex(itp::AbstractInterpolation{T,1}, i::Integer, j::Integer) where T
+    @boundscheck (j == 1 || Base.throw_boundserror(itp, (i, j)))
+    itp(i)
+end
+
+# deprecate getindex for other numeric indices
+@deprecate getindex(itp::AbstractInterpolation{T,N}, i::Vararg{Number,N}) where {T,N} itp(i...)
 
 include("nointerp/nointerp.jl")
 include("b-splines/b-splines.jl")
-include("gridded/gridded.jl")
-include("extrapolation/extrapolation.jl")
-include("scaling/scaling.jl")
+# include("gridded/gridded.jl")
+# include("extrapolation/extrapolation.jl")
+# include("scaling/scaling.jl")
 include("utils.jl")
 include("io.jl")
 include("convenience-constructors.jl")

src/b-splines/b-splines.jl

@@ -13,6 +13,10 @@ struct BSpline{D<:Degree} <: InterpolationType end
 BSpline(::D) where {D<:Degree} = BSpline{D}()
 
 bsplinetype(::Type{BSpline{D}}) where {D<:Degree} = D
+bsplinetype(::BS) where {BS<:BSpline} = bsplinetype(BS)
+
+degree(::BSpline{D}) where D<:Degree = D()
+degree(::NoInterp) = NoInterp()
 
 struct BSplineInterpolation{T,N,TCoefs<:AbstractArray,IT<:DimSpec{BSpline},GT<:DimSpec{GridType},Axs<:Tuple{Vararg{AbstractUnitRange,N}}} <: AbstractInterpolation{T,N,IT,GT}
     coefs::TCoefs
@@ -37,53 +41,44 @@ function BSplineInterpolation(::Type{TWeights}, A::AbstractArray{Tel,N}, it::IT,
     BSplineInterpolation{T,N,typeof(A),IT,GT,typeof(axs)}(A, fix_axis.(axs), it, gt)
 end
 
-# Utilities for working either with scalars or tuples/tuple-types
-iextract(::Type{T}, d) where {T<:BSpline} = T
-iextract(t, d) = t.parameters[d]
-padding(::Type{BSplineInterpolation{T,N,TCoefs,IT,GT,pad}}) where {T,N,TCoefs,IT,GT,pad} = pad
-padding(itp::AbstractInterpolation) = padding(typeof(itp))
-padextract(pad::Integer, d) = pad
-padextract(pad::Tuple{Vararg{Integer}}, d) = pad[d]
-
-lbound(itp::BSplineInterpolation{T,N,TCoefs,IT,OnGrid}, d::Integer) where {T,N,TCoefs,IT} =
-    first(axes(itp, d))
-ubound(itp::BSplineInterpolation{T,N,TCoefs,IT,OnGrid}, d::Integer) where {T,N,TCoefs,IT} =
-    last(axes(itp, d))
-lbound(itp::BSplineInterpolation{T,N,TCoefs,IT,OnCell}, d::Integer) where {T,N,TCoefs,IT} =
-    first(axes(itp, d)) - 0.5
-ubound(itp::BSplineInterpolation{T,N,TCoefs,IT,OnCell}, d::Integer) where {T,N,TCoefs,IT} =
-    last(axes(itp, d))+0.5
-
-lbound(itp::BSplineInterpolation{T,N,TCoefs,IT,OnGrid}, d, inds) where {T,N,TCoefs,IT} =
-    first(inds)
-ubound(itp::BSplineInterpolation{T,N,TCoefs,IT,OnGrid}, d, inds) where {T,N,TCoefs,IT} =
-    last(inds)
-lbound(itp::BSplineInterpolation{T,N,TCoefs,IT,OnCell}, d, inds) where {T,N,TCoefs,IT} =
-    first(inds) - 0.5
-ubound(itp::BSplineInterpolation{T,N,TCoefs,IT,OnCell}, d, inds) where {T,N,TCoefs,IT} =
-    last(inds)+0.5
-
-count_interp_dims(::Type{BSplineInterpolation{T,N,TCoefs,IT,GT,pad}}, n) where {T,N,TCoefs,IT<:DimSpec{InterpolationType},GT<:DimSpec{GridType},pad} = count_interp_dims(IT, n)
-
-function size(itp::BSplineInterpolation{T,N,TCoefs,IT,GT,pad}, d) where {T,N,TCoefs,IT,GT,pad}
-    d <= N ? size(itp.coefs, d) - 2*padextract(pad, d) : 1
-end
+coefficients(itp::BSplineInterpolation) = itp.coefs
+interpdegree(itp::BSplineInterpolation) = interpdegree(itpflag(itp))
+interpdegree(::BSpline{T}) where T = T()
+interpdegree(it::Tuple{Vararg{Union{BSpline,NoInterp},N}}) where N = interpdegree.(it)
 
-@inline axes(itp::BSplineInterpolation{T,N,TCoefs,IT,GT,pad}) where {T,N,TCoefs,IT,GT,pad} =
-    indices_removepad.(axes(itp.coefs), pad)
+# The unpadded defaults
+lbound(ax, ::BSpline, ::OnCell) = first(ax) - 0.5
+ubound(ax, ::BSpline, ::OnCell) = last(ax) + 0.5
+lbound(ax, ::BSpline, ::OnGrid) = first(ax)
+ubound(ax, ::BSpline, ::OnGrid) = last(ax)
 
 fix_axis(r::Base.OneTo) = r
 fix_axis(r::Base.Slice) = r
 fix_axis(r::UnitRange) = Base.Slice(r)
 fix_axis(r::AbstractUnitRange) = fix_axis(UnitRange(r))
-function axes(itp::BSplineInterpolation{T,N,TCoefs,IT,GT,pad}, d) where {T,N,TCoefs,IT,GT,pad}
-    d <= N ? indices_removepad(axes(itp.coefs, d), padextract(pad, d)) : axes(itp.coefs, d)
-end
+
+count_interp_dims(::Type{BSI}, n) where BSI<:BSplineInterpolation = count_interp_dims(itptype(BSI), n)
 
 function interpolate(::Type{TWeights}, ::Type{TC}, A, it::IT, gt::GT) where {TWeights,TC,IT<:DimSpec{BSpline},GT<:DimSpec{GridType}}
     Apad = prefilter(TWeights, TC, A, it, gt)
     BSplineInterpolation(TWeights, Apad, it, gt, axes(A))
 end
+
+"""
+    itp = interpolate(A, interpmode, gridstyle)
+
+Interpolate an array `A` in the mode determined by `interpmode` and `gridstyle`.
+`interpmode` may be one of
+
+- `BSpline(NoInterp())`
+- `BSpline(Linear())`
+- `BSpline(Quadratic(BC()))` (see [`BoundaryCondition`](@ref))
+- `BSpline(Cubic(BC()))`
+
+It may also be a tuple of such values, if you want to use different interpolation schemes along each axis.
+
+`gridstyle` should be one of `OnGrid()` or `OnCell()`.
+"""
 function interpolate(A::AbstractArray, it::IT, gt::GT) where {IT<:DimSpec{BSpline},GT<:DimSpec{GridType}}
     interpolate(tweight(A), tcoef(A), A, it, gt)
 end
@@ -109,24 +104,7 @@ function interpolate!(A::AbstractArray, it::IT, gt::GT) where {IT<:DimSpec{BSpli
     interpolate!(tweight(A), A, it, gt)
 end
 
-offsetsym(off, d) = off == -1 ? Symbol("ixm_", d) :
-                    off ==  0 ? Symbol("ix_", d) :
-                    off ==  1 ? Symbol("ixp_", d) :
-                    off ==  2 ? Symbol("ixpp_", d) : error("offset $off not recognized")
-
-# Ideally we might want to shift the indices symmetrically, but this
-# would introduce an inconsistency, so we just append on the right
-@inline indices_removepad(inds::Base.OneTo, pad) = Base.OneTo(length(inds) - 2*pad)
-@inline indices_removepad(inds, pad) = oftype(inds, first(inds):last(inds) - 2*pad)
-@inline indices_addpad(inds::Base.OneTo, pad) = Base.OneTo(length(inds) + 2*pad)
-@inline indices_addpad(inds, pad) = oftype(inds, first(inds):last(inds) + 2*pad)
-@inline indices_interior(inds, pad) = first(inds)+pad:last(inds)-pad
-
-@static if VERSION < v"0.7.0-DEV.3449"
-    lut!(dl, d, du) = lufact!(Tridiagonal(dl, d, du), Val{false})
-else
-    lut!(dl, d, du) = lu!(Tridiagonal(dl, d, du), Val(false))
-end
+lut!(dl, d, du) = lu!(Tridiagonal(dl, d, du), Val(false))
 
 include("constant.jl")
 include("linear.jl")
@@ -137,4 +115,5 @@ include("prefiltering.jl")
 include("../filter1d.jl")
 
 Base.parent(A::BSplineInterpolation{T,N,TCoefs,UT}) where {T,N,TCoefs,UT<:Union{BSpline{Linear},BSpline{Constant}}} = A.coefs
-Base.parent(A::BSplineInterpolation{T,N,TCoefs,UT}) where {T,N,TCoefs,UT} = throw(ArgumentError("The given BSplineInterpolation does not serve as a \"view\" for a parent array. This would only be true for Constant and Linear b-splines."))
+Base.parent(A::BSplineInterpolation{T,N,TCoefs,UT}) where {T,N,TCoefs,UT} =
+    throw(ArgumentError("The given BSplineInterpolation does not serve as a \"view\" for a parent array. This would only be true for Constant and Linear b-splines."))