Skip to content

RFC: More extrapolation behaviors #50

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
tomasaschan merged 3 commits into from
Sep 29, 2015
Merged

RFC: More extrapolation behaviors #50

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
2197549
Refactor extrapolation, and add behaviors Periodic and Reflect
tomasaschan Jul 22, 2015
e37dc88
Fix splatting-related perf problem
tomasaschan Sep 29, 2015
e763408
Improve docs for extrap_prep
tomasaschan Sep 29, 2015
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 6 additions & 1 deletion src/Interpolations.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -67,8 +67,13 @@ lbound{T,N}(itp::AbstractInterpolation{T,N}, d) = convert(T, 1)
ubound{T,N}(itp::AbstractInterpolation{T,N}, d) = convert(T, size(itp, d))
itptype{T,N,IT,GT}(itp::AbstractInterpolation{T,N,IT,GT}) = IT
gridtype{T,N,IT,GT}(itp::AbstractInterpolation{T,N,IT,GT}) = GT
count_interp_dims{T<:AbstractInterpolation}(::Type{T}, N) = N

@inline gradient{T,N}(itp::AbstractInterpolation{T,N}, xs...) = gradient!(Array(T,N), itp, xs...)
@generated function gradient{T,N}(itp::AbstractInterpolation{T,N}, xs...)
n = count_interp_dims(itp, N)
Tg = promote_type(T, [x <: AbstractArray ? eltype(x) : x for x in xs]...)
:(gradient!(Array($Tg,$n), itp, xs...))
Copy link
Member

@timholy timholy Sep 29, 2015

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Might need to insert a $(Expr(:meta, :inline)) here if you don't want to pay a splatting penalty. (Worth testing.)

Copy link
Contributor Author

@tomasaschan tomasaschan Sep 29, 2015

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Definitely worth testing.

I'm not sure if I did this correctly, so happy for any comments on the following benchmark (mostly if this would be the correct way to implement the inlining):

function bar(A, xs...)
    A
end

@generated function foo1(xs...)
    n = length(xs)
    T = promote_type(xs...)
    :(bar(Array($T,$n), xs...))
end

@generated function foo2(xs...)
    n=length(xs)
    T = promote_type(xs...)
    quote
        $(Expr(:meta, :inline))
        bar(Array($T,$n), xs...)
    end
end

# after warmup

julia> gc(); @time for _ in 1:1e5; foo1(1,2,3); end
  0.004200 seconds (100.00 k allocations: 9.155 MB)

julia> gc(); @time for _ in 1:1e5; foo2(1,2,3); end
  0.003687 seconds (100.00 k allocations: 9.155 MB)

Repeating the final two incantations give results somewhere in the range 0.003 to 0.005 seconds for both functions (I'd have to start collecting larger samples and look at the distributions of timing results to see a difference). If my benchmarking strategy is valid, I think it's safe to assume that splatting will be negligible compared to array allocation.

Copy link
Member

@timholy timholy Sep 29, 2015

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For a test this simple, bar is being inlined automatically:

julia> @code_typed foo1(1,2,3)
1-element Array{Any,1}:
 :($(Expr(:lambda, Any[:(xs::Any...)], Any[Any[Any[:xs,Tuple{Int64,Int64,Int64},0],Any[symbol("##xs#6818"),Tuple{Int64,Int64,Int64},0]],Any[],Any[Int64,Array{Int64,1}],Any[]], :(begin  # none, line 2:
        GenSym(1) = (top(ccall))(:jl_alloc_array_1d,(top(apply_type))(Base.Array,Int64,1)::Type{Array{Int64,1}},(top(svec))(Base.Any,Base.Int)::SimpleVector,Array{Int64,1},0,3,0)::Array{Int64,1}
        return GenSym(1)
    end::Array{Int64,1}))))

You can see the absence of a top(call)(:bar, ... expression. You can fix that by adding @noinline in front of the definition of bar.

That said, there are several other non-ideal aspects of this test, including (1) you're not generating a version that elides the splat in the call to bar, and (2) you're doing something really expensive in your test which can mask the impact of spaltting. Here's a better one:

@noinline function bar1(A, xs...)
    A[xs...]
end

@inline function bar2(A, xs...)
    A[xs...]
end

function call_bar1a(A, n, xs...)
    s = zero(eltype(A))
    for i = 1:n
        s += bar1(A, xs...)
    end
    s
end

function call_bar2a(A, n, xs...)
    s = zero(eltype(A))
    for i = 1:n
        s += bar2(A, xs...)
    end
    s
end

@generated function call_bar1b(A, n, xs...)
    xargs = [:(xs[$d]) for d = 1:length(xs)]
    quote
        s = zero(eltype(A))
        for i = 1:n
            s += bar1(A, $(xargs...))
        end
        s
    end
end

@generated function call_bar2b(A, n, xs...)
    xargs = [:(xs[$d]) for d = 1:length(xs)]
    quote
        s = zero(eltype(A))
        for i = 1:n
            s += bar2(A, $(xargs...))
        end
        s
    end
end

A = rand(3,3,3)
call_bar1a(A, 1, 1, 2, 3)
call_bar2a(A, 1, 1, 2, 3)
call_bar1b(A, 1, 1, 2, 3)
call_bar2b(A, 1, 1, 2, 3)

@time 1
@time call_bar1a(A, 10^6, 1, 2, 3)
@time call_bar2a(A, 10^6, 1, 2, 3)
@time call_bar1b(A, 10^6, 1, 2, 3)
@time call_bar2b(A, 10^6, 1, 2, 3)

Results:

julia> include("/tmp/test_splat.jl")
  0.000001 seconds (3 allocations: 144 bytes)
  0.760468 seconds (7.00 M allocations: 137.329 MB, 3.72% gc time)
  0.761979 seconds (7.00 M allocations: 137.329 MB, 3.62% gc time)
  0.563945 seconds (5.00 M allocations: 106.812 MB, 4.30% gc time)
  0.003080 seconds (6 allocations: 192 bytes)

As you can see, the difference is not small 😄. I was half-expecting call_bar2a to be fast, so even call-site splatting is deadly.

In this case, eliding the splatting in call_bar* is totally unimportant because it's only called once; the analog of what I was suggesting with the :meta expression are the @noinline/@inline statements in front of bar1 and bar2.

Copy link
Member

@timholy timholy Sep 29, 2015

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

As a further general point, avoid deliberate allocations when you're benchmarking stuff. That way any allocations that happen because of type-instability, splatting, etc, show up very clearly.

Copy link
Contributor Author

@tomasaschan tomasaschan Sep 29, 2015

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah, thanks a lot! As always, you don't only help me write the code I need, you also help me understand why I need it :)

The main reason I didn't write a benchmark without the array allocation is that in the actual use case has it, so I wanted to see which effect was important. But the trick with xargs seems much cleaner anyway, so I'll go with that and be happy.

end

include("nointerp/nointerp.jl")
include("b-splines/b-splines.jl")
Expand Down
19 changes: 2 additions & 17 deletions src/b-splines/b-splines.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -41,6 +41,8 @@ ubound{T,N,TCoefs,IT}(itp::BSplineInterpolation{T,N,TCoefs,IT,OnGrid}, d) = conv
lbound{T,N,TCoefs,IT}(itp::BSplineInterpolation{T,N,TCoefs,IT,OnCell}, d) = convert(T, .5)
ubound{T,N,TCoefs,IT}(itp::BSplineInterpolation{T,N,TCoefs,IT,OnCell}, d) = convert(T, size(itp, d) + .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)

@generated function size{T,N,TCoefs,IT,GT,pad}(itp::BSplineInterpolation{T,N,TCoefs,IT,GT,pad}, d)
quote
d <= $N ? size(itp.coefs, d) - 2*padextract($pad, d) : 1
Expand Down Expand Up @@ -71,23 +73,6 @@ function interpolate!{IT<:DimSpec{BSpline},GT<:DimSpec{GridType}}(A::AbstractArr
interpolate!(tweight(A), A, IT, GT)
end

define_indices{IT}(::Type{IT}, N, pad) = Expr(:block, Expr[define_indices_d(iextract(IT, d), d, padextract(pad, d)) for d = 1:N]...)

coefficients{IT}(::Type{IT}, N) = Expr(:block, Expr[coefficients(iextract(IT, d), N, d) for d = 1:N]...)

function gradient_coefficients{IT<:DimSpec{BSpline}}(::Type{IT}, N, dim)
exs = Expr[d==dim ? gradient_coefficients(iextract(IT, dim), d) :
coefficients(iextract(IT, d), N, d) for d = 1:N]
Expr(:block, exs...)
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")
Expand Down
13 changes: 12 additions & 1 deletion src/b-splines/indexing.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,18 @@ import Base.getindex

Base.linearindexing{T<:AbstractInterpolation}(::Type{T}) = Base.LinearSlow()

define_indices{IT}(::Type{IT}, N, pad) = Expr(:block, Expr[define_indices_d(iextract(IT, d), d, padextract(pad, d)) for d = 1:N]...)

coefficients{IT}(::Type{IT}, N) = Expr(:block, Expr[coefficients(iextract(IT, d), N, d) for d = 1:N]...)

function gradient_coefficients{IT<:DimSpec{BSpline}}(::Type{IT}, N, dim)
exs = Expr[d==dim ? gradient_coefficients(iextract(IT, dim), d) :
coefficients(iextract(IT, d), N, d) for d = 1:N]
Expr(:block, exs...)
end

index_gen{IT}(::Type{IT}, N::Integer, offsets...) = index_gen(iextract(IT, min(length(offsets)+1, N)), IT, N, offsets...)

function getindex_impl{T,N,TCoefs,IT<:DimSpec{BSpline},GT<:DimSpec{GridType},Pad}(itp::Type{BSplineInterpolation{T,N,TCoefs,IT,GT,Pad}})
meta = Expr(:meta, :inline)
quote
Expand Down Expand Up @@ -57,7 +69,6 @@ function gradient_impl{T,N,TCoefs,IT<:DimSpec{BSpline},GT<:DimSpec{GridType},Pad
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)
Expand Down
18 changes: 0 additions & 18 deletions src/extrapolation/constant.jl
View file
Open in desktop

This file was deleted.

15 changes: 0 additions & 15 deletions src/extrapolation/error.jl
View file
Open in desktop

This file was deleted.

138 changes: 138 additions & 0 deletions src/extrapolation/extrap_prep.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,138 @@
"""
Copy link
Member

@timholy timholy Sep 29, 2015

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The commentsdocs are helpful, but it might be even easier to understand with some kind of overarching comment at the top. Otherwise these help phrases don't mean a lot.

Copy link
Contributor Author

@tomasaschan tomasaschan Sep 29, 2015

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I had an implementation that didn't use dispatch so much but rather had a bunch of (read: ton of) compile-time if-clauses. Would that be easier to understand, do you think?

Copy link
Member

@timholy timholy Sep 29, 2015

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I could see either one.

The bigger issue is this: presumably these are non-exported functions. I presume you've written help for them as encouragement for others to contribute to the development of Interpolations, which I definitely support. However, these are useful only if there's a little bit of context provided somewhere to aid in understanding the overarching purpose of these methods. In other words, "1-dimensional, same lo/hi schemes" (as the first lines of text in the file) is nearly useless on its own.

Copy link
Contributor Author

@tomasaschan tomasaschan Sep 29, 2015

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I now removed all the one-liners that explained each method separately (I figure that understanding can be reconstructed with @which, if one really wants it) and added a more general overview of how these methods are constructed to work together, and what one needs to implement for new schemes. Is this understandable for someone who didn't write the code?

Copy link
Member

@timholy timholy Sep 29, 2015

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

👍 Very nice indeed!

`extrap_prep(T, Val{1}())`

1-dimensional, same lo/hi schemes
"""
extrap_prep{T}(::Type{T}, n::Val{1}) = extrap_prep(T, n, Val{1}())

"""
`extrap_prep(Tuple{T}, Val{1}())`

1-dimensional, same lo/hi schemes but specified as if N-dimensional

Equivalent with `extrap_prep(T, Val{1}())`
"""
extrap_prep{T}(::Type{Tuple{T}}, n::Val{1}) = extrap_prep(T, n)

"""
`extrap_prep(Tuple{T,T}, Val{1}())`

1-dimensional, same lo/hi schemes but specified as tuple

Equivalent with `extrap_prep(T, Val{1}())`
"""
extrap_prep{T}(::Type{Tuple{T,T}}, n::Val{1}) = extrap_prep(T, n)

"""
`extrap_prep(Tuple{Tuple{T,T}}, Val{1}())`

1-dimensional, same lo/hi schemes but specified as tuple and as if N-dimensional

Equivalent with `extrap_prep(T, Val{1}())`
"""
extrap_prep{T}(::Type{Tuple{Tuple{T,T}}}, n::Val{1}) = extrap_prep(T, n)

"""
`extrap_prep(Tuple{S,T}, Val{1}())`

1-dimensional, different lo/hi schemes
"""
function extrap_prep{S,T}(::Type{Tuple{S,T}}, n::Val{1})
quote
$(extrap_prep(S, n, Val{1}(), Val{:lo}()))
$(extrap_prep(T, n, Val{1}(), Val{:hi}()))
end
end

"""
`extrap_prep(Tuple{Tuple{S,T}}, Val{1}())`

1-dimensional, different lo/hi schemes but specified as if N-dimensional

Equivalent with `extrap_prep(Tuple{S,T}, Val{1}())`
"""
extrap_prep{S,T}(::Type{Tuple{Tuple{S,T}}}, n::Val{1}) = extrap_prep(Tuple{S,T}, n)

"""
`extrap_prep(Tuple{...}, Val{1}())`

1-dimensional, but incorrect tuple spec
""" # needed for ambiguity resolution
extrap_prep{T<:Tuple}(::Type{T}, ::Val{1}) = :(throw(ArgumentError("The 1-dimensional extrap configuration $T is not supported")))

"""
`extrap_prep(T, Val{N}())`

N-dimensional, all schemes same
"""
function extrap_prep{T,N}(::Type{T}, n::Val{N})
exprs = Expr[]
for d in 1:N
push!(exprs, extrap_prep(T, n, Val{d}()))
end
return Expr(:block, exprs...)
end

"""
`extrap_prep(Tuple{...}, Val{N}())`

N-dimensional, different specs in different dimensions.

Note that the tuple must be of length N.
"""
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 $(lenght(T.parameters)))")))
exprs = Expr[]
for d in 1:N
Tdim = T.parameters[d]
if Tdim <: Tuple
length(Tdim.parameters) == 2 || return :(throw(ArgumentError("The extrap configuration $Tdim for dimension $d is not supported - must be a tuple of length 2 or a simple configuration type")))
if Tdim.parameters[1] != Tdim.parameters[2]
push!(exprs, extrap_prep(Tdim, n, Val{d}()))
else
push!(exprs, extrap_prep(Tdim.parameters[1], n, Val{d}()))
end
else
push!(exprs, extrap_prep(Tdim, n, Val{d}()))
end
end
return Expr(:block, exprs...)
end

"""
`extrap_prep(T, Val{N}(), Val{d}())`

N-dimensional fallback for expanding the same scheme lo/hi in a single dimension
"""
function extrap_prep{T,N,d}(::Type{T}, n::Val{N}, dim::Val{d})
quote
$(extrap_prep(T, n, dim, Val{:lo}()))
$(extrap_prep(T, n, dim, Val{:hi}()))
end
end

"""
`extrap_prep(Tuple{S,T}, Val{N}(), Val{d}())`

N-dimensional fallback for expanding the different schemes lo/hi in a single dimension
"""
function extrap_prep{S,T,N,d}(::Type{Tuple{S,T}}, n::Val{N}, dim::Val{d})
quote
$(extrap_prep(S, n, dim, Val{:lo}()))
$(extrap_prep(T, n, dim, Val{:hi}()))
end
end

"""
`extrap_prep(Tuple{T,T}, Val{N}(), Val{d}())`

N-dimensional fallback for expanding the same schemes lo/hi in a single dimension

Equivalent with `extrap_prep(T, Val{N}(), Val{d}())`
"""
function extrap_prep{T,N,d}(::Type{Tuple{T,T}}, n::Val{N}, dim::Val{d})
quote
$(extrap_prep(T, n, dim, Val{:lo}()))
$(extrap_prep(T, n, dim, Val{:hi}()))
end
end
115 changes: 115 additions & 0 deletions src/extrapolation/extrap_prep_gradient.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,115 @@
"""
`extrap_prep(Val{:gradient}(), T, Val{1}())`

1-dimensional for gradient, same lo/hi schemes
"""
extrap_prep{T}(g::Val{:gradient}, ::Type{T}, n::Val{1}) = extrap_prep(g, T, n, Val{1}())

"""
`extrap_prep(Val{:gradient}(), Tuple{T}, Val{1}())`

1-dimensional for gradient, same lo/hi schemes but specified as if N-dimensional

Equivalent with `extrap_prep(Val{:gradient}(), T, Val{1}())`
"""
extrap_prep{T}(g::Val{:gradient}, ::Type{Tuple{T}}, n::Val{1}) = extrap_prep(g, T, n)

"""
`extrap_prep(Val{:gradient}(), Tuple{T,T}, Val{1}())`

1-dimensional for gradient, same lo/hi schemes, but specified as tuple

Equivalent with `extrap_prep(Val{:gradient}(), T, Val{1}())`
"""
extrap_prep{T}(g::Val{:gradient}, ::Type{Tuple{T,T}}, n::Val{1}) = extrap_prep(g, T, n)

"""
`extrap_prep(Val{:gradient}(), Tuple{Tuple{T}}, Val{1}())`

1-dimensional for gradient, same lo/hi schemes, but specified as tuple and N-dimensional

Equivalent with `extrap_prep(Val{:gradient}(), T, Val{1}())`
"""
extrap_prep{T}(g::Val{:gradient}, ::Type{Tuple{Tuple{T,T}}}, n::Val{1}) = extrap_prep(g, T, n)


"""
`extrap_prep(Val{:gradient}(), Tuple{S,T}, Val{1}())`

1-dimensional for gradient, different lo/hi schemes
"""
function extrap_prep{S,T}(g::Val{:gradient}, ::Type{Tuple{S,T}}, ::Val{1})
quote
$(extrap_prep(g, S, n, Val{1}(), Val{:lo}()))
$(extrap_prep(g, T, n, Val{1}(), Val{:hi}()))
end
end

"""
`extrap_prep(Val{:gradient}(), Tuple{Tuple{S,T}}, Val{1}())`

1-dimensional for gradient, different lo/hi schemes but specified as if N-dimensional

Equivalent with `extrap_prep(Val{:gradient}(), Tuple{S,T}, Val{1}())`
"""
extrap_prep{S,T}(g::Val{:gradient}, ::Type{Tuple{Tuple{S,T}}}, n::Val{1}) = extrap_prep(g, Tuple{S,T}, n)

"""
`extrap_prep(Val{:gradient}(), Tuple{...}, Val{1}())`

1-dimensional for gradient, but incorrect tuple spec
""" # needed for ambiguity resolution
extrap_prep{T<:Tuple}(::Val{:gradient}, ::Type{T}, ::Val{1}) = :(throw(ArgumentError("The 1-dimensional extrap configuration $T is not supported")))


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]...)
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[]
for d in 1:N
Tdim = T.parameters[d]
if Tdim <: Tuple
length(Tdim.parameters) == 2 || return :(throw(ArgumentError("The extrap configuration $Tdim for dimension $d is not supported - must be a tuple of length 2 or a simple configuration type"))
)
if Tdim.parameters[1] != Tdim.parameters[2]
push!(exprs, extrap_prep(g, Tdim, n, Val{d}()))
else
push!(exprs, extrap_prep(g, Tdim.parameters[1], n, Val{d}()))
end
else
push!(exprs, extrap_prep(g, Tdim, n, Val{d}()))
end
end
return Expr(:block, exprs...)
end

function extrap_prep{S,T,N,d}(g::Val{:gradient}, ::Type{Tuple{S,T}}, n::Val{N}, dim::Val{d})
quote
$(extrap_prep(g, S, n, dim, Val{:lo}()))
$(extrap_prep(g, T, n, dim, Val{:hi}()))
end
end

function extrap_prep{T,N,d}(g::Val{:gradient}, ::Type{Tuple{T,T}}, n::Val{N}, dim::Val{d})
quote
$(extrap_prep(g, T, n, dim, Val{:lo}()))
$(extrap_prep(g, T, n, dim, Val{:hi}()))
end
end

function extrap_prep{T,N,d}(g::Val{:gradient}, ::Type{T}, n::Val{N}, dim::Val{d})
quote
$(extrap_prep(g, T, n, dim, Val{:lo}()))
$(extrap_prep(g, T, n, dim, Val{:hi}()))
end
end

"""
`extrap_prep(Val{:gradient}(), args...)`

Fallback that defaults to the same behavior as for value extrapolation (works e.g. for all schemes that are coordinate transformations).
"""
extrap_prep(g::Val{:gradient}, args...) = extrap_prep(args...)
Loading