WIP: complete refactor

tomasaschan · tomasaschan · commit be88fca87cc1 · 2015-09-22T18:00:42.000+02:00
diff --git a/src/b-splines/b-splines.jl b/src/b-splines/b-splines.jl
@@ -83,11 +83,6 @@ end
 
 index_gen{IT}(::Type{IT}, N::Integer, offsets...) = index_gen(iextract(IT, min(length(offsets)+1, N)), IT, N, offsets...)
 
-@generated function gradient{T,N,TCoefs,IT,GT,pad}(itp::BSplineInterpolation{T,N,TCoefs,IT,GT,pad}, xs...)
-    n = count_interp_dims(IT, N)
-    :(gradient!(Array(T,$n), itp, xs...))
-end
-
 include("constant.jl")
 include("linear.jl")
 include("quadratic.jl")
diff --git a/src/b-splines/indexing.jl b/src/b-splines/indexing.jl
@@ -58,6 +58,14 @@ function gradient_impl{T,N,TCoefs,IT<:DimSpec{BSpline},GT<:DimSpec{GridType},Pad
 end
 
 
+@generated function gradient{T,N,TCoefs,IT,GT,pad}(itp::BSplineInterpolation{T,N,TCoefs,IT,GT,pad}, xs...)
+    n = count_interp_dims(IT, N)
+    quote
+        @which gradient!(Array(T, $n), itp, xs...)
+        gradient!(Array(T,$n), itp, xs...)
+    end
+end
+
 @generated function gradient!{T,N}(g::AbstractVector, itp::BSplineInterpolation{T,N}, xs::Number...)
     length(xs) == N || error("Can only be called with $N indexes")
     gradient_impl(itp)
diff --git a/src/extrapolation/constant.jl b/src/extrapolation/constant.jl
@@ -1,15 +1,14 @@
-type ConstantExtrapolation{T,N,ITP,IT,GT} <: AbstractExtrapolation{T,N,ITP,IT,GT}
-    itp::ITP
+function extrap_prep(N, ::Type{Tuple{Flat,Flat}})
+    :(@nexprs $N d->(xs_d = clamp(xs_d, lbound(etp,d), ubound(etp,d))))
 end
-ConstantExtrapolation{T,ITP,IT,GT}(::Type{T}, N, itp::ITP, ::Type{IT}, ::Type{GT}) =
-    ConstantExtrapolation{T,N,ITP,IT,GT}(itp)
 
-extrapolate{T,N,IT,GT}(itp::AbstractInterpolation{T,N,IT,GT}, ::Type{Flat}) =
-    ConstantExtrapolation(T,N,itp,IT,GT)
-
-function extrap_prep{T,N,ITP,IT,GT}(etp::Type{ConstantExtrapolation{T,N,ITP,IT,GT}}, xs...)
-    :(@nexprs $N d->(xs_d = clamp(xs[d], lbound(etp,d), ubound(etp,d))))
+function extrap_prep{T}(N, ::Type{Tuple{Flat,T}})
+	quote
+		@nexprs $N d->(xs_d = max(xs_d, lbound(etp,d)))
+		$(extrap_prep_hi(N, T))
+	end
 end
 
-lbound(etp::ConstantExtrapolation, d) = lbound(etp.itp, d)
-ubound(etp::ConstantExtrapolation, d) = ubound(etp.itp, d)
+function extrap_prep_hi(N, ::Type{Flat})
+	:(@nexprs $N d->(xs_d = min(xs_d, ubound(etp,d))))
+end
diff --git a/src/extrapolation/error.jl b/src/extrapolation/error.jl
@@ -1,15 +1,16 @@
 immutable Throw end
-type ErrorExtrapolation{T,N,ITPT,IT,GT} <: AbstractExtrapolation{T,N,ITPT,IT,GT}
-    itp::IT
+
+function extrap_prep(N, ::Type{Tuple{Throw,Throw}})
+    :(@nexprs $N d->(lbound(etp, d) <= xs_d <= ubound(etp, d) || throw(BoundsError())))
 end
-ErrorExtrapolation{T,ITPT,IT,GT}(::Type{T}, N, itp::ITPT, ::Type{IT}, ::Type{GT}) =
-    ErrorExtrapolation{T,N,ITPT,IT,GT}(itp)
-extrapolate{T,N,IT,GT}(itp::AbstractInterpolation{T,N,IT,GT}, ::Type{Throw}) = 
-    ErrorExtrapolation(T,N,itp,IT,GT)
 
-function extrap_prep{T,N,ITPT,IT,GT}(etp::Type{ErrorExtrapolation{T,N,ITPT,IT,GT}}, xs...)
-    :(@nexprs $N d->(lbound(etp,d) <= xs[d] <= ubound(etp,d) || throw(BoundsError())))
+function extrap_prep{T}(N, ::Type{Tuple{Throw,T}})
+	quote
+		@nexprs $N d->(lbound(etp, d) <= xs_d || throw(BoundsError()))
+		extrap_prep_hi(T, xs...)
+	end
 end
 
-lbound(etp::ErrorExtrapolation, d) = lbound(etp.itp, d)
-ubound(etp::ErrorExtrapolation, d) = ubound(etp.itp, d)
+function extrap_prep_hi(N, ::Type{Throw})
+	:(@nexprs $N d->(xs_d <= ubound(etp, d) || throw(BoundsError())))
+end
diff --git a/src/extrapolation/extrapolation.jl b/src/extrapolation/extrapolation.jl
@@ -1,11 +1,41 @@
 export Throw,
        FilledExtrapolation   # for direct control over typeof(fillvalue)
 
+type Extrapolation{T,N,ITPT,IT,GT,ET} <: AbstractInterpolationWrapper{T,N,ITPT,IT,GT}
+	itp::ITPT
+end
+Extrapolation{T,ITPT,IT,GT,ET}(::Type{T}, N, itp::ITPT, ::Type{IT}, ::Type{GT}, ::Type{ET}) =
+	Extrapolation{T,N,ITPT,IT,GT,ET}(itp)
+extrapolate{T,N,IT,GT,ET}(itp::AbstractInterpolation{T,N,IT,GT}, ::Type{ET}) =
+	Extrapolation(T,N,itp,IT,GT,ET)
+
+extrap_prep{T}(N, ::Type{T}) = extrap_prep(N, Tuple{T,T})
+
+
 include("error.jl")
 include("constant.jl")
 include("linear.jl")
 include("reflect.jl")
 include("periodic.jl")
-include("indexing.jl")
+
+@generated function getindex{T,N,ITPT,IT,GT,ET}(etp::Extrapolation{T,N,ITPT,IT,GT,ET}, xs...)
+    quote
+    	@nexprs $N d->(xs_d = xs[d])
+        $(extrap_prep(N, ET))
+        itp = etp.itp
+        @nref $N itp xs
+    end
+end
+
+@generated function gradient!{T,N,ITPT,IT,GT,ET}(g::AbstractVector, etp::Extrapolation{T,N,ITPT,IT,GT,ET}, xs...)
+	quote
+		$(extrap_prep_gradient(N, ET, xs...))
+		gradient!(g)
+	end
+end
+
+lbound(etp::Extrapolation, d) = lbound(etp.itp, d)
+ubound(etp::Extrapolation, d) = ubound(etp.itp, d)
+size(etp::Extrapolation, d) = size(etp.itp, d)
 
 include("filled.jl")
diff --git a/src/extrapolation/indexing.jl b/src/extrapolation/indexing.jl
@@ -1,7 +0,0 @@
-@generated function getindex{T,N}(etp::AbstractExtrapolation{T,N}, xs...)
-    quote
-        $(extrap_prep(etp, xs...))
-        itp = etp.itp
-        @nref $N itp xs
-    end
-end
diff --git a/src/extrapolation/linear.jl b/src/extrapolation/linear.jl
@@ -1,43 +1,32 @@
-type LinearExtrapolation{T,N,ITP,IT,GT} <: AbstractExtrapolation{T,N,ITP,IT,GT}
-    itp::ITP
+function extrap_prep{T}(N, ::Type{Tuple{Linear,T}}, xs...)
+    quote
+        $(extrap_prep_lo(N, Linear, xs...))
+        $(extrap_prep_hi(N, T, xs...))
+    end
 end
-LinearExtrapolation{T,ITP,IT,GT}(::Type{T}, N, itp::ITP, ::Type{IT}, ::Type{GT}) =
-    LinearExtrapolation{T,N,ITP,IT,GT}(itp)
 
-extrapolate{T,N,IT,GT}(itp::AbstractInterpolation{T,N,IT,GT}, ::Type{Linear}) =
-    LinearExtrapolation(T,N,itp,IT,GT)
-
-function extrap_prep{T,N,ITP,IT}(etp::Type{LinearExtrapolation{T,N,ITP,IT,OnGrid}}, xs...)
+function extrap_prep_lo(N, ::Type{Linear})
+    coords = [symbol("xs_", k) for k in 1:N]
     quote
-        xc = collect(xs...)
         @nexprs $N d->begin
-            if xc[d] < 1
-                xc[d] = one(eltype(xc))
-                k = gradient(etp.itp, xc...)[d]
-                return etp.itp[xc...] + k * (xs[d] - one(typeof(xc[d])))
-            elseif xc[d] > size(etp.itp, d)
-                xc[d] = size(etp.itp, d)
-                k = gradient(etp.itp, xc...)[d]
-                return etp.itp[xc...] + k * (xs[d] - size(etp.itp, d))
+            if xs_d < lbound(etp.itp, d)
+                xs_d = lbound(etp.itp, d)
+                k = gradient(etp.itp, $(coords...))[d]
+                return (@nref $N etp xs) + k * (xs[d] - lbound(etp.itp, d))
             end
+            
         end
     end
 end
-
-function extrap_prep{T,N,ITP,IT}(etp::Type{LinearExtrapolation{T,N,ITP,IT,OnCell}}, xs...)
+function extrap_prep_hi(N, ::Type{Linear})
+    coords = [symbol("xs_", k) for k in 1:N]
     quote
-        xc = collect(xs...)
         @nexprs $N d->begin
-            if xc[d] < 1//2
-                xc[d] = convert(typeof(xc[d], 1//2))
-                k = gradient(etp.itp, xc...)[d]
-                return etp.itp[xc...] + k * (xs[d] - convert(typeof(xc[d], 1//2)))
-            elseif xc[d] > size(etp.itp, d) + 1//2
-                xc[d] = size(etp.itp, d) + convert(typeof(xc[d], 1//2))
-                k = gradient(etp.itp, xc...)[d]
-                return etp.itp[xc...] + k * (xs[d] - size(etp.itp, d) - convert(typeof(xc[d], 1//2)))
+            if xs_d > ubound(etp.itp, d)
+                xs_d = ubound(etp.itp, d)
+                k = gradient(etp.itp, $(coords...))[d]
+                return (@nref $N etp xs) + k * (xs[d] - ubound(etp.itp, d))
             end
         end
-        @nexprs $N d->(x_d = xs[d])
     end
 end
diff --git a/src/extrapolation/periodic.jl b/src/extrapolation/periodic.jl
@@ -1,12 +1,39 @@
-type PeriodicExtrapolation{T,N,ITP,IT,GT} <: AbstractExtrapolation{T,N,ITP,IT,GT}
-    itp::ITP
+function extrap_prep(N, ::Type{Tuple{Periodic, Periodic}})
+    :(@nexprs $N d->$(extrap_prep_dim(:d, Periodic)))
 end
-PeriodicExtrapolation{T,ITP,IT,GT}(::Type{T}, N, itp::ITP, ::Type{IT}, ::Type{GT}) =
-    PeriodicExtrapolation{T,N,ITP,IT,GT}(itp)
 
-extrapolate{T,N,IT,GT}(itp::AbstractInterpolation{T,N,IT,GT}, ::Type{Periodic}) =
-    PeriodicExtrapolation(T,N,itp,IT,GT)
+function extrap_prep_dim(d, ::Type{Periodic})
+    xs_d = symbol("xs_", d)
+    :($xs_d = mod(xs[$d] - lbound(etp.itp, $d), ubound(etp.itp, $d) - lbound(etp.itp, $d)) + lbound(etp.itp, $d))
+end
 
-function extrap_prep{T,N,ITP,IT,GT}(etp::Type{PeriodicExtrapolation{T,N,ITP,IT,GT}}, xs...)
-    :(@nexprs $N d->(x_d = mod(xs[d]-one($xs[d]), size(etp.itp,d)-one($xs[d]))+one($xs[d])))
+function extrap_prep{T}(N, ::Type{Tuple{Periodic, T}})
+    quote
+        @nexprs $N d -> begin
+            xs_d = xs[d]
+            if xs_d < lbound(etp.itp, d)
+                $(extrap_prep_dim(:d, Periodic))
+            end
+        end
+        $(extrap_prep_hi(N, T))
+    end
+end
+function extrap_prep_hi(N, ::Type{Periodic})
+    quote
+        @nexprs $N d -> begin
+            xs_d = xs[d]
+            if xs_d > ubound(etp.itp, d)
+                $(extrap_prep_dim(d, Periodic))
+            end
+        end
+    end
 end
+#=
+
+0 <= x - lbound(itp) < ubound(itp) - lbound(itp)
+
+=> N = ubound(itp) - lbound(itp)
+
+mod(x - lbound(itp), ubound(itp) - lbound(itp)) + lbound(itp) = 
+
+=#
diff --git a/src/extrapolation/reflect.jl b/src/extrapolation/reflect.jl
@@ -1,30 +1,25 @@
-type ReflectingExtrapolation{T,N,ITP,IT,GT} <: AbstractExtrapolation{T,N,ITP,IT,GT}
-    itp::ITP
+function extrap_prep(N, ::Type{Tuple{Reflect,Reflect}})
+    :(@nexprs $N d->$(extrap_prep_dim(:d, Reflect)))
 end
-ReflectingExtrapolation{T,ITP,IT,GT}(::Type{T}, N, itp::ITP, ::Type{IT}, ::Type{GT}) =
-    ReflectingExtrapolation{T,N,ITP,IT,GT}(itp)
 
-extrapolate{T,N,IT,GT}(itp::AbstractInterpolation{T,N,IT,GT}, ::Type{Reflect}) =
-    ReflectingExtrapolation(T,N,itp,IT,GT)
 
-function extrap_prep{T,N,ITP,IT,GT}(etp::Type{ReflectingExtrapolation{T,N,ITP,IT,GT}}, xs...)
-    # First, translate x into the domain over [1,2N-1) (OnGrid) or [.5,2N+.5) (OnCell)
-    # i.e. into twice the size of the domain.
-    # Next, if x is now in the upper part of this "double-domain", reflect over the middle
-    # to obtain a new value x' for which f(x') == f(x), but where x' is inside the domain
-    if GT == OnGrid
-        return quote
-            @nexprs $N d->begin
-                x_d = mod(xs[d] - one($xs[d]), 2 * size(etp.itp,d) - 2) + one($xs[d])
-                x_d > size(etp.itp,d) && (x_d = convert($xs[d], 2) * size(etp.itp,d) - x_d)
-            end
-        end
-    else # GT == OnCell
-        return quote
-            @nexprs $N d->begin
-                x_d = mod(xs[d] - 1//2, 2*size(etp.itp,d) + 1//2)
-                x_d > size(etp.itp,d) + 1//2 && (x = convert($xs[d],2) * size(etp.itp,d) + one($xs[d] - x_d))
-            end
-        end
+"""
+`extrap_prep_dim(d, ::Type{Reflect})`
+
+First, translate x into the domain over [lbound, 2(ubound-lbound)) i.e. into twice the size of the domain.
+Next, if x is now in the upper part of this ''double-domain´´, reflect over the middle to obtain a new value x' for which f(x') == f(x), but where x' is inside the domain
+"""
+function extrap_prep_dim(d, ::Type{Reflect})
+    xs_d = symbol("xs_", d)
+    quote
+        start = lbound(etp.itp, $d)
+        width = ubound(etp.itp, $d) - start
+
+        $xs_d = mod($xs_d - start, 2width) + start
+        $xs_d > start + width && (xs_d = start + width - $xs_d)
     end
 end
+
+function extrap_prep_hi(N, ::Type{Reflect})
+    :(@nexprs $N d->(xs_d > ubound(etp.itp, d) && $(extrap_prep_dim(:d, Periodic))))
+end
diff --git a/test/extrapolation/runtests.jl b/test/extrapolation/runtests.jl
@@ -1,6 +1,6 @@
 module ExtrapTests
 
-using Base.Test
+using Base.Test, DualNumbers
 using Interpolations
 
 f(x) = sin((x-3)*2pi/9 - 1)
@@ -39,4 +39,16 @@ k_hi = A[end] - A[end-1]
 x_hi = xmax + 5.7
 @test_approx_eq etpl[x_hi] A[end]+k_hi*(x_hi-xmax)
 
+
+xmax, ymax = 8,8
+g(x, y) = (x^2 + 3x - 8) * (-2y^2 + y + 1)
+
+itp2g = interpolate(Float64[g(x,y) for x in 1:xmax, y in 1:ymax], Tuple{BSpline(Quadratic(Free)), BSpline(Linear)}, OnGrid)
+etp2g = extrapolate(itp2g, Tuple{Linear, Flat})
+
+@test_approx_eq @inferred(getindex(etp2g, -.5, 4)) itp2g[1,4]-1.5*epsilon(g(dual(1,1),4))
+@test_approx_eq @inferred(getindex(etp2g, 5, 100)) itp2g[5,ymax]
+
 end
+
+include("type-stability.jl")
diff --git a/test/extrapolation/type-stability.jl b/test/extrapolation/type-stability.jl
@@ -0,0 +1,68 @@
+module ExtrapTypeStability
+
+using Base.Test, Interpolations, DualNumbers
+
+# Test type-stability of 1-dimensional extrapolation
+f(x) = sin((x-3)*2pi/9 - 1)
+xmax = 10
+A = Float64[f(x) for x in 1:xmax]
+itpg = interpolate(A, BSpline(Linear), OnGrid)
+
+schemes = (Flat, Linear, Reflect, Periodic)
+
+
+for E in schemes
+    etp1 = extrapolate(itpg, E)
+    println(E)
+    
+    # In-bounds evaluation
+    @inferred(getindex(etp1, 3.4))
+    @inferred(getindex(etp1, 3))
+    @inferred(getindex(etp1, dual(3,1)))
+
+    # Out-of-bounds evaluation
+    @inferred(getindex(etp1, -3.4))
+    @inferred(getindex(etp1, -3))
+    @inferred(getindex(etp1, dual(-3,1)))
+
+    @inferred(getindex(etp1, 13.4))
+    @inferred(getindex(etp1, 13))
+    @inferred(getindex(etp1, dual(13,1)))
+end
+
+# Test type-stability of 2-dimensional extrapolation with homogeneous scheme
+g(y) = (y/100)^3
+ymax = 4
+A = Float64[f(x)*g(y) for x in 1:xmax, y in 1:ymax]
+itp2 = interpolate(A, BSpline(Linear), OnGrid)
+
+for E in schemes
+    etp2 = extrapolate(itp2, E)
+
+    # In-bounds evaluation
+    @inferred(getindex(etp2, 3.4, 2.1))
+    @inferred(getindex(etp2, 3,2))
+    @inferred(getindex(etp2, dual(3,1), 2.3))
+
+    # Out-of-bounds evaluation
+    @inferred(getindex(etp2, -3.4, 2.1))
+    @inferred(getindex(etp2, 3.4, -2.1))
+    @inferred(getindex(etp2, -3, 2))
+    @inferred(getindex(etp2, 3, -2))
+    @inferred(getindex(etp2, dual(-3,1), 2.3))
+    @inferred(getindex(etp2, dual(3,1), -2.3))
+    @inferred(getindex(etp2, -3.1, dual(-2.3,0)))
+    @inferred(getindex(etp2, 3,1, dual(-2.3,1)))
+
+    @inferred(getindex(etp2, 13.4,2.1))
+    @inferred(getindex(etp2, 13,2))
+    @inferred(getindex(etp2, dual(13,1),2.1))
+    @inferred(getindex(etp2, 3.4,12.1))
+    @inferred(getindex(etp2, 3,12))
+    @inferred(getindex(etp2, dual(3,1),12.1))
+    @inferred(getindex(etp2, dual(13,1),2.1))
+    @inferred(getindex(etp2, 3.1,dual(12.1,1)))
+    @inferred(getindex(etp2, 13,dual(12,1)))
+end
+
+end
