-
Notifications
You must be signed in to change notification settings - Fork 109
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.
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
Break all the things #226
Break all the things #226
Conversation
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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!
src/Interpolations.jl
Outdated
- `InPlace()` | ||
- `InPlaceQ()` | ||
""" | ||
abstract type BoundaryCondition <: Flag end |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
👍 👍 👍
src/Interpolations.jl
Outdated
axes(itp::AbstractInterpolation) = itp.parentaxes | ||
axes(exp::AbstractExtrapolation) = axes(exp.itp) | ||
|
||
bounds(itp::AbstractInterpolation) = map(twotuple, lbounds(itp), ubounds(itp)) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
It makes the code easier to read if twotuple
is defined before its usage here.
src/Interpolations.jl
Outdated
# Commented out because you can't add methods to an abstract type. | ||
# @inline function (itp::AbstractInterpolation)(x::Vararg{UnexpandedIndexTypes}) | ||
# itp(to_indices(itp, x)) | ||
# end |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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?
return :(itp.coefs[$(indices...)]) | ||
end | ||
end | ||
padded_axis(ax::AbstractUnitRange, ::BSpline{Constant}) = ax |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
❗️
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?
src/b-splines/cubic.jl
Outdated
c_d = q''(fx_d) | ||
cp_d = q''(1-fx_d) | ||
cpp_d = p''(1-fx_d) | ||
hessian_weights(::Cubic, δx) = (1-δx, 3*δx-2, 3*(1-δx)-2, δx) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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!
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think I'm going to add tests for Hessians and then insert deliberate throw
s 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.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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.
Discussed at JuliaCon with @ChrisRackauckas (now is the time to complain about anything else!) |
f862ac9
to
cc4c6de
Compare
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. |
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. |
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 :) |
OK, this is now extremely close: technically the only thing holding it up is an inference failure, which triggers failures in the new
In the OP I've updated the list of bugs that this fixes, it's getting to be quite a long list! |
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 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 A final visible change is that I changed the extrapolation boundary condition from I believe I have added deprecation warnings for all of these changes, so nothing should break old code. |
I just tried out this branch. Great work. I have a question though. I noticed I can still use the getindex syntax |
It's useful to have these still be Keep in mind that the coefficients array may not match the input, if prefiltering was used (quadratic and cubic). Consequently |
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 |
Ah, yes, you're commenting on how it behaves when you have knots (that's an important detail, since |
Storing parentaxes will allow us to do padding via OffsetArrays, which has the major advantage of keeping the indices of the padded array in-register with those of the original array. This should simplify a lot of code.
This makes several major changes: - it eliminates the generated functions and switches emphasis to the value domain rather than the type domain - it uses OffsetArrays to handle the padding. When using interpolate!, the axes may be narrowed to account for the boundary conditions. - now, the axes of the output correspond to the region where we can guarantee that we reconstruct the original array values - switches to the bounds-checking framework in base (may not be working yet)
Also add an analysis & plotting function
This makes it backwards-incompatible but will be more consistent for the future.
39f6dbf
to
144dc00
Compare
🎉 great! |
Here's
an early preview of where I'm heading with my attempteda rewrite. Some of the goals are laid out in #224. Despite adding more tests and more docstrings,so farthis isa tiny bitshorter than its predecessor, which may be a good sign.A couple things I'm particularly happy about:
BSpline(Quadratic(BC()))
) is considerably simpler and smaller than before. This may help attract user/developers.interpolate
, this always agrees with the parent, but withinterpolate!
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.interpolate
; useextrapolate
if you want to evaluate out-of-bounds.BSplineInterpolation
object. The good news is that when they are all empty, there is no storage overhead:There are a couple of broken tests, and I've temporarily commented outGridded
,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:
Related issues:
[]
#212