Implement remaining boundary conditions for linear interpolation

This also:
- compares against the correct answer
- provides a number of fixes to many of the boundary conditions
- add a BCinbounds type that avoids any bounds-checking

src/Grid.jl

@@ -28,6 +28,7 @@ export
     BCperiodic,
     BCnearest,
     BCfill,
+    BCinbounds,
     Counter,
     AbstractInterpGrid,
     InterpGrid,

src/boundaryconditions.jl

@@ -10,6 +10,7 @@ type BCreflect <: BoundaryCondition; end # Reflecting boundary conditions
 type BCperiodic <: BoundaryCondition; end # Periodic boundary conditions
 type BCnearest <: BoundaryCondition; end # Return closest edge element
 type BCfill <: BoundaryCondition; end # Use specified fill value
+type BCinbounds <: BoundaryCondition; end # user guarantees that no out-of-bounds values will be accessed
 
 # Note: for interpolation, BCna is currently defined to be identical
 # to BCnan. Other applications might define different behavior,

src/newinterp.jl

@@ -29,9 +29,9 @@ size(G::InterpGridNew) = size(G.coefs)
 function body_gen(::Type{BCnil}, ::Type{InterpLinear}, N::Integer)
     ex = interplinear_gen(N)
     quote
-        @nexprs $N d->(1 <= x_d < size(G,d) || throw(BoundsError()))
+        @nexprs $N d->(1 <= x_d <= size(G,d) || throw(BoundsError()))
         @nexprs $N d->(ix_d = ifloor(x_d); fx_d = x_d - convert(typeof(x_d), ix_d))
-        @nexprs $N d->(ixp_d = ix_d + 1)
+        @nexprs $N d->(ixp_d = ix_d == size(G,d) ? ix_d : ix_d + 1)
         A = G.coefs
         @inbounds ret = $ex
         ret
@@ -41,9 +41,9 @@ end
 function body_gen(::Type{BCnan}, ::Type{InterpLinear}, N::Integer)
     ex = interplinear_gen(N)
     quote
-        @nexprs $N d->(1 <= x_d < size(G,d) || return(nan(eltype(G))))
+        @nexprs $N d->(1 <= x_d <= size(G,d) || return(nan(eltype(G))))
         @nexprs $N d->(ix_d = ifloor(x_d); fx_d = x_d - convert(typeof(x_d), ix_d))
-        @nexprs $N d->(ixp_d = ix_d + 1)
+        @nexprs $N d->(ixp_d = fx_d == 0 ? ix_d : ix_d + 1)
         A = G.coefs
         @inbounds ret = $ex
         ret
@@ -68,12 +68,36 @@ function body_gen(::Type{BCna}, ::Type{InterpLinear}, N::Integer)
     end
 end
 
+# mod is slow, so this goes to some effort to avoid two calls to mod
+function body_gen(::Type{BCreflect}, ::Type{InterpLinear}, N::Integer)
+    ex = interplinear_gen(N)
+    quote
+        @nexprs $N d->(ix_d = ifloor(x_d); fx_d = x_d - convert(typeof(x_d), ix_d))
+        @nexprs $N d->( len = size(G,d);
+                        if !(1 <= ix_d <= len)
+                            ix_d = mod(ix_d-1, 2len)
+                            if ix_d < len
+                                ix_d = ix_d+1
+                                ixp_d = ix_d < len ? ix_d+1 : len
+                            else
+                                ix_d = 2len-ix_d
+                                ixp_d = ix_d > 1 ? ix_d-1 : 1
+                            end
+                        else
+                            ixp_d = ix_d == len ? len : ix_d + 1
+                        end)
+        A = G.coefs
+        @inbounds ret = $ex
+        ret
+    end
+end
+
 function body_gen(::Type{BCperiodic}, ::Type{InterpLinear}, N::Integer)
     ex = interplinear_gen(N)
     quote
         @nexprs $N d->(ix_d = ifloor(x_d); fx_d = x_d - convert(typeof(x_d), ix_d))
-        @nexprs $N d->(ixp_d = (ix_d % size(G,d)) + 1)
-        @nexprs $N d->(ix_d = ((ix_d - 1) % size(G,d)) + 1)
+        @nexprs $N d->(ix_d = mod(ix_d-1, size(G,d)) + 1)
+        @nexprs $N d->(ixp_d = ix_d < size(G,d) ? ix_d+1 : 1)
         A = G.coefs
         @inbounds ret = $ex
         ret
@@ -92,6 +116,30 @@ function body_gen(::Type{BCnearest}, ::Type{InterpLinear}, N::Integer)
     end
 end
 
+function body_gen(::Type{BCfill}, ::Type{InterpLinear}, N::Integer)
+    ex = interplinear_gen(N)
+    quote
+        @nexprs $N d->(1 <= x_d <= size(G,d) || return(G.fillval))
+        @nexprs $N d->(ix_d = ifloor(x_d); fx_d = x_d - convert(typeof(x_d), ix_d))
+        @nexprs $N d->(ixp_d = ix_d == size(G,d) ? ix_d : ix_d + 1)
+        A = G.coefs
+        @inbounds ret = $ex
+        ret
+    end
+end
+
+
+function body_gen(::Type{BCinbounds}, ::Type{InterpLinear}, N::Integer)
+    ex = interplinear_gen(N)
+    quote
+        @nexprs $N d->(ix_d = ifloor(x_d); fx_d = x_d - convert(typeof(x_d), ix_d))
+        ixp_d = ix_d + 1
+        A = G.coefs
+        @inbounds ret = $ex
+        ret
+    end
+end
+
 # Here's the function that does the indexing magic. It works recursively.
 # Try this:
 #   interplinear_gen(1)
@@ -110,6 +158,7 @@ function interplinear_gen(N::Integer, offsets...)
     end
 end
 
+
 # For type inference
 promote_type_grid(T, x...) = promote_type(T, typeof(x)...)
 
@@ -121,7 +170,7 @@ promote_type_grid(T, x...) = promote_type(T, typeof(x)...)
 # This creates getindex
 for IT in (InterpLinear,)
     # for BC in subtypes(BoundaryCondition)
-    for BC in (BCnil, BCnan, BCna, BCperiodic, BCnearest)
+    for BC in subtypes(BoundaryCondition)
         eval(ngenerate(:N, :(promote_type_grid(T, x...)), :(getindex{T,N}(G::InterpGridNew{T,N,$BC,$IT}, x::NTuple{N,Real}...)),
                       N->body_gen(BC, IT, N)))
     end

test/newinterp.jl

@@ -1,50 +1,59 @@
 using Grid, Base.Test
 include(Pkg.dir("Grid", "src/newinterp.jl"))
 
-function runinterp(n, G, x)
-    s = 0.0
-    for j = 1:n
-        for i = 1:length(x)
-            s += G[x[i]]
-        end
-    end
-    s
+## Construct ground-truth values for 1d interpolation
+A1 = float([1,2,3,4])
+xpos = [-1.5:0.1:6.5]
+inbounds = 1 .<= xpos .<= length(A1)
+A1inbounds = xpos[inbounds]
+A1nil = fill(Inf, length(xpos)); A1nil[inbounds] = A1inbounds  # use Inf as a sentinel for BoundsError
+A1nan = fill(NaN, length(xpos)); A1nan[inbounds] = A1inbounds
+A1na = fill(NaN, length(xpos)); A1na[inbounds] = A1inbounds; A1na[0 .< xpos .< 1] = 1; A1na[4 .< xpos .< 5] = 4
+A1reflect = zeros(length(xpos)); A1reflect[inbounds] = A1inbounds; A1reflect[0 .<= xpos .<= 1] = 1; A1reflect[4 .<= xpos .<= 5] = 4
+l = findfirst(xpos .>= 0); A1reflect[1:l] = A1inbounds[l:-1:1]; l = findfirst(xpos .>= 5); A1reflect[l:end] = A1inbounds[end:-1:length(A1reflect)-l+1]
+xirange = ifloor(minimum(xpos))-1:iceil(maximum(xpos))+1
+Aiper = A1[Int[mod(x-1,length(A1))+1 for x in xirange]]
+A1periodic = similar(A1na)
+for i = 1:length(xpos)
+    x = xpos[i]
+    ix = ifloor(x)
+    j = find(xirange .== ix)[1]
+    fx = x - ix
+    A1periodic[i] = (1-fx)*Aiper[j]+fx*Aiper[j+1]
 end
-
-function runinterp(n, G, x, y)
-    s = 0.0
-    for j = 1:n
-        for i = 1:length(x)
-            s += G[x[i], y[i]]
-        end
-    end
-    s
+A1nearest = zeros(length(xpos)); A1nearest[inbounds] = A1inbounds; A1nearest[xpos .< 1] = 1; A1nearest[4 .< xpos] = 4
+A1fill = zeros(length(xpos)); A1fill[inbounds] = A1inbounds
+# Plot the values, to make sure we've done it right
+if isdefined(Main, :Gadfly)
+    println("Plotting")
+    Aall = vcat(
+        DataFrame(x = xpos, y = A1nil, t = "BCnil"),
+        DataFrame(x = xpos, y = A1nan, t = "BCnan"),
+        DataFrame(x = xpos, y = A1na,  t = "BCna"),
+        DataFrame(x = xpos, y = A1reflect, t = "BCreflect"),
+        DataFrame(x = xpos, y = A1periodic, t = "BCperiodic"),
+        DataFrame(x = xpos, y = A1nearest, t = "BCnearest"),
+        DataFrame(x = xpos, y = A1fill, t = "BCfill"))
+    set_default_plot_size(30cm, 7.5cm)
+    plot(Aall, x = :x, y = :y, xgroup = :t, Geom.subplot_grid(Geom.line))
 end
 
-A1 = float([1,2,3,4])
-A2 = rand(4,4)
-x = rand(1.0:3.999, 10^6)
-y = rand(1.0:3.999, 10^6)
 
-for BC in (BCnil, BCnan, BCna, BCperiodic, BCnearest)
+for (BC,correct) in ((BCnil,A1nil), (BCnan,A1nan), (BCna,A1na), (BCreflect, A1reflect), (BCperiodic, A1periodic), (BCnearest, A1nearest), (BCfill,A1fill))
     println(BC)
-
-    G = InterpGrid(A1, BC, InterpLinear)
-    Gnew = InterpNew.InterpGridNew{Float64,1,BC,InterpLinear}(A1, 0.0);
-
-    runinterp(1, G, [1.1])
-    runinterp(1, Gnew, [1.1])
-    @time s = runinterp(10, G, x)
-    @time snew = runinterp(10, Gnew, x)
-    @test_approx_eq s snew
-
-
-    G = InterpGrid(A2, BC, InterpLinear)
-    Gnew = InterpNew.InterpGridNew{Float64,2,BC,InterpLinear}(A2, 0.0);
-
-    runinterp(1, G, [1.1], [1.2])
-    runinterp(1, Gnew, [1.1], [1.2])
-    @time s = runinterp(10, G, x, y)
-    @time snew = runinterp(10, Gnew, x, y)
-    @test_approx_eq s snew
+    G = InterpNew.InterpGridNew{Float64,1,BC,InterpLinear}(A1, 0.0);
+    for i = 1:length(xpos)
+        x = xpos[i]
+        y = correct[i]
+        if !isinf(y)
+            try
+                @test_approx_eq y G[x]
+            catch err
+                @show x, y
+                rethrow(err)
+            end
+        else
+            @test_throws BoundsError G[x]
+        end
+    end
 end