Merge pull request #154 from JuliaMath/teh/faster_etp

Faster extrapolation
JuliaMath · Apr 28, 2017 · 70e4128 · 70e4128
2 parents b86fad2 + 87964b1
commit 70e4128
Show file tree

Hide file tree

Showing 12 changed files with 74 additions and 39 deletions.
diff --git a/src/b-splines/b-splines.jl b/src/b-splines/b-splines.jl
@@ -38,15 +38,24 @@ padding(itp::AbstractInterpolation) = padding(typeof(itp))
 padextract(pad::Integer, d) = pad
 padextract(pad::Tuple{Vararg{Integer}}, d) = pad[d]
 
-lbound{T,N,TCoefs,IT}(itp::BSplineInterpolation{T,N,TCoefs,IT,OnGrid}, d) =
+lbound{T,N,TCoefs,IT}(itp::BSplineInterpolation{T,N,TCoefs,IT,OnGrid}, d::Integer) =
     first(indices(itp, d))
-ubound{T,N,TCoefs,IT}(itp::BSplineInterpolation{T,N,TCoefs,IT,OnGrid}, d) =
+ubound{T,N,TCoefs,IT}(itp::BSplineInterpolation{T,N,TCoefs,IT,OnGrid}, d::Integer) =
     last(indices(itp, d))
-lbound{T,N,TCoefs,IT}(itp::BSplineInterpolation{T,N,TCoefs,IT,OnCell}, d) =
+lbound{T,N,TCoefs,IT}(itp::BSplineInterpolation{T,N,TCoefs,IT,OnCell}, d::Integer) =
     first(indices(itp, d)) - 0.5
-ubound{T,N,TCoefs,IT}(itp::BSplineInterpolation{T,N,TCoefs,IT,OnCell}, d) =
+ubound{T,N,TCoefs,IT}(itp::BSplineInterpolation{T,N,TCoefs,IT,OnCell}, d::Integer) =
     last(indices(itp, d))+0.5
 
+lbound{T,N,TCoefs,IT}(itp::BSplineInterpolation{T,N,TCoefs,IT,OnGrid}, d, inds) =
+    first(inds)
+ubound{T,N,TCoefs,IT}(itp::BSplineInterpolation{T,N,TCoefs,IT,OnGrid}, d, inds) =
+    last(inds)
+lbound{T,N,TCoefs,IT}(itp::BSplineInterpolation{T,N,TCoefs,IT,OnCell}, d, inds) =
+    first(inds) - 0.5
+ubound{T,N,TCoefs,IT}(itp::BSplineInterpolation{T,N,TCoefs,IT,OnCell}, d, inds) =
+    last(inds)+0.5
+
 count_interp_dims{T,N,TCoefs,IT<:DimSpec{InterpolationType},GT<:DimSpec{GridType},pad}(::Type{BSplineInterpolation{T,N,TCoefs,IT,GT,pad}}, n) = count_interp_dims(IT, n)
 
 function size{T,N,TCoefs,IT,GT,pad}(itp::BSplineInterpolation{T,N,TCoefs,IT,GT,pad}, d)

diff --git a/src/extrapolation/extrap_prep.jl b/src/extrapolation/extrap_prep.jl
@@ -38,12 +38,16 @@ coordinate transformations, but for anything else methods for the low- and high-
 value cases need to be implemented for each scheme.
 """ extrap_prep
 
-extrap_prep{T}(::Type{T}, n::Val{1}) = extrap_prep(T, n, Val{1}())
+extrap_prep{T}(::Type{T}, n::Val{1}) = quote
+    inds_etp = indices(etp)
+    $(extrap_prep(T, n, Val{1}()))
+end
 extrap_prep{T}(::Type{Tuple{T}}, n::Val{1}) = extrap_prep(T, n)
 extrap_prep{T}(::Type{Tuple{T,T}}, n::Val{1}) = extrap_prep(T, n)
 extrap_prep{T}(::Type{Tuple{Tuple{T,T}}}, n::Val{1}) = extrap_prep(T, n)
 function extrap_prep{S,T}(::Type{Tuple{S,T}}, n::Val{1})
     quote
+        inds_etp = indices(etp)
         $(extrap_prep(S, n, Val{1}(), Val{:lo}()))
         $(extrap_prep(T, n, Val{1}(), Val{:hi}()))
     end
@@ -54,15 +58,15 @@ extrap_prep{S,T}(::Type{Tuple{Tuple{S,T}}}, n::Val{1}) = extrap_prep(Tuple{S,T},
 extrap_prep{T<:Tuple}(::Type{T}, ::Val{1}) = :(throw(ArgumentError("The 1-dimensional extrap configuration $T is not supported")))
 
 function extrap_prep{T,N}(::Type{T}, n::Val{N})
-    exprs = Expr[]
+    exprs = [:(inds_etp = indices(etp))]
     for d in 1:N
         push!(exprs, extrap_prep(T, n, Val{d}()))
     end
     return Expr(:block, exprs...)
 end
 function extrap_prep{N,T<:Tuple}(::Type{T}, n::Val{N})
     length(T.parameters) == N || return :(throw(ArgumentError("The $N-dimensional extrap configuration $T is not supported - must be a tuple of length $N (was length $(length(T.parameters)))")))
-    exprs = Expr[]
+    exprs = [:(inds_etp = indices(etp))]
     for d in 1:N
         Tdim = T.parameters[d]
         if Tdim <: Tuple

diff --git a/src/extrapolation/extrap_prep_gradient.jl b/src/extrapolation/extrap_prep_gradient.jl
@@ -1,11 +1,15 @@
 # See ?extrap_prep for documentation for all these methods
 
-extrap_prep{T}(g::Val{:gradient}, ::Type{T}, n::Val{1}) = extrap_prep(g, T, n, Val{1}())
+extrap_prep{T}(g::Val{:gradient}, ::Type{T}, n::Val{1}) = quote
+    inds_etp = indices(etp)
+    $(extrap_prep(g, T, n, Val{1}()))
+end
 extrap_prep{T}(g::Val{:gradient}, ::Type{Tuple{T}}, n::Val{1}) = extrap_prep(g, T, n)
 extrap_prep{T}(g::Val{:gradient}, ::Type{Tuple{T,T}}, n::Val{1}) = extrap_prep(g, T, n)
 extrap_prep{T}(g::Val{:gradient}, ::Type{Tuple{Tuple{T,T}}}, n::Val{1}) = extrap_prep(g, T, n)
 function extrap_prep{S,T}(g::Val{:gradient}, ::Type{Tuple{S,T}}, ::Val{1})
     quote
+        inds_etp = indices(etp)
         $(extrap_prep(g, S, n, Val{1}(), Val{:lo}()))
         $(extrap_prep(g, T, n, Val{1}(), Val{:hi}()))
     end
@@ -16,12 +20,12 @@ extrap_prep{T<:Tuple}(::Val{:gradient}, ::Type{T}, ::Val{1}) = :(throw(ArgumentE
 
 
 function extrap_prep{T,N}(g::Val{:gradient}, ::Type{T}, n::Val{N})
-    Expr(:block, [extrap_prep(g, T, n, Val{d}()) for d in 1:N]...)
+    Expr(:block, :(inds_etp = indices(etp)), [extrap_prep(g, T, n, Val{d}()) for d in 1:N]...)
 end
 
 function extrap_prep{T<:Tuple,N}(g::Val{:gradient}, ::Type{T}, n::Val{N})
     length(T.parameters) == N || return :(throw(ArgumentError("The $N-dimensional extrap configuration $T is not supported")))
-    exprs = Expr[]
+    exprs = [:(inds_etp = indices(etp))]
     for d in 1:N
         Tdim = T.parameters[d]
         if Tdim <: Tuple
@@ -47,4 +51,4 @@ function extrap_prep{S,T,N,d}(g::Val{:gradient}, ::Type{Tuple{S,T}}, n::Val{N},
 end
 
 extrap_prep{T,N,d}(g::Val{:gradient}, ::Type{T}, n::Val{N}, dim::Val{d}) = extrap_prep(g, Tuple{T,T}, n, dim)
-extrap_prep(g::Val{:gradient}, args...) = extrap_prep(args...)
+extrap_prep(g::Val{:gradient}, args...) = extrap_prep(args...)
diff --git a/src/extrapolation/extrapolation.jl b/src/extrapolation/extrapolation.jl
@@ -50,6 +50,7 @@ and indices. The heavy lifting is done by the `extrap_prep` function; see
 function getindex_impl{T,N,ITPT,IT,GT,ET}(etp::Type{Extrapolation{T,N,ITPT,IT,GT,ET}}, xs...)
     coords = [Symbol("xs_",d) for d in 1:N]
     quote
+        $(Expr(:meta, :inline))
         @nexprs $N d->(xs_d = xs[d])
         $(extrap_prep(ET, Val{N}()))
         etp.itp[$(coords...)]
@@ -66,6 +67,7 @@ checkbounds(::AbstractExtrapolation,I...) = nothing
 function gradient!_impl{T,N,ITPT,IT,GT,ET}(g, etp::Type{Extrapolation{T,N,ITPT,IT,GT,ET}}, xs...)
     coords = [Symbol("xs_", d) for d in 1:N]
     quote
+        $(Expr(:meta, :inline))
         @nexprs $N d->(xs_d = xs[d])
         $(extrap_prep(Val{:gradient}(), ET, Val{N}()))
         gradient!(g, etp.itp, $(coords...))
@@ -79,8 +81,10 @@ end
 
 lbound(etp::Extrapolation, d) = lbound(etp.itp, d)
 ubound(etp::Extrapolation, d) = ubound(etp.itp, d)
+lbound(etp::Extrapolation, d, inds) = lbound(etp.itp, d, inds)
+ubound(etp::Extrapolation, d, inds) = ubound(etp.itp, d, inds)
 size(etp::Extrapolation, d) = size(etp.itp, d)
-indices(etp::AbstractExtrapolation) = indices(etp.itp)
+@inline indices(etp::AbstractExtrapolation) = indices(etp.itp)
 indices(etp::AbstractExtrapolation, d) = indices(etp.itp, d)
 
 include("filled.jl")
diff --git a/src/extrapolation/filled.jl b/src/extrapolation/filled.jl
@@ -26,7 +26,8 @@ function getindex_impl{T,N,ITP,IT,GT,FT}(fitp::Type{FilledExtrapolation{T,N,ITP,
        $meta
        # Check to see if we're in the extrapolation region, i.e.,
        # out-of-bounds in an index
-       @nexprs $N d->((args[d] < lbound(fitp,d) || args[d] > ubound(fitp, d))) && return convert($Tret, fitp.fillvalue)::$Tret
+       inds_etp = indices(fitp)
+       @nexprs $N d->((args[d] < lbound(fitp, d, inds_etp[d]) || args[d] > ubound(fitp, d, inds_etp[d]))) && return convert($Tret, fitp.fillvalue)::$Tret
        # In the interpolation region
        return convert($Tret, getindex(fitp.itp,args...))::$Tret
    end
@@ -41,3 +42,5 @@ getindex{T}(fitp::FilledExtrapolation{T,1}, x::Number, y::Int) = y == 1 ? fitp[x
 
 lbound(etp::FilledExtrapolation, d) = lbound(etp.itp, d)
 ubound(etp::FilledExtrapolation, d) = ubound(etp.itp, d)
+lbound(etp::FilledExtrapolation, d, inds) = lbound(etp.itp, d, inds)
+ubound(etp::FilledExtrapolation, d, inds) = ubound(etp.itp, d, inds)
diff --git a/src/extrapolation/flat.jl b/src/extrapolation/flat.jl
@@ -1,25 +1,25 @@
 function extrap_prep{N,d}(::Type{Flat}, ::Val{N}, ::Val{d})
     xs_d = Symbol("xs_", d)
-    :($xs_d = clamp($xs_d, lbound(etp, $d), ubound(etp, $d)))
+    :($xs_d = clamp($xs_d, lbound(etp, $d, inds_etp[$d]), ubound(etp, $d, inds_etp[$d])))
 end
 
 function extrap_prep{N,d}(::Type{Flat}, ::Val{N}, ::Val{d}, ::Val{:lo})
     xs_d = Symbol("xs_", d)
-    :($xs_d = max($xs_d, lbound(etp, $d)))
+    :($xs_d = max($xs_d, lbound(etp, $d, inds_etp[$d])))
 end
 
 function extrap_prep{N,d}(::Type{Flat}, ::Val{N}, ::Val{d}, ::Val{:hi})
     xs_d = Symbol("xs_", d)
-    :($xs_d = min($xs_d, ubound(etp, $d)))
+    :($xs_d = min($xs_d, ubound(etp, $d, inds_etp[$d])))
 end
 
 function extrap_prep{N,d}(::Val{:gradient}, ::Type{Flat}, ::Val{N}, ::Val{d}, ::Val{:lo})
     coords = [Symbol("xs_", k) for k in 1:N]
     xs_d = coords[d]
 
     quote
-        if $xs_d < lbound(etp, $d)
-            $xs_d = lbound(etp, $d)
+        if $xs_d < lbound(etp, $d, inds_etp[$d])
+            $xs_d = lbound(etp, $d, inds_etp[$d])
             gradient!(g, etp.itp, $(coords...))
             g[$d] = 0
             return g
@@ -32,11 +32,11 @@ function extrap_prep{N,d}(::Val{:gradient}, ::Type{Flat}, ::Val{N}, ::Val{d}, ::
     xs_d = coords[d]
 
     quote
-        if $xs_d > ubound(etp, $d)
-            $xs_d = ubound(etp, $d)
+        if $xs_d > ubound(etp, $d, inds_etp[$d])
+            $xs_d = ubound(etp, $d, inds_etp[$d])
             gradient!(g, etp.itp, $(coords...))
             g[$d] = 0
             return g
         end
     end
-end
+end
diff --git a/src/extrapolation/linear.jl b/src/extrapolation/linear.jl
@@ -2,8 +2,8 @@ function extrap_prep{N,d}(::Type{Linear}, ::Val{N}, ::Val{d}, ::Val{:lo})
     coords = [Symbol("xs_", k) for k in 1:N]
     xs_d = coords[d]
     quote
-        if $xs_d < lbound(etp.itp, $d)
-            $xs_d = lbound(etp.itp, $d)
+        if $xs_d < lbound(etp.itp, $d, inds_etp[$d])
+            $xs_d = lbound(etp.itp, $d, inds_etp[$d])
             return etp[$(coords...)] + gradient(etp, $(coords...))[$d] * (xs[$d] - $xs_d)
         end
     end
@@ -12,8 +12,8 @@ function extrap_prep{N,d}(::Type{Linear}, ::Val{N}, ::Val{d}, ::Val{:hi})
     coords = [Symbol("xs_", k) for k in 1:N]
     xs_d = coords[d]
     quote
-        if $xs_d > ubound(etp, $d)
-            $xs_d = ubound(etp, $d)
+        if $xs_d > ubound(etp, $d, inds_etp[$d])
+            $xs_d = ubound(etp, $d, inds_etp[$d])
             return etp[$(coords...)] + gradient(etp, $(coords...))[$d] * (xs[$d] - $xs_d)
         end
     end
@@ -23,4 +23,4 @@ extrap_prep{N,d,l}(::Val{:gradient}, ::Type{Linear}, n::Val{N}, dim::Val{d}, loh
     extrap_prep(Flat, n, dim, lohi)
 
 extrap_prep{N,d}(::Val{:gradient}, ::Type{Linear}, n::Val{N}, dim::Val{d}) =
-    extrap_prep(Flat, n, dim)
+    extrap_prep(Flat, n, dim)
diff --git a/src/extrapolation/periodic.jl b/src/extrapolation/periodic.jl
@@ -5,23 +5,23 @@ Translate x into the domain [lbound, ubound] my means of `mod()`
 """
 function extrap_prep_dim(::Type{Periodic}, d)
     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))
+    :($xs_d = mod(xs[$d] - lbound(etp.itp, $d, inds_etp[$d]), ubound(etp.itp, $d, inds_etp[$d]) - lbound(etp.itp, $d, inds_etp[$d])) + lbound(etp.itp, $d, inds_etp[$d]))
 end
 
 extrap_prep{N,d}(::Type{Periodic}, ::Val{N}, ::Val{d}) = extrap_prep_dim(Periodic, d)
 function extrap_prep{N,d}(::Type{Periodic}, ::Val{N}, ::Val{d}, ::Val{:lo})
     xs_d = Symbol("xs_", d)
     quote
-        if $xs_d < lbound(etp.itp, $d)
+        if $xs_d < lbound(etp.itp, $d, inds_etp[$d])
             $(extrap_prep_dim(Periodic, d))
         end
     end
 end
 function extrap_prep{N,d}(::Type{Periodic}, ::Val{N}, ::Val{d}, ::Val{:hi})
     xs_d = Symbol("xs_", d)
     quote
-        if $xs_d > ubound(etp.itp, d)
+        if $xs_d > ubound(etp.itp, d, inds_etp[$d])
             $(extrap_prep_dim(Periodic, d))
         end
     end
-end
+end
diff --git a/src/extrapolation/reflect.jl b/src/extrapolation/reflect.jl
@@ -7,8 +7,8 @@ Next, if x is now in the upper part of this ''double-domain´´, reflect over th
 function extrap_prep_dim(::Type{Reflect}, d)
     xs_d = Symbol("xs_", d)
     quote
-        start = lbound(etp.itp, $d)
-        width = ubound(etp.itp, $d) - start
+        start = lbound(etp.itp, $d, inds_etp[$d])
+        width = ubound(etp.itp, $d, inds_etp[$d]) - start
 
         $xs_d = mod($xs_d - start, 2width) + start
         $xs_d > start + width && (xs_d = start + width - $xs_d)
@@ -18,9 +18,9 @@ end
 extrap_prep{N,d}(::Type{Reflect}, ::Val{N}, ::Val{d}) = extrap_prep_dim(Reflect, d)
 function extrap_prep{N,d}(::Type{Reflect}, ::Val{N}, ::Val{d}, ::Val{:lo})
     xs_d = Symbol("xs_", d)
-    :($xs_d < lbound(etp.itp, $d) && $(extrap_prep_dim(Reflect, d)))
+    :($xs_d < lbound(etp.itp, $d, inds_etp[$d]) && $(extrap_prep_dim(Reflect, d)))
 end
 function extrap_prep{N,d}(::Type{Reflect}, ::Val{N}, ::Val{d}, ::Val{:hi})
     xs_d = Symbol("xs_", d)
-    :($xs_d > ubound(etp.itp, $d) && $(extrap_prep_dim(Reflect, d)))
-end
+    :($xs_d > ubound(etp.itp, $d, inds_etp[$d]) && $(extrap_prep_dim(Reflect, d)))
+end
diff --git a/src/extrapolation/throw.jl b/src/extrapolation/throw.jl
@@ -1,14 +1,14 @@
 function extrap_prep{N,d}(::Type{Throw}, ::Val{N}, ::Val{d})
     xsym = Symbol("xs_", d)
-    :(lbound(etp, $d) <= $xsym <= ubound(etp, $d) || throw(BoundsError()))
+    :(lbound(etp, $d, inds_etp[$d]) <= $xsym <= ubound(etp, $d, inds_etp[$d]) || throw(BoundsError()))
 end
 
 function extrap_prep{N,d}(::Type{Throw}, ::Val{N}, ::Val{d}, ::Val{:lo})
     xsym = Symbol("xs_", d)
-    :(lbound(etp, $d) <= $xsym || throw(BoundsError()))
+    :(lbound(etp, $d, inds_etp[$d]) <= $xsym || throw(BoundsError()))
 end
 
 function extrap_prep{N,d}(::Type{Throw}, ::Val{N}, ::Val{d}, ::Val{:hi})
     xsym = Symbol("xs_", d)
-    :($xsym <= ubound(etp, $d) || throw(BoundsError()))
-end
+    :($xsym <= ubound(etp, $d, inds_etp[$d]) || throw(BoundsError()))
+end
diff --git a/src/gridded/gridded.jl b/src/gridded/gridded.jl
@@ -64,6 +64,8 @@ end
 
 lbound(itp::GriddedInterpolation, d) = itp.knots[d][1]
 ubound(itp::GriddedInterpolation, d) = itp.knots[d][end]
+lbound(itp::GriddedInterpolation, d, inds) = itp.knots[d][1]
+ubound(itp::GriddedInterpolation, d, inds) = itp.knots[d][end]
 
 include("constant.jl")
 include("linear.jl")

diff --git a/src/scaling/scaling.jl b/src/scaling/scaling.jl
@@ -50,6 +50,15 @@ lbound{T,N,ITPT,IT}(sitp::ScaledInterpolation{T,N,ITPT,IT,OnCell}, d) = 1 <= d <
 ubound{T,N,ITPT,IT}(sitp::ScaledInterpolation{T,N,ITPT,IT,OnGrid}, d) = 1 <= d <= N ? sitp.ranges[d][end] : throw(BoundsError())
 ubound{T,N,ITPT,IT}(sitp::ScaledInterpolation{T,N,ITPT,IT,OnCell}, d) = 1 <= d <= N ? sitp.ranges[d][end] + boundstep(sitp.ranges[d]) : throw(BoundsError())
 
+lbound{T,N,ITPT,IT}(sitp::ScaledInterpolation{T,N,ITPT,IT,OnGrid}, d, inds) =
+    sitp.ranges[d][1]
+lbound{T,N,ITPT,IT}(sitp::ScaledInterpolation{T,N,ITPT,IT,OnCell}, d, inds) =
+    sitp.ranges[d][1] - boundstep(sitp.ranges[d])
+ubound{T,N,ITPT,IT}(sitp::ScaledInterpolation{T,N,ITPT,IT,OnGrid}, d, inds) =
+    sitp.ranges[d][end]
+ubound{T,N,ITPT,IT}(sitp::ScaledInterpolation{T,N,ITPT,IT,OnCell}, d, inds) =
+    sitp.ranges[d][end] + boundstep(sitp.ranges[d])
+
 boundstep(r::StepRange) = r.step / 2
 boundstep(r::UnitRange) = 1//2