Merge pull request #8 from tlycken/constant-extrap

WIP features (so far: Flat BC, Constant extrap)

src/Interpolations.jl

@@ -11,9 +11,11 @@ export
     Quadratic,
     ExtrapError,
     ExtrapNaN,
+    ExtrapConstant,
     OnCell,
     OnGrid,
-    ExtendInner
+    ExtendInner,
+    Flat
 
 abstract Degree{N}
 
@@ -24,6 +26,7 @@ type OnCell <: GridRepresentation end
 abstract BoundaryCondition
 type None <: BoundaryCondition end
 type ExtendInner <: BoundaryCondition end
+type Flat <: BoundaryCondition end
 
 abstract InterpolationType{D<:Degree,BC<:BoundaryCondition,GR<:GridRepresentation}
 
@@ -63,8 +66,9 @@ include("quadratic.jl")
 
 
 # This creates getindex methods for all supported combinations
-for IT in (Constant{OnCell},Linear{OnGrid},Quadratic{ExtendInner,OnCell})
-    for EB in (ExtrapError,ExtrapNaN)
+for IT in (Constant{OnCell},Linear{OnGrid},Quadratic{ExtendInner,OnCell},Quadratic{Flat,OnCell})
+    for EB in (ExtrapError,ExtrapNaN,ExtrapConstant)
+
         eval(:(function getindex{T}(itp::Interpolation{T,1,$IT,$EB}, x::Real, d)
             d == 1 || throw(BoundsError())
             itp[x]
@@ -90,10 +94,10 @@ for IT in (Constant{OnCell},Linear{OnGrid},Quadratic{ExtendInner,OnCell})
                 # indices calculates the indices required for this interpolation degree,
                 # based on ix_d defined by bc_gen(), as well as the distance fx_d from 
                 # the cell index ix_d to the interpolation coordinate x_d
-                $(indices(degree(it), N))
+                $(indices(it, N))
 
                 # These coefficients determine the interpolation basis expansion
-                $(coefficients(degree(it), N))
+                $(coefficients(it, N))
 
                 @inbounds ret = $(index_gen(degree(it), N))
                 ret
@@ -103,7 +107,7 @@ for IT in (Constant{OnCell},Linear{OnGrid},Quadratic{ExtendInner,OnCell})
 end
 
 # This creates prefilter specializations for all interpolation types that need them
-for IT in (Quadratic{ExtendInner,OnCell},)
+for IT in (Quadratic{ExtendInner,OnCell},Quadratic{Flat,OnCell})
     @ngenerate N promote_type_grid(T, x...) function prefilter{T,N}(A::Array{T,N},it::IT)
         ret = similar(A)
         szs = collect(size(A))

src/constant.jl

@@ -9,13 +9,13 @@ function bc_gen{IT<:Constant}(::IT, N)
     end
 end
 
-function indices(::ConstantDegree, N)
+function indices(::Constant, N)
     quote
         # Constant interpolation doesn't need an fx_d
     end
 end
 
-function coefficients(::ConstantDegree, N)
+function coefficients(::Constant, N)
     quote
         @nexprs $N d->(c_d = one(typeof(x_d)))
     end

src/extrapolation.jl

@@ -17,3 +17,11 @@ function extrap_gen(::OnGrid, ::ExtrapNaN, N)
     end
 end
 extrap_gen(::OnCell, e::ExtrapNaN, N) = extrap_gen(OnGrid(), e, N)
+
+type ExtrapConstant <: ExtrapolationBehavior end
+function extrap_gen(::OnGrid, ::ExtrapConstant, N)
+    quote
+        @nexprs $N d->(x_d = clamp(x_d, 1, size(itp,d)))
+    end
+end
+extrap_gen(::OnCell, e::ExtrapConstant, N) = extrap_gen(OnGrid(), e, N)

src/linear.jl

@@ -12,7 +12,7 @@ function bc_gen(::Linear{OnGrid}, N)
     end
 end
 
-function indices(::LinearDegree, N)
+function indices(::Linear{OnGrid}, N)
     quote
         # fx_d is a parameter in [0,1] such that x_d = ix_d + fx_d
         @nexprs $N d->(fx_d = x_d - convert(typeof(x_d), ix_d))
@@ -21,7 +21,7 @@ function indices(::LinearDegree, N)
     end
 end
 
-function coefficients(::LinearDegree, N)
+function coefficients(::Linear{OnGrid}, N)
     quote
         @nexprs $N d->begin
             c_d = one(typeof(fx_d)) - fx_d

src/quadratic.jl

@@ -5,15 +5,16 @@ Quadratic{BC<:BoundaryCondition,GR<:GridRepresentation}(::BC, ::GR) = Quadratic{
 
 function bc_gen(::Quadratic{ExtendInner,OnCell}, N)
     quote
-        # We've already done extrapolation, but BC.ExtendInner mandates that
-        # the two outermost cells are treated specially. Therefore, we need to
+        # We've already done extrapolation, so if we're here, we know that
+        # 1 <= x_d <= size(itp,d), but ExtendInner mandates that the two
+        # outermost cells are treated specially. Therefore, we need to
         # further clamp the indices of the spline coefficients by one cell in
         # each direction.
         @nexprs $N d->(ix_d = clamp(iround(x_d), 2, size(itp,d)-1))
     end
 end
 
-function indices(::QuadraticDegree, N)
+function indices(::Quadratic{ExtendInner,OnCell}, N)
     quote
         @nexprs $N d->begin
             ixp_d = ix_d + 1
@@ -25,7 +26,7 @@ function indices(::QuadraticDegree, N)
     end
 end
 
-function coefficients(::QuadraticDegree, N)
+function coefficients(::Quadratic{ExtendInner,OnCell}, N)
     quote
         @nexprs $N d->begin
             cm_d = (fx_d-.5)^2 / 2
@@ -35,6 +36,44 @@ function coefficients(::QuadraticDegree, N)
     end
 end
 
+function bc_gen(::Quadratic{Flat,OnCell}, N)
+    quote
+        # After extrapolation has been handled, 1 <= x_d <= size(itp,d)
+        # The index is simply the closest integer.
+         @nexprs $N d->(ix_d = iround(x_d))
+        # end
+    end
+end
+
+function indices(::Quadratic{Flat,OnCell}, N)
+    quote
+        @nexprs $N d->begin
+            ixp_d = ix_d + 1
+            ixm_d = ix_d - 1
+
+            fx_d = x_d - convert(typeof(x_d), ix_d)
+
+            if ix_d == size(itp,d)
+                ixp_d = ixm_d
+            end
+            if ix_d == 1
+                ixm_d = ixp_d
+            end
+        end
+    end
+end
+
+
+function coefficients(::Quadratic{Flat,OnCell}, N)
+    quote
+        @nexprs $N d->begin
+            cm_d = (fx_d-.5)^2 / 2
+            c_d = .75 - fx_d^2
+            cp_d = (fx_d+.5)^2 / 2
+        end
+    end
+end
+
 # This assumes integral values ixm_d, ix_d, and ixp_d (typically ixm_d = ix_d-1, ixp_d = ix_d+1, except at boundaries),
 # coefficients cm_d, c_d, and cp_d, and an array itp.coefs
 function index_gen(degree::QuadraticDegree, N::Integer, offsets...)
@@ -49,12 +88,15 @@ function index_gen(degree::QuadraticDegree, N::Integer, offsets...)
     end
 end
 
-function prefiltering_system_matrix{T}(::Type{T}, n::Int, ::Quadratic{ExtendInner,OnCell})
+function unmodified_system_matrix{T}(::Type{T}, n::Int, ::Quadratic)
     du = fill(convert(T,1/8), n-1)
     d = fill(convert(T,3/4),n)
     dl = copy(du)
+    (dl,d,du)
+end
 
-    MT = lufact!(Tridiagonal(dl,d,du))
+function prefiltering_system_matrix{T}(::Type{T}, n::Int, q::Quadratic{ExtendInner,OnCell})
+    MT = lufact!(Tridiagonal(unmodified_system_matrix(T,n,q)...))
     U = zeros(T,n,2)
     V = zeros(T,2,n)
     C = zeros(T,2,2)
@@ -65,3 +107,10 @@ function prefiltering_system_matrix{T}(::Type{T}, n::Int, ::Quadratic{ExtendInne
 
     Woodbury(MT, U, C, V)
 end
+
+function prefiltering_system_matrix{T}(::Type{T}, n::Int, q::Quadratic{Flat,OnCell})
+    dl,d,du = unmodified_system_matrix(T,n,q)
+    du[1] += 1/8
+    dl[end] += 1/8
+    lufact!(Tridiagonal(dl, d, du))
+end

test/linear.jl

@@ -31,6 +31,17 @@ xlo, xhi = itp2[.9], itp2[xmax+.2]
 @test isnan(xlo)
 @test isnan(xhi)
 
+## ExtrapConstant
+
+itp3 = Interpolation(A, Linear(OnGrid()), ExtrapConstant())
+for x in [3.1:.2:4.3]
+    @test_approx_eq_eps f(x) itp3[x] abs(.1*f(x))
+end
+
+xlo, xhi = itp3[.9], itp3[xmax+.2]
+@test xlo == A[1]
+@test xhi == A[end]
+
 end
 
 module Linear2DTests

test/quadratic.jl

@@ -31,11 +31,33 @@ xlo, xhi = itp2[.9], itp2[xmax+.2]
 @test isnan(xlo)
 @test isnan(xhi)
 
+# Flat/ExtrapConstant
+
+itp3 = Interpolation(A, Quadratic(Flat(),OnCell()), ExtrapConstant())
+
+# Check inbounds and extrap values
+
+for x in [3.1:.2:4.3]
+    @test_approx_eq_eps f(x) itp1[x] abs(.1*f(x))
+end
+
+xlo, xhi = itp3[.9], itp3[xmax+.2]
+@test xlo == A[1]
+@test xhi == A[end]
+
+# Check continuity
+xs = [0:.1:length(A)+1]
+
+for i in 1:length(xs)-1
+    @test_approx_eq_eps itp3[xs[i]] itp3[xs[i+1]] .1
+end
+
 # Values on all data points except edges for ExtendInner
 
 for x in 2:xmax-1
     @test_approx_eq A[x] itp1[x]
     @test_approx_eq A[x] itp2[x]
+    @test_approx_eq A[x] itp3[x]
 end
 
 end