Merge branch 'master' into master

Author: mkitti

docs/src/NEWS.md

@@ -0,0 +1 @@
+../../NEWS.md

docs/src/index.md

@@ -5,7 +5,7 @@
 [![PkgEval Status](http://pkg.julialang.org/badges/Interpolations_0.4.svg)](http://pkg.julialang.org/?pkg=Interpolations)
 [![Interpolations](http://pkg.julialang.org/badges/Interpolations_0.5.svg)](http://pkg.julialang.org/?pkg=Interpolations)
 
-**NEWS** v0.9 was a breaking release. See the [news](../../NEWS.md) for details on how to update.
+**NEWS** v0.9 was a breaking release. See the [news](NEWS.md) for details on how to update.
 
 This package implements a variety of interpolation schemes for the
 Julia language.  It has the goals of ease-of-use, broad algorithmic

perf/benchmarks.jl

@@ -30,14 +30,14 @@ end
 strip_prefix(str::AbstractString) = replace(str, "Interpolations."=>"")
 benchstr(::Type{T}) where {T<:Interpolations.GridType} = strip_prefix(string(T))
 
-benchstr(::Type{Constant}) = "Constant()"
+benchstr(::Type{<:Constant}) = "Constant()"
 benchstr(::Type{Linear}) = "Linear()"
 benchstr(::Type{Quadratic{BC}}, ::Type{GT}) where {BC<:Interpolations.BoundaryCondition,GT<:Interpolations.GridType} =
     string("Quadratic(", strip_prefix(string(BC)), "(", strip_prefix(string(GT)), "()))")
 benchstr(::Type{Cubic{BC}}, ::Type{GT}) where {BC<:Interpolations.BoundaryCondition,GT<:Interpolations.GridType} =
     string("Cubic(", strip_prefix(string(BC)), "(", strip_prefix(string(GT)), "()))")
 
-groupstr(::Type{Constant}) = "constant"
+groupstr(::Type{<:Constant}) = "constant"
 groupstr(::Type{Linear}) = "linear"
 groupstr(::Type{Quadratic}) = "quadratic"
 groupstr(::Type{Cubic}) = "cubic"

src/b-splines/b-splines.jl

@@ -161,6 +161,11 @@ function interpolate(::Type{TWeights}, ::Type{TC}, A, it::IT) where {TWeights,TC
     BSplineInterpolation(TWeights, Apad, it, axes(A))
 end
 
+function interpolate(::Type{TWeights}, ::Type{TC}, A, it::IT, λ::Real, k::Int) where {TWeights,TC,IT<:DimSpec{BSpline}}
+    Apad = prefilter(TWeights, TC, A, it, λ, k)
+    BSplineInterpolation(TWeights, Apad, it, axes(A))
+end
+
 """
     itp = interpolate(A, interpmode)
 
@@ -179,6 +184,33 @@ function interpolate(A::AbstractArray, it::IT) where {IT<:DimSpec{BSpline}}
     interpolate(tweight(A), tcoef(A), A, it)
 end
 
+"""
+    itp = interpolate(A, interpmode, gridstyle, λ, k)
+
+Interpolate an array `A` in the mode determined by `interpmode` and `gridstyle`
+with regularization following [1], of order `k` and constant `λ`. 
+`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()`.
+
+`k` corresponds to the derivative to penalize. In the limit λ->∞, the interpolation function is
+a polynomial of order `k-1`. A value of 2 is the most common.
+
+`λ` is non-negative. If its value is zero, it falls back to non-regularized interpolation.
+
+[1] https://projecteuclid.org/euclid.ss/1038425655.
+"""
+function interpolate(A::AbstractArray, it::IT, λ::Real, k::Int) where {IT<:DimSpec{BSpline}}
+    interpolate(tweight(A), tcoef(A), A, it, λ, k)
+end
+
 # We can't just return a tuple-of-types due to julia #12500
 tweight(A::AbstractArray) = Float64
 tweight(A::AbstractArray{T}) where T<:AbstractFloat = T
@@ -201,6 +233,16 @@ function interpolate!(A::AbstractArray, it::IT) where {IT<:DimSpec{BSpline}}
     interpolate!(tweight(A), A, it)
 end
 
+function interpolate!(::Type{TWeights}, A::AbstractArray, it::IT, λ::Real, k::Int) 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, λ, k), it, fix_axis.(padinset.(axsA, axspad)))
+end
+function interpolate!(A::AbstractArray, it::IT, λ::Real, k::Int) where {IT<:DimSpec{BSpline}}
+    interpolate!(tweight(A), A, it, λ, k)
+end
+
 lut!(dl, d, du) = lu!(Tridiagonal(dl, d, du), Val(false))
 
 include("constant.jl")
@@ -211,6 +253,6 @@ include("indexing.jl")
 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<: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."))

src/b-splines/constant.jl

@@ -1,12 +1,33 @@
-struct Constant <: Degree{0} end
+export Nearest, Previous, Next
+
+abstract type ConstantInterpType end
+struct Nearest <: ConstantInterpType end
+struct Previous <: ConstantInterpType end
+struct Next <: ConstantInterpType end
+
+struct Constant{T<:ConstantInterpType} <: Degree{0}
+    Constant() = new{Nearest}()
+    Constant{Previous}() = new{Previous}()
+    Constant{Next}() = new{Next}()
+end
 
 """
 Constant b-splines are *nearest-neighbor* interpolations, and effectively
 return `A[round(Int,x)]` when interpolating.
 """
 Constant
 
-function positions(::Constant, ax, x)  # discontinuity occurs at half-integer locations
+function positions(c::Constant{Previous}, ax, x)  # discontinuity occurs at integer locations
+    xm = floorbounds(x, ax)
+    δx = x - xm
+    fast_trunc(Int, xm), δx
+end
+function positions(c::Constant{Next}, ax, x)  # discontinuity occurs at integer locations
+    xm = ceilbounds(x, ax)
+    δx = x - xm
+    fast_trunc(Int, xm), δx
+end
+function positions(c::Constant{Nearest}, ax, x)  # discontinuity occurs at half-integer locations
     xm = roundbounds(x, ax)
     δx = x - xm
     fast_trunc(Int, xm), δx
@@ -16,4 +37,4 @@ value_weights(::Constant, δx) = (1,)
 gradient_weights(::Constant, δx) = (0,)
 hessian_weights(::Constant, δx) = (0,)
 
-padded_axis(ax::AbstractUnitRange, ::BSpline{Constant}) = ax
+padded_axis(ax::AbstractUnitRange, ::BSpline{<:Constant}) = ax

src/b-splines/indexing.jl

@@ -199,6 +199,16 @@ function _floorbounds(x, ax::Union{Tuple{Real,Real}, AbstractUnitRange})
     ifelse(x < l, floor(x+h), floor(x+zero(h)))
 end
 
+ceilbounds(x::Integer, ax::Tuple{Real,Real}) = x
+ceilbounds(x::Integer, ax::AbstractUnitRange) = x
+ceilbounds(x, ax::Tuple{Real,Real}) = _ceilbounds(x, ax)
+ceilbounds(x, ax::AbstractUnitRange) = _ceilbounds(x, ax)
+function _ceilbounds(x, ax::Union{Tuple{Real,Real}, AbstractUnitRange})
+    u = last(ax)
+    h = half(x)
+    ifelse(x > u, ceil(x+h), ceil(x+zero(h)))
+end
+
 half(x) = oneunit(x)/2
 
 symmatrix(h::NTuple{1,Any}) = SMatrix{1,1}(h)

src/b-splines/prefiltering.jl

@@ -40,10 +40,18 @@ function prefilter(
     prefilter!(TWeights, ret, it)
 end
 
+function prefilter(
+    ::Type{TWeights}, ::Type{TC}, A::AbstractArray,
+    it::Union{BSpline,Tuple{Vararg{Union{BSpline,NoInterp}}}},
+    λ::Real, k::Int
+    ) where {TWeights,TC}
+    ret = copy_with_padding(TC, A, it)
+    prefilter!(TWeights, ret, it, λ, k)
+end
+
 function prefilter!(
     ::Type{TWeights}, ret::TCoefs, it::BSpline
     ) where {TWeights,TCoefs<:AbstractArray}
-    local buf, shape, retrs
     sz = size(ret)
     first = true
     for dim in 1:ndims(ret)
@@ -53,12 +61,57 @@ function prefilter!(
     ret
 end
 
+function diffop(::Type{TWeights}, n::Int, k::Int) where {TWeights}
+    D = spdiagm(0 => ones(TWeights,n))
+    for i in 1:k
+        D = diff(D; dims=1)
+    end
+    ## TODO: Normalize by n?
+    D' * D
+end
+diffop(n::Int, k::Int) = diffop(Float64, n, k)
+### TODO: add compiled constructor for most common operators of order k=1,2
+
+
+function prefilter!(
+    ::Type{TWeights}, ret::TCoefs, it::BSpline, λ::Real, k::Int
+    ) where {TWeights,TCoefs<:AbstractArray}
+
+    sz = size(ret)
+
+    # Test if regularizing is allowed
+    fallback = if ndims(ret) > 1
+        @warn "Smooth BSpline only available for Vectors, fallback to non-smoothed"
+        true
+    elseif λ < 0
+        @warn "Smooth BSpline require a non-negative λ, fallback to non-smoothed: $(λ)"
+        true
+    elseif λ == 0
+        true
+    else
+        false
+    end
+    
+    if fallback
+        return prefilter!(TWeights, ret, it)
+    end
+
+    M, b = prefiltering_system(TWeights, eltype(TCoefs), sz[1], degree(it))
+    ## Solve with regularization
+    # Convert to dense Matrix because `*(Woodbury{Adjoint}, Woodbury)` is not defined
+    Mt = Matrix(M)'
+    Q = diffop(TWeights, sz[1], k)
+    K = cholesky(Mt * Matrix(M) + λ * Q)
+    B = Mt * popwrapper(ret)
+    ldiv!(popwrapper(ret), K, B)
+
+    ret
+end
+
 function prefilter!(
     ::Type{TWeights}, ret::TCoefs, its::Tuple{Vararg{Union{BSpline,NoInterp}}}
     ) where {TWeights,TCoefs<:AbstractArray}
-    local buf, shape, retrs
     sz = size(ret)
-    first = true
     for dim in 1:ndims(ret)
         it = iextract(its, dim)
         if it != NoInterp

src/extrapolation/filled.jl

@@ -30,7 +30,7 @@ extrapolate(itp::AbstractInterpolation{T,N,IT}, fillvalue) where {T,N,IT} = Fill
 @inline function (etp::FilledExtrapolation{T,N})(x::Vararg{Number,N}) where {T,N}
     itp = parent(etp)
     coefs = coefficients(itp)
-    Tret = typeof(prod(x) * zero(eltype(coefs)))
+    Tret = promote_type(eltype(coefs), typeof.(x)...)
     if checkbounds(Bool, itp, x...)
         @inbounds itp(x...)
     else

src/gridded/gridded.jl

@@ -28,7 +28,6 @@ end
 
 function GriddedInterpolation(::Type{TWeights}, knots::NTuple{N,GridIndex}, A::AbstractArray{TCoefs,N}, it::IT) where {N,TCoefs,TWeights<:Real,IT<:DimSpec{Gridded},pad}
     isconcretetype(IT) || error("The b-spline type must be a leaf type (was $IT)")
-    isconcretetype(TCoefs) || @warn("For performance reasons, consider using an array of a concrete type (eltype(A) == $(eltype(A)))")
 
     check_gridded(it, knots, axes(A))
     c = zero(TWeights)

src/gridded/indexing.jl

@@ -40,6 +40,14 @@ function gridded_floorbounds(x, knotvec::AbstractVector)
     max(i, first(axes1(knotvec)))
 end
 
+ceilbounds(x::Integer, knotvec::AbstractVector) = gridded_ceilbounds(x, knotvec)
+ceilbounds(x, knotvec::AbstractVector) = gridded_ceilbounds(x, knotvec)
+function gridded_ceilbounds(x, knotvec::AbstractVector)
+    i = find_knot_index(knotvec, x)
+    iclamp = max(i, first(axes1(knotvec)))
+    min(iclamp+1, last(axes1(knotvec)))
+end
+
 @inline find_knot_index(knotv, x) = searchsortedfirst(knotv, x, Base.Order.ForwardOrdering()) - 1
 @inline find_knot_index(knotv, x::AbstractVector) = searchsortedfirst_vec(knotv, x) .- 1