Put the GridType into the BoundaryCondition (fixes #228)

Author: timholy

src/Interpolations.jl

@@ -44,37 +44,52 @@ struct OnCell <: GridType end
 
 const DimSpec{T} = Union{T,Tuple{Vararg{Union{T,NoInterp}}},NoInterp}
 
-abstract type AbstractInterpolation{T,N,IT<:DimSpec{InterpolationType},GT<:DimSpec{GridType}} <: AbstractArray{T,N} end
-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
+abstract type AbstractInterpolation{T,N,IT<:DimSpec{InterpolationType}} <: AbstractArray{T,N} end
+abstract type AbstractInterpolationWrapper{T,N,ITPT,IT} <: AbstractInterpolation{T,N,IT} end
+abstract type AbstractExtrapolation{T,N,ITPT,IT} <: AbstractInterpolationWrapper{T,N,ITPT,IT} 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()`
+- `Throw(gt)`
+- `Flat(gt)`
+- `Line(gt)`
+- `Free(gt)`
+- `Periodic(gt)`
+- `Reflect(gt)`
+- `InPlace(gt)`
+- `InPlaceQ(gt)`
+
+where `gt` is the grid type, e.g., `OnGrid()` or `OnCell()`. `OnGrid` means that the boundary
+condition "activates" at the first and/or last integer location within the interpolation region,
+`OnCell` means the interpolation extends a half-integer beyond the edge before
+activating the boundary condition.
 """
 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
+# Put the gridtype into the boundary condition, since that's all it affects (see issue #228)
+# Nothing is used for extrapolation
+struct Throw{GT<:Union{GridType,Nothing}} <: BoundaryCondition gt::GT end
+struct Flat{GT<:Union{GridType,Nothing}} <: BoundaryCondition gt::GT end
+struct Line{GT<:Union{GridType,Nothing}} <: BoundaryCondition gt::GT end
+struct Free{GT<:Union{GridType,Nothing}} <: BoundaryCondition gt::GT end
+struct Periodic{GT<:Union{GridType,Nothing}} <: BoundaryCondition gt::GT end
+struct Reflect{GT<:Union{GridType,Nothing}} <: BoundaryCondition gt::GT end
+struct InPlace{GT<:Union{GridType,Nothing}} <: BoundaryCondition gt::GT end
 # InPlaceQ is exact for an underlying quadratic. This is nice for ground-truth testing
 # of in-place (unpadded) interpolation.
-struct InPlaceQ <: BoundaryCondition end
+struct InPlaceQ{GT<:Union{GridType,Nothing}} <: BoundaryCondition gt::GT end
 const Natural = Line
 
+(::Type{BC})() where BC<:BoundaryCondition = BC(nothing)
+function Base.show(io::IO, bc::BoundaryCondition)
+    print(io, nameof(typeof(bc)), '(')
+    bc.gt === nothing || show(io, bc.gt)
+    print(io, ')')
+end
+
+
 Base.IndexStyle(::Type{<:AbstractInterpolation}) = IndexCartesian()
 
 size(exp::AbstractExtrapolation) = size(exp.itp)
@@ -85,16 +100,16 @@ twotuple(x, y) = (x, y)
 bounds(itp::AbstractInterpolation) = map(twotuple, lbounds(itp), ubounds(itp))
 bounds(itp::AbstractInterpolation, d) = bounds(itp)[d]
 
-itptype(::Type{AbstractInterpolation{T,N,IT,GT}}) where {T,N,IT<:DimSpec{InterpolationType},GT<:DimSpec{GridType}} = IT
+itptype(::Type{AbstractInterpolation{T,N,IT}}) where {T,N,IT<:DimSpec{InterpolationType}} = IT
 itptype(::Type{ITP}) where {ITP<:AbstractInterpolation} = itptype(supertype(ITP))
 itptype(itp::AbstractInterpolation) = itptype(typeof(itp))
-gridtype(::Type{AbstractInterpolation{T,N,IT,GT}}) where {T,N,IT<:DimSpec{InterpolationType},GT<:DimSpec{GridType}} = GT
+gridtype(::Type{AbstractInterpolation{T,N,IT}}) where {T,N,IT<:DimSpec{InterpolationType}} = GT
 gridtype(::Type{ITP}) where {ITP<:AbstractInterpolation} = gridtype(supertype(ITP))
 gridtype(itp::AbstractInterpolation) = gridtype(typeof(itp))
-ndims(::Type{AbstractInterpolation{T,N,IT,GT}}) where {T,N,IT<:DimSpec{InterpolationType},GT<:DimSpec{GridType}} = N
+ndims(::Type{AbstractInterpolation{T,N,IT}}) where {T,N,IT<:DimSpec{InterpolationType}} = 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{AbstractInterpolation{T,N,IT}}) where {T,N,IT<:DimSpec{InterpolationType}} = T
 eltype(::Type{ITP}) where {ITP<:AbstractInterpolation} = eltype(supertype(ITP))
 eltype(itp::AbstractInterpolation) = eltype(typeof(itp))
 
@@ -209,9 +224,6 @@ import Base: getindex
     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")
@@ -220,5 +232,6 @@ include("scaling/scaling.jl")
 include("utils.jl")
 include("io.jl")
 include("convenience-constructors.jl")
+include("deprecations.jl")
 
 end # module

src/b-splines/b-splines.jl

@@ -8,61 +8,89 @@ export
     Cubic
 
 abstract type Degree{N} <: Flag end
+abstract type DegreeBC{N} <: Degree{N} end  # degree type supporting a BoundaryCondition
 
-struct BSpline{D<:Degree} <: InterpolationType end
-BSpline(::D) where {D<:Degree} = BSpline{D}()
+struct BSpline{D<:Degree} <: InterpolationType
+    degree::D
+end
 
 bsplinetype(::Type{BSpline{D}}) where {D<:Degree} = D
 bsplinetype(::BS) where {BS<:BSpline} = bsplinetype(BS)
 
-degree(::BSpline{D}) where D<:Degree = D()
+degree(mode::BSpline) = mode.degree
 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}
+iscomplete(mode::BSpline) = iscomplete(degree(mode))
+iscomplete(deg::DegreeBC) = _iscomplete(deg.bc.gt)
+iscomplete(deg::Degree) = true
+_iscomplete(::Nothing) = false
+_iscomplete(::GridType) = true
+
+function Base.show(io::IO, bs::BSpline)
+    print(io, "BSpline(")
+    show(io, degree(bs))
+    print(io, ')')
+end
+
+function Base.show(io::IO, deg::DegreeBC)
+    print(io, nameof(typeof(deg)), '(')
+    show(io, deg.bc)
+    print(io, ')')
+end
+
+struct BSplineInterpolation{T,N,TCoefs<:AbstractArray,IT<:DimSpec{BSpline},Axs<:Tuple{Vararg{AbstractUnitRange,N}}} <: AbstractInterpolation{T,N,IT}
     coefs::TCoefs
     parentaxes::Axs
     it::IT
-    gt::GT
 end
-function BSplineInterpolation(::Type{TWeights}, A::AbstractArray{Tel,N}, it::IT, gt::GT, axs) where {N,Tel,TWeights<:Real,IT<:DimSpec{BSpline},GT<:DimSpec{GridType}}
+function BSplineInterpolation(::Type{TWeights}, A::AbstractArray{Tel,N}, it::IT, axs) where {N,Tel,TWeights<:Real,IT<:DimSpec{BSpline}}
     # String interpolation causes allocation, noinline avoids that unless they get called
     @noinline err_concrete(IT) = error("The b-spline type must be a concrete type (was $IT)")
     @noinline warn_concrete(A) = @warn("For performance reasons, consider using an array of a concrete type (typeof(A) == $(typeof(A)))")
+    @noinline err_incomplete(it) = error("OnGrid/OnCell is not supplied for some of the interpolation modes in $it")
 
     isconcretetype(IT) || err_concrete(IT)
     isconcretetype(typeof(A)) || warn_concrete(A)
+    iscomplete(it) || err_incomplete(it)
 
     # Compute the output element type when positions have type TWeights
     if isempty(A)
         T = Base.promote_op(*, TWeights, eltype(A))
     else
         T = typeof(zero(TWeights) * first(A))
     end
-    BSplineInterpolation{T,N,typeof(A),IT,GT,typeof(axs)}(A, fix_axis.(axs), it, gt)
+    BSplineInterpolation{T,N,typeof(A),IT,typeof(axs)}(A, fix_axis.(axs), it)
 end
 
+iscomplete(its::Tuple) = all(iscomplete, its)
+
 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)
 itpflag(itp::BSplineInterpolation) = itp.it
-gridflag(itp::BSplineInterpolation) = itp.gt
 
 size(itp::BSplineInterpolation) = map(length, itp.parentaxes)
 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, 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::AbstractUnitRange, ::BSpline, ::OnCell) = first(ax) - 0.5
-ubound(ax::AbstractUnitRange, ::BSpline, ::OnCell) = last(ax) + 0.5
-lbound(ax::AbstractUnitRange, ::BSpline, ::OnGrid) = first(ax)
-ubound(ax::AbstractUnitRange, ::BSpline, ::OnGrid) = last(ax)
+lbounds(itp::BSplineInterpolation) = _lbounds(itp.parentaxes, itpflag(itp))
+ubounds(itp::BSplineInterpolation) = _ubounds(itp.parentaxes, itpflag(itp))
+_lbounds(axs, itp) = (lbound(axs[1], getfirst(itp)), _lbounds(Base.tail(axs), getrest(itp))...)
+_ubounds(axs, itp) = (ubound(axs[1], getfirst(itp)), _ubounds(Base.tail(axs), getrest(itp))...)
+_lbounds(::Tuple{}, itp) = ()
+_ubounds(::Tuple{}, itp) = ()
+
+lbound(ax::AbstractUnitRange, bs::BSpline)   = lbound(ax, degree(bs))
+lbound(ax::AbstractUnitRange, deg::Degree)   = first(ax)
+lbound(ax::AbstractUnitRange, deg::DegreeBC) = lbound(ax, deg, deg.bc.gt)
+ubound(ax::AbstractUnitRange, bs::BSpline)   = ubound(ax, degree(bs))
+ubound(ax::AbstractUnitRange, deg::Degree)   = last(ax)
+ubound(ax::AbstractUnitRange, deg::DegreeBC) = ubound(ax, deg, deg.bc.gt)
+
+lbound(ax::AbstractUnitRange, ::DegreeBC, ::OnCell) = first(ax) - 0.5
+ubound(ax::AbstractUnitRange, ::DegreeBC, ::OnCell) = last(ax) + 0.5
+lbound(ax::AbstractUnitRange, ::DegreeBC, ::OnGrid) = first(ax)
+ubound(ax::AbstractUnitRange, ::DegreeBC, ::OnGrid) = last(ax)
 
 fix_axis(r::Base.OneTo) = r
 fix_axis(r::Base.Slice) = r
@@ -71,9 +99,9 @@ fix_axis(r::AbstractUnitRange) = fix_axis(UnitRange(r))
 
 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))
+function interpolate(::Type{TWeights}, ::Type{TC}, A, it::IT) where {TWeights,TC,IT<:DimSpec{BSpline}}
+    Apad = prefilter(TWeights, TC, A, it)
+    BSplineInterpolation(TWeights, Apad, it, axes(A))
 end
 
 """
@@ -91,8 +119,8 @@ It may also be a tuple of such values, if you want to use different interpolatio
 
 `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)
+function interpolate(A::AbstractArray, it::IT) where {IT<:DimSpec{BSpline}}
+    interpolate(tweight(A), tcoef(A), A, it)
 end
 
 # We can't just return a tuple-of-types due to julia #12500
@@ -106,14 +134,14 @@ tcoef(A::AbstractArray{Float32}) = Float32
 tcoef(A::AbstractArray{Rational{Int}}) = Rational{Int}
 tcoef(A::AbstractArray{T}) where {T<:Integer} = typeof(float(zero(T)))
 
-function interpolate!(::Type{TWeights}, A, it::IT, gt::GT) where {TWeights,IT<:DimSpec{BSpline},GT<:DimSpec{GridType}}
+function interpolate!(::Type{TWeights}, A::AbstractArray, it::IT) where {TWeights,IT<:DimSpec{BSpline}}
     # Set the bounds of the interpolant inward, if necessary
     axsA = axes(A)
     axspad = padded_axes(axsA, it)
-    BSplineInterpolation(TWeights, prefilter!(TWeights, A, it, gt), it, gt, fix_axis.(padinset.(axsA, axspad)))
+    BSplineInterpolation(TWeights, prefilter!(TWeights, A, it), it, fix_axis.(padinset.(axsA, axspad)))
 end
-function interpolate!(A::AbstractArray, it::IT, gt::GT) where {IT<:DimSpec{BSpline},GT<:DimSpec{GridType}}
-    interpolate!(tweight(A), A, it, gt)
+function interpolate!(A::AbstractArray, it::IT) where {IT<:DimSpec{BSpline}}
+    interpolate!(tweight(A), A, it)
 end
 
 lut!(dl, d, du) = lu!(Tridiagonal(dl, d, du), Val(false))

src/b-splines/cubic.jl

@@ -1,5 +1,8 @@
-struct Cubic{BC<:Flag} <: Degree{3} end
-Cubic(::BC) where {BC<:Flag} = Cubic{BC}()
+struct Cubic{BC<:BoundaryCondition} <: DegreeBC{3}
+    bc::BC
+end
+
+(deg::Cubic)(gt::GridType) = Cubic(deg.bc(gt))
 
 """
 Assuming uniform knots with spacing 1, the `i`th piece of cubic spline
@@ -32,7 +35,7 @@ function base_rem(::Cubic, bounds, x)
 end
 
 expand_index(::Cubic{BC}, xi::Number, ax::AbstractUnitRange, δx) where BC = (xi-1, xi, xi+1, xi+2)
-expand_index(::Cubic{Periodic}, xi::Number, ax::AbstractUnitRange, δx) =
+expand_index(::Cubic{Periodic{GT}}, xi::Number, ax::AbstractUnitRange, δx) where GT<:GridType =
     (modrange(xi-1, ax), modrange(xi, ax), modrange(xi+1, ax), modrange(xi+2, ax))
 
 # expand_coefs(::Type{BSpline{Cubic{BC}}}, δx) = cvcoefs(δx)
@@ -45,31 +48,31 @@ expand_index(::Cubic{Periodic}, xi::Number, ax::AbstractUnitRange, δx) =
 #     end
 # end
 
-function value_weights(::Cubic, δx)
+function value_weights(::BSpline{<:Cubic}, δx)
     x3, xcomp3 = cub(δx), cub(1-δx)
     (SimpleRatio(1,6) * xcomp3,
      SimpleRatio(2,3) - sqr(δx) + SimpleRatio(1,2)*x3,
      SimpleRatio(2,3) - sqr(1-δx) + SimpleRatio(1,2)*xcomp3,
      SimpleRatio(1,6) * x3)
 end
 
-function gradient_weights(::Cubic, δx)
+function gradient_weights(::BSpline{<:Cubic}, δx)
     x2, xcomp2 = sqr(δx), sqr(1-δx)
     (-SimpleRatio(1,2) * xcomp2,
      -2*δx + SimpleRatio(3,2)*x2,
      +2*(1-δx) - SimpleRatio(3,2)*xcomp2,
      SimpleRatio(1,2) * x2)
 end
 
-hessian_weights(::Cubic, δx) = (1-δx, 3*δx-2, 3*(1-δx)-2, δx)
+hessian_weights(::BSpline{<:Cubic}, δx) = (1-δx, 3*δx-2, 3*(1-δx)-2, δx)
 
 
 # ------------ #
 # Prefiltering #
 # ------------ #
 
 padded_axis(ax::AbstractUnitRange, ::BSpline{<:Cubic}) = first(ax)-1:last(ax)+1
-padded_axis(ax::AbstractUnitRange, ::BSpline{Cubic{Periodic}}) = ax
+padded_axis(ax::AbstractUnitRange, ::BSpline{Cubic{Periodic{GT}}}) where GT<:GridType = ax
 
 # # Due to padding we can extend the bounds
 # lbound(ax, ::BSpline{Cubic{BC}}, ::OnGrid) where BC = first(ax) - 0.5
@@ -97,7 +100,7 @@ Applying this condition yields
     -cm + cp = 0
 """
 function prefiltering_system(::Type{T}, ::Type{TC}, n::Int,
-                             degree::Cubic{Flat}, ::OnGrid) where {T,TC}
+                             degree::Cubic{Flat{OnGrid}}) where {T,TC}
     dl, d, du = inner_system_diags(T, n, degree)
     d[1] = d[end] = -oneunit(T)
     du[1] = dl[end] = zero(T)
@@ -123,7 +126,7 @@ were to use `y_0'(x)` we would have to introduce new coefficients, so that would
 close the system. Instead, we extend the outermost polynomial for an extra half-cell.)
 """
 function prefiltering_system(::Type{T}, ::Type{TC}, n::Int,
-                             degree::Cubic{Flat}, ::OnCell) where {T,TC}
+                             degree::Cubic{Flat{OnCell}}) where {T,TC}
     dl, d, du = inner_system_diags(T,n,degree)
     d[1] = d[end] = -9
     du[1] = dl[end] = 11
@@ -152,7 +155,7 @@ were to use `y_0'(x)` we would have to introduce new coefficients, so that would
 close the system. Instead, we extend the outermost polynomial for an extra half-cell.)
 """
 function prefiltering_system(::Type{T}, ::Type{TC}, n::Int,
-                             degree::Cubic{Line}, ::OnCell) where {T,TC}
+                             degree::Cubic{Line{OnCell}}) where {T,TC}
     dl,d,du = inner_system_diags(T,n,degree)
     d[1] = d[end] = 3
     du[1] = dl[end] = -7
@@ -177,7 +180,7 @@ condition gives:
     1 cm -2 c + 1 cp = 0
 """
 function prefiltering_system(::Type{T}, ::Type{TC}, n::Int,
-                             degree::Cubic{Line}, ::OnGrid) where {T,TC}
+                             degree::Cubic{Line{OnGrid}}) where {T,TC}
     dl,d,du = inner_system_diags(T,n,degree)
     d[1] = d[end] = 1
     du[1] = dl[end] = -2
@@ -201,7 +204,7 @@ as periodic, yielding
 where `N` is the number of data points.
 """
 function prefiltering_system(::Type{T}, ::Type{TC}, n::Int,
-                             degree::Cubic{Periodic}, ::GridType) where {T,TC}
+                             degree::Cubic{<:Periodic}) where {T,TC}
     dl, d, du = inner_system_diags(T,n,degree)
 
     specs = WoodburyMatrices.sparse_factors(T, n,
@@ -220,7 +223,7 @@ continuous derivative at the second-to-last cell boundary; this means
     1 cm -3 c + 3 cp -1 cpp = 0
 """
 function prefiltering_system(::Type{T}, ::Type{TC}, n::Int,
-                             degree::Cubic{Free}, ::GridType) where {T,TC}
+                             degree::Cubic{<:Free}) where {T,TC}
     dl, d, du = inner_system_diags(T,n,degree)
 
     specs = WoodburyMatrices.sparse_factors(T, n,

src/b-splines/linear.jl

@@ -1,4 +1,4 @@
-struct Linear <: Degree{1} end
+struct Linear <: Degree{1} end  # boundary conditions not supported
 
 """
 Assuming uniform knots with spacing 1, the `i`th peice of linear b-spline