Break all the things

[timholy]

Here's an early preview of where I'm heading with my attempted a rewrite. Some of the goals are laid out in #224. Despite adding more tests and more docstrings, so far this is a tiny bit shorter than its predecessor, which may be a good sign.

A couple things I'm particularly happy about:

• I got rid of most of the generated functions. These were necessary when we wrote the package, but since then Julia tuples have gotten so fast you can do everything that way. In general things just work better with ordinary functions (e.g., world age issues, tracking with Revise, etc.)
• while there is still some complexity in the "generic" code, the customization that has to be done for a given interpolation mode (e.g., BSpline(Quadratic(BC()))) is considerably simpler and smaller than before. This may help attract user/developers.
• I'm using the axes of the returned object to indicate the region at which we "exactly" (up to numerical precision) reconstruct the input array. If you use interpolate, this always agrees with the parent, but with interpolate! the result may lose its edges depending on the interpolation order and choice of boundary conditions. I am pleased that we can now signal to the user where the values can plausibly be trusted; formerly that was only something gurus knew.
• This turns on bounds-checking for interpolate; use extrapolate if you want to evaluate out-of-bounds.
• Long ago (RFC: Use instances instead of types when specifying options #74) we switched our API from types to values, with the major advantage that, in principle, it allows one to store data related to the interpolation mode or boundary conditions. However, because internally most of those settings were tossed and all choices were represented only via the type system, there was no way we could exploit this possibility. Consequently, I've switched to storing the initialization settings as fields in the BSplineInterpolation object. The good news is that when they are all empty, there is no storage overhead:

julia> struct Empty end

julia> struct Wrapper1
           x::Float64
       end

julia> struct Wrapper2{E}
           x::Float64
           y::E
       end

julia> wx = Wrapper1(3.2)
Wrapper1(3.2)

julia> wy = Wrapper2(3.2, Empty())
Wrapper2{Empty}(3.2, Empty())

julia> sizeof(wx)
8

julia> sizeof(wy)
8

julia> struct NotEmpty
           z::UInt16
       end

julia> wz = Wrapper2(3.2, NotEmpty(2))
Wrapper2{NotEmpty}(3.2, NotEmpty(0x0002))

julia> sizeof(wz)
16

There are a couple of broken tests, and I've temporarily commented out Gridded, Scaled, and extrapolated objects. These will be addressed soon. But this should give people a chance to start looking it over and deciding whether they like it or not.

Left to do:

• add some hessian tests, and errors for anything that doesn't actually work
• reimplement extrapolation
• reimplement scale
• reimplement gridded

Related issues:

• fixes Confused about the true meaning of OnCell and OnGrid #228
• fixes Interaction with ForwardDiff.jl #213
• fixes Deprecation of array indexing calls [] #212
• fixes Gradients of 3D extrapolations not working. #202
• fixes Scaled interpolations failing on vector with grid points below 1.0 #200
• fixes Calling with function syntax? #198
• fixes Functions with Arrays of AbstractInterpolations cannot be handled by ForwardDiff #191
• fixes Julia 0.6 breaks scaling with LinSpace grid? #165
• fixes Scaled interpolation does not automatically extrapolate when called with a vector #158
• fixes Quadratic/Cubic Interpolations Not Working for Dimensions > 3 #151
• fixes Throw not exported? #127
• fixes Hessian fails for function of two variables #122

[tomasaschan]

I really like what you've done here, especially how you cut down the implementation of specific B-spline degrees to defining the weights (and boundary conditions). Amazing job, very well done!

[tomasaschan]

👍 👍 👍

[tomasaschan]

It makes the code easier to read if twotuple is defined before its usage here.

[tomasaschan]

I assume this is intended to show how to implement this for concrete implementations? Maybe we could name this evaluate(itp::AbstractInterpolation, x::Vararg{UnexpandedIndexTypes}) = itp(to_indices(itp, x)) to have an actual implementation that concrete interpolation types can use to fall back on?

[tomasaschan]

❗️

This is friggin' awesome. You are a genious. You just took the most complicated parts of this library to get right, and made them into "copy-from-wikipedia".

Have I told you how much I love the code you write?

[tomasaschan]

I think there might be bugs here (and we know from #122 and #181 that there are bugs in the Hessian coefficients as they stood before). These should just be the derivatives of the gradient coefficients, but as far as I can tell they aren't here: the third weight expression should be 3δx - 5.

As we don't have any tests for Hessians this is difficult to verify, so leave that change out of this PR if you wish. But it's really nice that the mathematics is so much easier to check here, because the code is much closer to the results of the derivation. Well done!

[timholy]

I think I'm going to add tests for Hessians and then insert deliberate throws into code-paths that don't work, just to prevent users from relying on stuff that's broken. The error message will include a plea for help 😄.

I've noticed many of our tests follow a general pattern, so I'm going to rework our tests so that there's a common core of test utilities. That should make it much easier to add the Hessian tests. Probably the very next priority.

[timholy]

I may also ultimately fix the Hessians myself, but once the rewrite is done I have higher priorities (other packages to upgrade to 1.0) since I don't use the hessian code myself.

[timholy]

Discussed at JuliaCon with @ChrisRackauckas (now is the time to complain about anything else!)

[timholy]

I expect this to pass tests, but don't merge: I've still got a lot commented out. Of note, with the new testing infrastructure now almost 90k tests run just for BSpline & NoInterp! But it's still quite fast, because most of the tests run inside compiled functions.

[pkofod]

Nobody's going to challenge your word, but do you have benchmarks to back up that it is faster? :)

[timholy]

Nobody's going to challenge your word, but do you have benchmarks to back up that it is faster? :)

You should! I actually don't yet know that it will be. But I did a couple of small tests before on simple problems (e.g., 1d linear interpolation) and was able to get better generated code than that produced by Interpolations. I don't think it will be a big difference.

But truly, my only goal is to make sure I don't (appreciably) slow anything down. If I get approximate performance parity and complete feature parity (fixing a whole bunch of bugs in the process), with more maintainable code, I'll consider it a win.

[pkofod]

But truly, my only goal is to make sure I don't (appreciably) slow anything down. If I get approximate performance parity and complete feature parity (fixing a whole bunch of bugs in the process), with more maintainable code, I'll consider it a win.

A big win, I agree! I have always had a hard time reading the internals of this package, and I really appreciate what you're doing here! I just saw in the associated issue that you said it'd be faster :)

[timholy]

OK, this is now extremely close: technically the only thing holding it up is an inference failure, which triggers failures in the new test/core.jl and in test/b-splines/constant.jl. To make sure this doesn't get merged prematurely here's the remaining list:

• fix or work around Inference failure dependent on Tuple signature annotation JuliaLang/julia#29117
• run the benchmarks and fix any serious regressions that are discovered

In the OP I've updated the list of bugs that this fixes, it's getting to be quite a long list!

[timholy]

At long last, this is ready to merge. Here are the benchmark results (old version execution time is along the horizontal axis, new version execution time is along the vertical axis, and the dashed line represents equal performance; points below the line indicate improved performance, and points above the line represent worsened performance):

You can see we're considerably faster on construction and basically on par or better in terms of run-time performance.

From a user perspective, the main thing to note is that this implements #228 (comment), meaning that now one uses

interpolate(A, BSpline(Quadratic(Line(OnCell()))))

rather than

interpolate(A, BSpline(Quadratic(Line())), OnCell())

The reason is that the old API (and accompanying docs) really weren't correct, see #228. This implementation makes it clearer that all OnGrid and OnCell do is change the interpretation of the boundary condition. Note that since Constant and Linear don't need a boundary condition, you should no longer supply OnGrid or OnCell for them.

Another really important behavioral change is that interpolation objects now check bounds and throw errors if you exceed them. If you want to evaluate beyond the bounds, you must use extrapolate. Note that you can turn off bounds-checking with @inbounds; the performance benchmarks above were collected with bounds checking on, and would be even more favorable if the accesses were made using @inbounds.

A final visible change is that I changed the extrapolation boundary condition from Linear to Line. Given that all the other options are consistent with the interpolation boundary conditions, I believe that was a typo that got enshrined; Linear is now reserved for the order of interpolation, and Line is used both for boundary conditions and extrapolation.

I believe I have added deprecation warnings for all of these changes, so nothing should break old code.

[lstagner]

I just tried out this branch. Great work. I have a question though. I noticed I can still use the getindex syntax itp[xval] if xval is an integer. Is this required for your implementation? I was hoping we could use that syntax as a way to access the parent array e.g. itp[i] = parent(itp)[i] that way the itp could act both as an array and as a interpolating function.

[timholy]

It's useful to have these still be AbstractArrays and therefore we have to support indexing with integers.

Keep in mind that the coefficients array may not match the input, if prefiltering was used (quadratic and cubic). Consequently itp[i] for integer i still performs the calculations of interpolation so as to "undo" the prefiltering. Is that what you're asking?

[lstagner]

Perhaps an example would better explain

x = linspace(0,2pi,100)
y = sin.(x)
itp_y = CubicSplineInterpolation(x,y)
itp_y(2.0) ≊ sin(2.0)
itp_y[50] === y[50] # I want this to be true

Basically I would like itp_y to act like y when using the getindex [] syntax with an index. Currently to get this behavior I have use parent(itp_y) to get the underlying y array.

[timholy]

Ah, yes, you're commenting on how it behaves when you have knots (that's an important detail, since BSpline works as-is). Basically you're asking for the separation in #227. That's not fully part of this, though the part you're asking for isn't that many additional lines of code. Doable, in other words, though perhaps worth a separate PR.

[floswald]

🎉 great!