WIP: metaprogramming-based interpolation

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/counter.jl

@@ -18,7 +18,7 @@ import Base: done, next, start
 # However, with sz appropriately defined, this version works for
 # arbitrary dimensions.
 
-type Counter
+immutable Counter
     max::Vector{Int}
 end
 Counter(sz::Tuple) = Counter(Int[sz...])

src/newinterp.jl

@@ -0,0 +1,179 @@
+module InterpNew
+
+using Grid
+import Grid: BoundaryCondition, InterpType
+import Base: getindex, size
+
+if VERSION.minor < 3
+    using Cartesian
+else
+    using Base.Cartesian
+end
+
+type InterpGridNew{T, N, BC<:BoundaryCondition, IT<:InterpType} <: AbstractArray{T,N}
+    coefs::Array{T,N}
+    fillval::T
+end
+size(G::InterpGridNew, d::Integer) = size(G.coefs, d)
+size(G::InterpGridNew) = size(G.coefs)
+
+# Create bodies like these:
+#    1 <= x_1 < size(G,1) || throw(BoundsError()))
+#    1 <= x_2 < size(G,2) || throw(BoundsError()))
+#    ix_1 = ifloor(x_1); fx_1 = x_1 - convert(typeof(x_1), ix_1)
+#    ix_2 = ifloor(x_2); fx_2 = x_2 - convert(typeof(x_2), ix_2)
+#    @inbounds ret =
+#       (1-fx_1)*((1-fx_2)*A[ix_1, ix_2] + fx_2*A[ix_1, ix_2+1]) +
+#       fx_1*((1-fx_2)*A[ix_1+1, ix_2] + fx_2*A[ix_1+1, ix_2+1])
+#    ret
+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->(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{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->(ix_d = ifloor(x_d); fx_d = x_d - convert(typeof(x_d), ix_d))
+        @nexprs $N d->(ixp_d = fx_d == 0 ? ix_d : ix_d + 1)
+        A = G.coefs
+        @inbounds ret = $ex
+        ret
+    end
+end
+
+function body_gen(::Type{BCna}, ::Type{InterpLinear}, N::Integer)
+    ex = interplinear_gen(N)
+    quote
+        @nexprs $N d->(0 < x_d < size(G,d) + 1 || 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->( if ix_d < 1
+                            ix_d, ixp_d = 1, 1
+                        elseif ix_d >= size(G,d)
+                            ix_d, ixp_d = size(G,d), size(G,d)
+                        else
+                            ixp_d = ix_d+1
+                        end)
+        A = G.coefs
+        @inbounds ret = $ex
+        ret
+    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->(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
+    end
+end
+
+function body_gen(::Type{BCnearest}, ::Type{InterpLinear}, N::Integer)
+    ex = interplinear_gen(N)
+    quote
+        @nexprs $N d->(x_d = clamp(x_d, 1, size(G,d)))
+        @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{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)
+#   interplinear_gen(2)
+# and you'll see what it does.
+# indexfunc = (dimension, offset) -> expression,
+#   for example  indexfunc(3, 0) returns  :ix_3
+#                indexfunc(3, 1) returns  :(ix_3 + 1)
+function interplinear_gen(N::Integer, offsets...)
+    if length(offsets) < N
+        sym = symbol("fx_"*string(length(offsets)+1))
+        return :((one($sym)-$sym)*$(interplinear_gen(N, offsets..., 0)) + $sym*$(interplinear_gen(N, offsets..., 1)))
+    else
+        indexes = [offsets[d] == 0 ? symbol("ix_"*string(d)) : symbol("ixp_"*string(d))  for d = 1:N]
+        return :(A[$(indexes...)])
+    end
+end
+
+
+# For type inference
+promote_type_grid(T, x...) = promote_type(T, typeof(x)...)
+
+# Because interpolation expressions are naturally generated by recursion, it's best to
+# have a body-generation function. That makes creating the functions a little more awkward,
+# but not too bad. If you want to see what this does, just copy-and-paste what's inside
+# the `eval` call at the command line (after importing all relevant names)
+
+# This creates getindex
+for IT in (InterpLinear,)
+    # for BC in subtypes(BoundaryCondition)
+    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
+end
+
+end

test/newinterp.jl

@@ -0,0 +1,59 @@
+using Grid, Base.Test
+include(Pkg.dir("Grid", "src/newinterp.jl"))
+
+## 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
+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
+
+
+for (BC,correct) in ((BCnil,A1nil), (BCnan,A1nan), (BCna,A1na), (BCreflect, A1reflect), (BCperiodic, A1periodic), (BCnearest, A1nearest), (BCfill,A1fill))
+    println(BC)
+    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