RFC: More extrapolation behaviors

[tomasaschan]

This PR re-implements most of the extrapolation behaviours that were defined before #31 merged.

Todo:

• Reflecting
• Periodic
• Linear
• make sure it works with DimSpecs
• Tests
• Squash
• Documentation

If you're missing something, let me know.

[timholy]

+1 on the variable rename

[tomasaschan]

I made a few too many typos for the history to be even remotely coherent, so I took the liberty to squash them away. I've fixed the looping - reflection is now done using modulo operations into a "double-domain", and if x is in the upper half it is reflected over the upper "single-domain" boundary, much like you suggested (and like Grid does).

The only other comment you had was a +1 on a parameter rename on a method that I've now been able to remove completely 😄

[timholy]

LGTM. I haven't really looked at the extrapolation stuff yet, but I will soon. I want fill (e.g., with NaN).

One question: do we clamp twice? Once in the itp and once in the etp?

[tomasaschan]

No, the itp doesn't do any bounds checking at all - it just uses the outermost polynomial as extrapolation. I've been thinking a little about what the most reasonable thing to do is, an I'm leaning towards throwing a bounds error unless in an @inbounds block, but I haven't had time to invest in trying something out.

[timholy]

If you grep for "clamp" you get 12 hits. For example, https://github.com/tlycken/Interpolations.jl/blob/6ce4b2a6cacd9bfefac5d8701a00fa157e849508/src/b-splines/constant.jl#L5

[tomasaschan]

Huh, yeah you're right. I'll have to take a look at that, but it will be at least September before I can do it.

[tomasaschan]

@timholy I think the only thing I want to do before merging this is adding some tests to make sure this doesn't break the DimSpec stuff. However, I'm a little unsure on exactly how to best test that (mainly because I haven't had time to look carefully at the implementation, nor seen any actual usage of them in other code). What's a good test case, and the expected behavior, for extrapolating DimSpecced interpolation objects?

[timholy]

Sorry I forgot about this question. Some examples are here.

[tomasaschan]

I've realized this approach might have a few limitations; since we instantiate different types of extrapolation objects for different extrapolation schemes, it's not possible to have different behaviors in different directions. That's kind of dull.

Unless anyone protests, I'm going to attempt a rewrite on the entire extrapolation implementation, that lets us do things like extrapolate(itp, Tuple{Flat, Linear}). In the long run, I would like to even support extrapolate(itp, Tuple{Flat, Tuple{Flat, Linear}}) to do constant extrapolation in the x-direction, and constant extrapolation for y < lbound(itp, 2) but linear extrapolation for y > ubound(itp, 2).

[timholy]

Sounds cool to me.

[tomasaschan]

I merged some new cool stuff, but periodic and reflecting extrapolation rely on mod, which isn't implemented for dual numbers yet. (I could turn the flag green by just removing those tests until that PR merges, but I'm afraid I'll forget to re-activate them...)

[tomasaschan]

@timholy, I'm really excited about this 💃

julia> using Interpolations

help?> extrapolate
search: extrapolate FilledExtrapolation AbstractExtrapolation

  extrapolate(itp, fillvalue) creates an extrapolation object that returns the fillvalue any
  time the indexes in itp[x1,x2,...] are out-of-bounds.

  extrapolate(itp, scheme) adds extrapolation behavior to an interpolation object, according
  to the provided scheme.

  The scheme can take any of these values:

    •  Throw - throws a BoundsError for out-of-bounds indices
    •  Flat - for constant extrapolation, taking the closest in-bounds value
    •  Linear - linear extrapolation (the wrapped interpolation object must support gradient)
    •  Reflect - reflecting extrapolation
    •  Periodic - periodic extrapolation

 You can also combine schemes in tuples. For example, the scheme Tuple{Linear, Flat} will
 use linear extrapolation in the first dimension, and constant in the second.

 Finally, you can specify different extrapolation behavior in different direction.
 Tuple{Tuple{Linear,Flat}, Flat} will extrapolate linearly in the first dimension if the index
  is too small, but use constant etrapolation if it is too large, and always use constant
  extrapolation in the second dimension.

[timholy]

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

[tomasaschan]

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.

[timholy]

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.

[timholy]

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.

[tomasaschan]

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.

[timholy]

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.

[tomasaschan]

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?

[timholy]

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.

[tomasaschan]

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?

[timholy]

👍 Very nice indeed!

[timholy]

I'm excited, too! I think this is the last step for matching Grid feature-for-feature? (And of course it exceeds Grid in so many ways.)

[tomasaschan]

The only thing lacking is cubic interpolation and evaluation of Hessians, but that never quite worked in Grid anyway, so yeah, I think all features that people actually used are there now (along with so much more, as you said).

[tomasaschan]

Woop woop! 🎉 🎈 🍻

[simonp0420]

As an interested observer and admirer of your work, please note that some of us make essential use of the Hessian calculations available in Grid. Hoping that eventually you add this functionality to Interpolations as well. In the meantime, thanks for all the great work you've put into this package.

[tomasaschan]

@simonp0420 Thanks for the kind words!

Cubic splines and Hessian evaluation is definitely on the roadmap; I won't be satisfied with this package until Grid.jl can be replaced entirely (even though the API is obviously quite different).