Skip to content

Commit 6ea564b

Browse files
timholytimholy
committed
Speed up extrapolation by eliminating branches in indices/lbound/ubound
1 parent 0c656b3 commit 6ea564b

File tree

12 files changed

+72
-39
lines changed

12 files changed

+72
-39
lines changed

‎src/b-splines/b-splines.jl‎

Lines changed: 13 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -38,15 +38,24 @@ padding(itp::AbstractInterpolation) = padding(typeof(itp))
3838
padextract(pad::Integer, d) = pad
3939
padextract(pad::Tuple{Vararg{Integer}}, d) = pad[d]
4040

41-
lbound{T,N,TCoefs,IT}(itp::BSplineInterpolation{T,N,TCoefs,IT,OnGrid}, d) =
41+
lbound{T,N,TCoefs,IT}(itp::BSplineInterpolation{T,N,TCoefs,IT,OnGrid}, d::Integer) =
4242
first(indices(itp, d))
43-
ubound{T,N,TCoefs,IT}(itp::BSplineInterpolation{T,N,TCoefs,IT,OnGrid}, d) =
43+
ubound{T,N,TCoefs,IT}(itp::BSplineInterpolation{T,N,TCoefs,IT,OnGrid}, d::Integer) =
4444
last(indices(itp, d))
45-
lbound{T,N,TCoefs,IT}(itp::BSplineInterpolation{T,N,TCoefs,IT,OnCell}, d) =
45+
lbound{T,N,TCoefs,IT}(itp::BSplineInterpolation{T,N,TCoefs,IT,OnCell}, d::Integer) =
4646
first(indices(itp, d)) - 0.5
47-
ubound{T,N,TCoefs,IT}(itp::BSplineInterpolation{T,N,TCoefs,IT,OnCell}, d) =
47+
ubound{T,N,TCoefs,IT}(itp::BSplineInterpolation{T,N,TCoefs,IT,OnCell}, d::Integer) =
4848
last(indices(itp, d))+0.5
4949

50+
lbound{T,N,TCoefs,IT}(itp::BSplineInterpolation{T,N,TCoefs,IT,OnGrid}, d, inds) =
51+
first(inds)
52+
ubound{T,N,TCoefs,IT}(itp::BSplineInterpolation{T,N,TCoefs,IT,OnGrid}, d, inds) =
53+
last(inds)
54+
lbound{T,N,TCoefs,IT}(itp::BSplineInterpolation{T,N,TCoefs,IT,OnCell}, d, inds) =
55+
first(inds) - 0.5
56+
ubound{T,N,TCoefs,IT}(itp::BSplineInterpolation{T,N,TCoefs,IT,OnCell}, d, inds) =
57+
last(inds)+0.5
58+
5059
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)
5160

5261
function size{T,N,TCoefs,IT,GT,pad}(itp::BSplineInterpolation{T,N,TCoefs,IT,GT,pad}, d)

‎src/extrapolation/extrap_prep.jl‎

Lines changed: 7 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -38,12 +38,16 @@ coordinate transformations, but for anything else methods for the low- and high-
3838
value cases need to be implemented for each scheme.
3939
""" extrap_prep
4040

41-
extrap_prep{T}(::Type{T}, n::Val{1}) = extrap_prep(T, n, Val{1}())
41+
extrap_prep{T}(::Type{T}, n::Val{1}) = quote
42+
inds_etp = indices(etp)
43+
$(extrap_prep(T, n, Val{1}()))
44+
end
4245
extrap_prep{T}(::Type{Tuple{T}}, n::Val{1}) = extrap_prep(T, n)
4346
extrap_prep{T}(::Type{Tuple{T,T}}, n::Val{1}) = extrap_prep(T, n)
4447
extrap_prep{T}(::Type{Tuple{Tuple{T,T}}}, n::Val{1}) = extrap_prep(T, n)
4548
function extrap_prep{S,T}(::Type{Tuple{S,T}}, n::Val{1})
4649
quote
50+
inds_etp = indices(etp)
4751
$(extrap_prep(S, n, Val{1}(), Val{:lo}()))
4852
$(extrap_prep(T, n, Val{1}(), Val{:hi}()))
4953
end
@@ -54,15 +58,15 @@ extrap_prep{S,T}(::Type{Tuple{Tuple{S,T}}}, n::Val{1}) = extrap_prep(Tuple{S,T},
5458
extrap_prep{T<:Tuple}(::Type{T}, ::Val{1}) = :(throw(ArgumentError("The 1-dimensional extrap configuration $T is not supported")))
5559

5660
function extrap_prep{T,N}(::Type{T}, n::Val{N})
57-
exprs = Expr[]
61+
exprs = [:(inds_etp = indices(etp))]
5862
for d in 1:N
5963
push!(exprs, extrap_prep(T, n, Val{d}()))
6064
end
6165
return Expr(:block, exprs...)
6266
end
6367
function extrap_prep{N,T<:Tuple}(::Type{T}, n::Val{N})
6468
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)))")))
65-
exprs = Expr[]
69+
exprs = [:(inds_etp = indices(etp))]
6670
for d in 1:N
6771
Tdim = T.parameters[d]
6872
if Tdim <: Tuple

‎src/extrapolation/extrap_prep_gradient.jl‎

Lines changed: 8 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1,11 +1,15 @@
11
# See ?extrap_prep for documentation for all these methods
22

3-
extrap_prep{T}(g::Val{:gradient}, ::Type{T}, n::Val{1}) = extrap_prep(g, T, n, Val{1}())
3+
extrap_prep{T}(g::Val{:gradient}, ::Type{T}, n::Val{1}) = quote
4+
inds_etp = indices(etp)
5+
$(extrap_prep(g, T, n, Val{1}()))
6+
end
47
extrap_prep{T}(g::Val{:gradient}, ::Type{Tuple{T}}, n::Val{1}) = extrap_prep(g, T, n)
58
extrap_prep{T}(g::Val{:gradient}, ::Type{Tuple{T,T}}, n::Val{1}) = extrap_prep(g, T, n)
69
extrap_prep{T}(g::Val{:gradient}, ::Type{Tuple{Tuple{T,T}}}, n::Val{1}) = extrap_prep(g, T, n)
710
function extrap_prep{S,T}(g::Val{:gradient}, ::Type{Tuple{S,T}}, ::Val{1})
811
quote
12+
inds_etp = indices(etp)
913
$(extrap_prep(g, S, n, Val{1}(), Val{:lo}()))
1014
$(extrap_prep(g, T, n, Val{1}(), Val{:hi}()))
1115
end
@@ -16,12 +20,12 @@ extrap_prep{T<:Tuple}(::Val{:gradient}, ::Type{T}, ::Val{1}) = :(throw(ArgumentE
1620

1721

1822
function extrap_prep{T,N}(g::Val{:gradient}, ::Type{T}, n::Val{N})
19-
Expr(:block, [extrap_prep(g, T, n, Val{d}()) for d in 1:N]...)
23+
Expr(:block, :(inds_etp = indices(etp)), [extrap_prep(g, T, n, Val{d}()) for d in 1:N]...)
2024
end
2125

2226
function extrap_prep{T<:Tuple,N}(g::Val{:gradient}, ::Type{T}, n::Val{N})
2327
length(T.parameters) == N || return :(throw(ArgumentError("The $N-dimensional extrap configuration $T is not supported")))
24-
exprs = Expr[]
28+
exprs = [:(inds_etp = indices(etp))]
2529
for d in 1:N
2630
Tdim = T.parameters[d]
2731
if Tdim <: Tuple
@@ -47,4 +51,4 @@ function extrap_prep{S,T,N,d}(g::Val{:gradient}, ::Type{Tuple{S,T}}, n::Val{N},
4751
end
4852

4953
extrap_prep{T,N,d}(g::Val{:gradient}, ::Type{T}, n::Val{N}, dim::Val{d}) = extrap_prep(g, Tuple{T,T}, n, dim)
50-
extrap_prep(g::Val{:gradient}, args...) = extrap_prep(args...)
54+
extrap_prep(g::Val{:gradient}, args...) = extrap_prep(args...)

‎src/extrapolation/extrapolation.jl‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -81,8 +81,10 @@ end
8181

8282
lbound(etp::Extrapolation, d) = lbound(etp.itp, d)
8383
ubound(etp::Extrapolation, d) = ubound(etp.itp, d)
84+
lbound(etp::Extrapolation, d, inds) = lbound(etp.itp, d, inds)
85+
ubound(etp::Extrapolation, d, inds) = ubound(etp.itp, d, inds)
8486
size(etp::Extrapolation, d) = size(etp.itp, d)
85-
indices(etp::AbstractExtrapolation) = indices(etp.itp)
87+
@inline indices(etp::AbstractExtrapolation) = indices(etp.itp)
8688
indices(etp::AbstractExtrapolation, d) = indices(etp.itp, d)
8789

8890
include("filled.jl")

‎src/extrapolation/filled.jl‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -26,7 +26,8 @@ function getindex_impl{T,N,ITP,IT,GT,FT}(fitp::Type{FilledExtrapolation{T,N,ITP,
2626
$meta
2727
# Check to see if we're in the extrapolation region, i.e.,
2828
# out-of-bounds in an index
29-
@nexprs $N d->((args[d] < lbound(fitp,d) || args[d] > ubound(fitp, d))) && return convert($Tret, fitp.fillvalue)::$Tret
29+
inds_etp = indices(fitp)
30+
@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
3031
# In the interpolation region
3132
return convert($Tret, getindex(fitp.itp,args...))::$Tret
3233
end
@@ -41,3 +42,5 @@ getindex{T}(fitp::FilledExtrapolation{T,1}, x::Number, y::Int) = y == 1 ? fitp[x
4142

4243
lbound(etp::FilledExtrapolation, d) = lbound(etp.itp, d)
4344
ubound(etp::FilledExtrapolation, d) = ubound(etp.itp, d)
45+
lbound(etp::FilledExtrapolation, d, inds) = lbound(etp.itp, d, inds)
46+
ubound(etp::FilledExtrapolation, d, inds) = ubound(etp.itp, d, inds)

‎src/extrapolation/flat.jl‎

Lines changed: 8 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -1,25 +1,25 @@
11
function extrap_prep{N,d}(::Type{Flat}, ::Val{N}, ::Val{d})
22
xs_d = Symbol("xs_", d)
3-
:($xs_d = clamp($xs_d, lbound(etp, $d), ubound(etp, $d)))
3+
:($xs_d = clamp($xs_d, lbound(etp, $d, inds_etp[$d]), ubound(etp, $d, inds_etp[$d])))
44
end
55

66
function extrap_prep{N,d}(::Type{Flat}, ::Val{N}, ::Val{d}, ::Val{:lo})
77
xs_d = Symbol("xs_", d)
8-
:($xs_d = max($xs_d, lbound(etp, $d)))
8+
:($xs_d = max($xs_d, lbound(etp, $d, inds_etp[$d])))
99
end
1010

1111
function extrap_prep{N,d}(::Type{Flat}, ::Val{N}, ::Val{d}, ::Val{:hi})
1212
xs_d = Symbol("xs_", d)
13-
:($xs_d = min($xs_d, ubound(etp, $d)))
13+
:($xs_d = min($xs_d, ubound(etp, $d, inds_etp[$d])))
1414
end
1515

1616
function extrap_prep{N,d}(::Val{:gradient}, ::Type{Flat}, ::Val{N}, ::Val{d}, ::Val{:lo})
1717
coords = [Symbol("xs_", k) for k in 1:N]
1818
xs_d = coords[d]
1919

2020
quote
21-
if $xs_d < lbound(etp, $d)
22-
$xs_d = lbound(etp, $d)
21+
if $xs_d < lbound(etp, $d, inds_etp[$d])
22+
$xs_d = lbound(etp, $d, inds_etp[$d])
2323
gradient!(g, etp.itp, $(coords...))
2424
g[$d] = 0
2525
return g
@@ -32,11 +32,11 @@ function extrap_prep{N,d}(::Val{:gradient}, ::Type{Flat}, ::Val{N}, ::Val{d}, ::
3232
xs_d = coords[d]
3333

3434
quote
35-
if $xs_d > ubound(etp, $d)
36-
$xs_d = ubound(etp, $d)
35+
if $xs_d > ubound(etp, $d, inds_etp[$d])
36+
$xs_d = ubound(etp, $d, inds_etp[$d])
3737
gradient!(g, etp.itp, $(coords...))
3838
g[$d] = 0
3939
return g
4040
end
4141
end
42-
end
42+
end

‎src/extrapolation/linear.jl‎

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -2,8 +2,8 @@ function extrap_prep{N,d}(::Type{Linear}, ::Val{N}, ::Val{d}, ::Val{:lo})
22
coords = [Symbol("xs_", k) for k in 1:N]
33
xs_d = coords[d]
44
quote
5-
if $xs_d < lbound(etp.itp, $d)
6-
$xs_d = lbound(etp.itp, $d)
5+
if $xs_d < lbound(etp.itp, $d, inds_etp[$d])
6+
$xs_d = lbound(etp.itp, $d, inds_etp[$d])
77
return etp[$(coords...)] + gradient(etp, $(coords...))[$d] * (xs[$d] - $xs_d)
88
end
99
end
@@ -12,8 +12,8 @@ function extrap_prep{N,d}(::Type{Linear}, ::Val{N}, ::Val{d}, ::Val{:hi})
1212
coords = [Symbol("xs_", k) for k in 1:N]
1313
xs_d = coords[d]
1414
quote
15-
if $xs_d > ubound(etp, $d)
16-
$xs_d = ubound(etp, $d)
15+
if $xs_d > ubound(etp, $d, inds_etp[$d])
16+
$xs_d = ubound(etp, $d, inds_etp[$d])
1717
return etp[$(coords...)] + gradient(etp, $(coords...))[$d] * (xs[$d] - $xs_d)
1818
end
1919
end
@@ -23,4 +23,4 @@ extrap_prep{N,d,l}(::Val{:gradient}, ::Type{Linear}, n::Val{N}, dim::Val{d}, loh
2323
extrap_prep(Flat, n, dim, lohi)
2424

2525
extrap_prep{N,d}(::Val{:gradient}, ::Type{Linear}, n::Val{N}, dim::Val{d}) =
26-
extrap_prep(Flat, n, dim)
26+
extrap_prep(Flat, n, dim)

‎src/extrapolation/periodic.jl‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -5,23 +5,23 @@ Translate x into the domain [lbound, ubound] my means of `mod()`
55
"""
66
function extrap_prep_dim(::Type{Periodic}, d)
77
xs_d = Symbol("xs_", d)
8-
:($xs_d = mod(xs[$d] - lbound(etp.itp, $d), ubound(etp.itp, $d) - lbound(etp.itp, $d)) + lbound(etp.itp, $d))
8+
:($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]))
99
end
1010

1111
extrap_prep{N,d}(::Type{Periodic}, ::Val{N}, ::Val{d}) = extrap_prep_dim(Periodic, d)
1212
function extrap_prep{N,d}(::Type{Periodic}, ::Val{N}, ::Val{d}, ::Val{:lo})
1313
xs_d = Symbol("xs_", d)
1414
quote
15-
if $xs_d < lbound(etp.itp, $d)
15+
if $xs_d < lbound(etp.itp, $d, inds_etp[$d])
1616
$(extrap_prep_dim(Periodic, d))
1717
end
1818
end
1919
end
2020
function extrap_prep{N,d}(::Type{Periodic}, ::Val{N}, ::Val{d}, ::Val{:hi})
2121
xs_d = Symbol("xs_", d)
2222
quote
23-
if $xs_d > ubound(etp.itp, d)
23+
if $xs_d > ubound(etp.itp, d, inds_etp[$d])
2424
$(extrap_prep_dim(Periodic, d))
2525
end
2626
end
27-
end
27+
end

‎src/extrapolation/reflect.jl‎

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -7,8 +7,8 @@ Next, if x is now in the upper part of this ''double-domain´´, reflect over th
77
function extrap_prep_dim(::Type{Reflect}, d)
88
xs_d = Symbol("xs_", d)
99
quote
10-
start = lbound(etp.itp, $d)
11-
width = ubound(etp.itp, $d) - start
10+
start = lbound(etp.itp, $d, inds_etp[$d])
11+
width = ubound(etp.itp, $d, inds_etp[$d]) - start
1212

1313
$xs_d = mod($xs_d - start, 2width) + start
1414
$xs_d > start + width && (xs_d = start + width - $xs_d)
@@ -18,9 +18,9 @@ end
1818
extrap_prep{N,d}(::Type{Reflect}, ::Val{N}, ::Val{d}) = extrap_prep_dim(Reflect, d)
1919
function extrap_prep{N,d}(::Type{Reflect}, ::Val{N}, ::Val{d}, ::Val{:lo})
2020
xs_d = Symbol("xs_", d)
21-
:($xs_d < lbound(etp.itp, $d) && $(extrap_prep_dim(Reflect, d)))
21+
:($xs_d < lbound(etp.itp, $d, inds_etp[$d]) && $(extrap_prep_dim(Reflect, d)))
2222
end
2323
function extrap_prep{N,d}(::Type{Reflect}, ::Val{N}, ::Val{d}, ::Val{:hi})
2424
xs_d = Symbol("xs_", d)
25-
:($xs_d > ubound(etp.itp, $d) && $(extrap_prep_dim(Reflect, d)))
26-
end
25+
:($xs_d > ubound(etp.itp, $d, inds_etp[$d]) && $(extrap_prep_dim(Reflect, d)))
26+
end

‎src/extrapolation/throw.jl‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1,14 +1,14 @@
11
function extrap_prep{N,d}(::Type{Throw}, ::Val{N}, ::Val{d})
22
xsym = Symbol("xs_", d)
3-
:(lbound(etp, $d) <= $xsym <= ubound(etp, $d) || throw(BoundsError()))
3+
:(lbound(etp, $d, inds_etp[$d]) <= $xsym <= ubound(etp, $d, inds_etp[$d]) || throw(BoundsError()))
44
end
55

66
function extrap_prep{N,d}(::Type{Throw}, ::Val{N}, ::Val{d}, ::Val{:lo})
77
xsym = Symbol("xs_", d)
8-
:(lbound(etp, $d) <= $xsym || throw(BoundsError()))
8+
:(lbound(etp, $d, inds_etp[$d]) <= $xsym || throw(BoundsError()))
99
end
1010

1111
function extrap_prep{N,d}(::Type{Throw}, ::Val{N}, ::Val{d}, ::Val{:hi})
1212
xsym = Symbol("xs_", d)
13-
:($xsym <= ubound(etp, $d) || throw(BoundsError()))
14-
end
13+
:($xsym <= ubound(etp, $d, inds_etp[$d]) || throw(BoundsError()))
14+
end

0 commit comments

Comments
 (0)