Merge pull request #428 from jaemolihm/periodic

Periodic boundary conditions for Linear and Constant

Author: mkitti

docs/src/extrapolation.md

@@ -12,3 +12,28 @@ itp = interpolate(1:7, BSpline(Linear()))
 etpf = extrapolate(itp, Flat())   # gives 1 on the left edge and 7 on the right edge
 etp0 = extrapolate(itp, 0)        # gives 0 everywhere outside [1,7]
 ```
+
+### Periodic extrapolation
+
+For uniformly sampled periodic data, one can perform periodic extrapolation for all types of
+B-Spline interpolations. By using the `Periodic(OnCell())` boundary condition in `interpolate`,
+one does not need to include the periodic image of the starting sample point.
+
+Examples:
+
+```julia
+f(x) = sin((x-3)*2pi/7 - 1)
+A = Float64[f(x) for x in 1:7] # Does not include the periodic image
+
+# Constant(Periodic())) is an alias for Constant{Nearest}(Periodic(OnCell()))
+itp0 = interpolate(A, BSpline(Constant(Periodic())))
+# Linear(Periodic())) is an alias for Linear(Periodic(OnCell()))
+itp1 = interpolate(A, BSpline(Linear(Periodic())))
+itp2 = interpolate(A, BSpline(Quadratic(Periodic(OnCell()))))
+itp3 = interpolate(A, BSpline(Cubic(Periodic(OnCell()))))
+
+etp0 = extrapolate(itp0, Periodic())
+etp1 = extrapolate(itp1, Periodic())
+etp2 = extrapolate(itp2, Periodic())
+etp3 = extrapolate(itp3, Periodic())
+```

src/b-splines/b-splines.jl

@@ -254,6 +254,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

@@ -5,8 +5,25 @@ struct Nearest <: ConstantInterpType end
 struct Previous <: ConstantInterpType end
 struct Next <: ConstantInterpType end
 
-struct Constant{T<:ConstantInterpType} <: Degree{0} end
-Constant() = Constant{Nearest}()
+struct Constant{T<:ConstantInterpType,BC<:Union{Throw{OnGrid},Periodic{OnCell}}} <: DegreeBC{0}
+    bc::BC
+end
+
+# Default to Nearest and Throw{OnGrid}
+Constant(args...) = Constant{Nearest}(args...)
+Constant{T}() where {T<:ConstantInterpType} = Constant{T,Throw{OnGrid}}(Throw(OnGrid()))
+Constant{T}(bc::BC) where {T<:ConstantInterpType,BC<:BoundaryCondition} = Constant{T,BC}(bc)
+Constant{T}(::Periodic{Nothing}) where {T<:ConstantInterpType} = Constant{T,Periodic{OnCell}}(Periodic(OnCell()))
+
+function Base.show(io::IO, deg::Constant)
+    print(io, nameof(typeof(deg)), '{', typeof(deg).parameters[1], '}', '(')
+    show(io, deg.bc)
+    print(io, ')')
+end
+
+function Base.show(io::IO, deg::Constant{T,Throw{OnGrid}}) where {T <: ConstantInterpType}
+    print(io, nameof(typeof(deg)), '{', typeof(deg).parameters[1], '}', '(', ')')
+end
 
 """
 Constant b-splines are *nearest-neighbor* interpolations, and effectively
@@ -30,6 +47,19 @@ function positions(c::Constant{Nearest}, ax, x)  # discontinuity occurs at half-
     fast_trunc(Int, xm), δx
 end
 
+function positions(c::Constant{Previous,Periodic{OnCell}}, ax, x)
+    # We do not use floorbounds because we do not want to add a half at
+    # the lowerbound to round up.
+    xm = floor(x)
+    δx = x - xm
+    modrange(fast_trunc(Int, xm), ax), δx
+end
+function positions(c::Constant{Next,Periodic{OnCell}}, ax, x)  # discontinuity occurs at integer locations
+    xm = ceilbounds(x, ax)
+    δx = x - xm
+    modrange(fast_trunc(Int, xm), ax), δx
+end
+
 value_weights(::Constant, δx) = (1,)
 gradient_weights(::Constant, δx) = (0,)
 hessian_weights(::Constant, δx) = (0,)

src/b-splines/linear.jl

@@ -1,4 +1,13 @@
-struct Linear <: Degree{1} end  # boundary conditions not supported
+struct Linear{BC<:Union{Throw{OnGrid},Periodic{OnCell}}} <: DegreeBC{1}
+    bc::BC
+end
+
+Linear() = Linear(Throw(OnGrid()))
+Linear(::Periodic{Nothing}) = Linear(Periodic(OnCell()))
+
+function Base.show(io::IO, deg::Linear{Throw{OnGrid}})
+    print(io, nameof(typeof(deg)), '(', ')')
+end
 
 """
     Linear()
@@ -27,16 +36,19 @@ a piecewise linear function connecting each pair of neighboring data points.
 """
 Linear
 
-function positions(::Linear, ax::AbstractUnitRange{<:Integer}, x)
+function positions(deg::Linear, ax::AbstractUnitRange{<:Integer}, x)
     f = floor(x)
     # When x == last(ax) we want to use the x-1, x pair
     f = ifelse(x == last(ax), f - oneunit(f), f)
     fi = fast_trunc(Int, f)
-    return fi, x-f
+    expand_index(deg, fi, ax), x-f
 end
+expand_index(::Linear{Throw{OnGrid}}, fi::Number, ax::AbstractUnitRange) = fi
+expand_index(::Linear{Periodic{OnCell}}, fi::Number, ax::AbstractUnitRange) =
+    (modrange(fi, ax), modrange(fi+1, ax))
 
 value_weights(::Linear, δx) = (1-δx, δx)
 gradient_weights(::Linear, δx) = (-oneunit(δx), oneunit(δx))
 hessian_weights(::Linear, δx) = (zero(δx), zero(δx))
 
-padded_axis(ax::AbstractUnitRange, ::BSpline{Linear}) = ax
+padded_axis(ax::AbstractUnitRange, ::BSpline{<:Linear}) = ax

test/b-splines/constant.jl

@@ -116,5 +116,55 @@
             end
         end
     end
-end
 
+    @testset "Constant periodic" begin
+        # Constructors
+        @test Constant() === Constant(Throw(OnGrid()))
+        @test Constant() isa Constant{Nearest,Throw{OnGrid}}
+        @test Constant(Periodic()) === Constant(Periodic(OnCell()))
+        @test Constant(Periodic()) isa Constant{Nearest,Periodic{OnCell}}
+        for T in (Nearest, Previous, Next)
+            it = Constant{T}()
+            @test it isa Constant{T, Throw{OnGrid}}
+            @test "$it" == "Constant{$T}()"
+            it = Constant{T}(Periodic())
+            @test it isa Constant{T, Periodic{OnCell}}
+            @test "$it" == "Constant{$T}(Periodic(OnCell()))"
+        end
+
+        for (constructor, copier) in ((interpolate, x -> x), (interpolate!, copy))
+            isinplace = constructor == interpolate!
+            itp_periodic = @inferred(constructor(copier(A1), BSpline(Constant(Periodic()))))
+            itp_previous = @inferred(constructor(copier(A1), BSpline(Constant{Previous}(Periodic()))))
+            itp_next = @inferred(constructor(copier(A1), BSpline(Constant{Next}(Periodic()))))
+
+            for itp in (itp_periodic, itp_previous, itp_next)
+                @test parent(itp) === itp.coefs
+                @test all(Interpolations.lbounds(itp) .≈ (0.5,))
+                @test all(Interpolations.ubounds(itp) .≈ (N1 + 0.5,))
+
+                check_axes(itp, A1, isinplace)
+                check_inbounds_values(itp, A1)
+                check_oob(itp)
+                can_eval_near_boundaries(itp)
+            end
+
+            # Evaluation between data points (tests constancy)
+            for i in 2:N1-1
+                @test A1[i] == itp_periodic(i + .3) == itp_periodic(i - .3)
+                @test A1[i] == itp_previous(i + .3) == itp_previous(i + .6)
+                @test A1[i - 1] == itp_previous(i - .3) == itp_previous(i - .6)
+                @test A1[i + 1] == itp_next(i + .3) == itp_next(i + .6)
+                @test A1[i] == itp_next(i - .3) == itp_next(i - .6)
+            end
+
+            # Evaluation between data points in [0.5, 1.5], [N1 - 0.5, N1 + 0.5].
+            @test A1[1] == itp_periodic(1 - .3) == itp_periodic(1 + .3)
+            @test A1[N1] == itp_periodic(N1 - .3) == itp_periodic(N1 + .3)
+            @test A1[1] == itp_previous(1 + .3)
+            @test A1[N1] == itp_previous(N1 + .3) == itp_previous(1 - .3)
+            @test A1[1] == itp_next(1 - .3) == itp_next(N1 + .3)
+            @test A1[N1] == itp_next(N1 - .3)
+        end
+    end
+end

test/b-splines/linear.jl

@@ -61,4 +61,39 @@
         etp = extrapolate(itp, Flat())
         @test sum(isnan.(etp)) == 0
     end
+
+    @testset "Linear periodic" begin
+        # Constructors
+        @test Linear() === Linear(Throw(OnGrid()))
+        @test Linear() isa Linear{Throw{OnGrid}}
+        @test Linear(Periodic()) === Linear(Periodic(OnCell()))
+        @test Linear(Periodic()) isa Linear{Periodic{OnCell}}
+
+        xmax = 10
+        f2(x) = sin((x-3)*2pi/xmax - 1) # Periodicity is xmax, not xmax-1
+        A1 = Float64[f2(x) for x in 1:xmax]
+
+        for (constructor, copier) in ((interpolate, identity), (interpolate!, copy))
+            isinplace = constructor == interpolate!
+            itp = @inferred(constructor(copier(A1), BSpline(Linear())))
+            itp_periodic = @inferred(constructor(copier(A1), BSpline(Linear(Periodic()))))
+
+            @test all(Interpolations.lbounds(itp_periodic) .≈ (0.5,))
+            @test all(Interpolations.ubounds(itp_periodic) .≈ (10.5,))
+
+            for x in 1:.2:xmax
+                @test itp(x) ≈ itp_periodic(x)
+            end
+
+            # Interpolation range for periodic should be [0.5, xmax + 0.5]
+            for x in 0.5:.1:0.9
+                @test f2(x) ≈ itp_periodic(x) atol=abs(0.1 * f2(x))
+                @test_throws BoundsError itp(x)
+            end
+            for x in xmax+0.1:.1:xmax+0.5
+                @test f2(x) ≈ itp_periodic(x) atol=abs(0.1 * f2(x))
+                @test_throws BoundsError itp(x)
+            end
+        end
+    end
 end

test/extrapolation/periodic.jl

@@ -0,0 +1,23 @@
+
+@testset "Periodic extrapolation" begin
+    xmax = 7
+    f(x) = sin((x-3)*2pi/xmax - 1)
+    A = Float64[f(x) for x in 1:xmax] # Does not include the periodic image
+
+    itp0 = interpolate(A, BSpline(Constant(Periodic())))
+    itp1 = interpolate(A, BSpline(Linear(Periodic())))
+    itp2 = interpolate(A, BSpline(Quadratic(Periodic(OnCell()))))
+    itp3 = interpolate(A, BSpline(Cubic(Periodic(OnCell()))))
+
+    etp0 = extrapolate(itp0, Periodic())
+    etp1 = extrapolate(itp1, Periodic())
+    etp2 = extrapolate(itp2, Periodic())
+    etp3 = extrapolate(itp3, Periodic())
+
+    for x in -xmax:.4:2*xmax
+        @test etp0(x) ≈ f(x) atol=0.5
+        @test etp1(x) ≈ f(x) atol=0.1
+        @test etp2(x) ≈ f(x) atol=0.01
+        @test etp3(x) ≈ f(x) atol=0.003
+    end
+end

test/extrapolation/runtests.jl

@@ -107,4 +107,5 @@ using Test
 
     include("type-stability.jl")
     include("non1.jl")
+    include("periodic.jl")
 end