Support arbitrary indices for interpolation & extrapolation #146
Conversation
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
sglyon
reviewed
Apr 11, 2017
| ubound{T,N,TCoefs,IT}(itp::BSplineInterpolation{T,N,TCoefs,IT,OnGrid}, d) = size(itp, d) | ||
| lbound{T,N,TCoefs,IT}(itp::BSplineInterpolation{T,N,TCoefs,IT,OnCell}, d) = 0.5 | ||
| ubound{T,N,TCoefs,IT}(itp::BSplineInterpolation{T,N,TCoefs,IT,OnCell}, d) = size(itp, d)+0.5 | ||
| lbound{T,N,TCoefs,IT}(itp::BSplineInterpolation{T,N,TCoefs,IT,OnGrid}, d) = |
There was a problem hiding this comment.
Could it potentially help here to have something like firstindex instead of first(indices?
I suppose since this will usually be constructing Base.OneTo or a <:Range it probably isn't going to make a performance difference, but I'd be curious to see if special casing first and last indices would help us overcome the performance regression with 3d Constant BSplines
|
Hey @timholy this looks really cool, thanks for putting it together! |
|
Also RE clearness and performance My take is that code clarity is always a plus, but given the level of performance-driven metaprogramming we are doing here that ship might have already sailed... What helped me most to wrap my head around how this library works was helpful comments and docstrings. |
Prefiltering currently requires extending A_ldiv_B_md! for non1-array types. This seems necessary because our linear algebra infrastructure all assumes 1-based indices.
|
OK, I made julia> regs = regressions(judge(minimum(pr), minimum(master)))
1-element BenchmarkTools.BenchmarkGroup:
tags: []
"bsplines" => 3-element BenchmarkGroup([])
julia> leaves(regs)
15-element Array{Any,1}:
(Any["bsplines","constant","1d_Constant()_OnGrid_construct"],TrialJudgement(-0.86% => invariant))
(Any["bsplines","constant","1d_Constant()_OnGrid_use"],TrialJudgement(+9.85% => regression))
(Any["bsplines","constant","3d_Constant()_OnCell_use"],TrialJudgement(+15.51% => regression))
(Any["bsplines","constant","2d_Constant()_OnGrid_construct"],TrialJudgement(-1.48% => invariant))
(Any["bsplines","constant","1d_Constant()_OnCell_use"],TrialJudgement(+9.84% => regression))
(Any["bsplines","constant","2d_Constant()_OnCell_use"],TrialJudgement(+9.19% => regression))
(Any["bsplines","constant","2d_Constant()_OnGrid_use"],TrialJudgement(+8.91% => regression))
(Any["bsplines","constant","1d_Constant()_OnCell_construct"],TrialJudgement(-0.85% => invariant))
(Any["bsplines","constant","3d_Constant()_OnGrid_use"],TrialJudgement(+8.58% => regression))
(Any["bsplines","constant","2d_Constant()_OnCell_construct"],TrialJudgement(-0.27% => invariant))
(Any["bsplines","cubic","3d_Quadratic(Periodic())_OnGrid_use"],TrialJudgement(+5.29% => regression))
(Any["bsplines","linear","1d_Linear()_OnGrid_construct"],TrialJudgement(-1.05% => invariant))
(Any["bsplines","linear","1d_Linear()_OnCell_construct"],TrialJudgement(-0.01% => invariant))
(Any["bsplines","linear","2d_Linear()_OnCell_construct"],TrialJudgement(-1.99% => invariant))
(Any["bsplines","linear","2d_Linear()_OnGrid_construct"],TrialJudgement(-1.67% => invariant))
julia> imps = improvements(judge(minimum(pr), minimum(master)))
1-element BenchmarkTools.BenchmarkGroup:
tags: []
"bsplines" => 2-element BenchmarkGroup([])
julia> leaves(imps)
32-element Array{Any,1}:
(Any["bsplines","cubic","1d_Quadratic(Line())_OnGrid_construct"],TrialJudgement(-3.22% => invariant))
(Any["bsplines","cubic","1d_Quadratic(Flat())_OnGrid_construct"],TrialJudgement(-5.76% => improvement))
(Any["bsplines","cubic","2d_Quadratic(Free())_OnCell_use"],TrialJudgement(-6.13% => improvement))
(Any["bsplines","cubic","1d_Quadratic(Periodic())_OnCell_construct"],TrialJudgement(-2.10% => invariant))
(Any["bsplines","cubic","1d_Quadratic(Line())_OnCell_use"],TrialJudgement(-10.77% => improvement))
(Any["bsplines","cubic","1d_Quadratic(Line())_OnCell_construct"],TrialJudgement(-2.26% => invariant))
(Any["bsplines","cubic","1d_Quadratic(Free())_OnCell_construct"],TrialJudgement(-2.28% => invariant))
(Any["bsplines","cubic","1d_Quadratic(Free())_OnGrid_construct"],TrialJudgement(-2.11% => invariant))
(Any["bsplines","cubic","1d_Quadratic(Periodic())_OnGrid_construct"],TrialJudgement(-1.48% => invariant))
(Any["bsplines","cubic","1d_Quadratic(Flat())_OnCell_construct"],TrialJudgement(-1.61% => invariant))
(Any["bsplines","quadratic","1d_Quadratic(Line())_OnGrid_construct"],TrialJudgement(-3.17% => invariant))
(Any["bsplines","quadratic","3d_Quadratic(Line())_OnCell_construct"],TrialJudgement(-13.75% => improvement))
(Any["bsplines","quadratic","1d_Quadratic(Line())_OnCell_construct"],TrialJudgement(-4.16% => invariant))
(Any["bsplines","quadratic","2d_Quadratic(Reflect())_OnCell_construct"],TrialJudgement(-1.51% => invariant))
(Any["bsplines","quadratic","1d_Quadratic(Flat())_OnCell_construct"],TrialJudgement(-13.07% => improvement))
(Any["bsplines","quadratic","1d_Quadratic(Periodic())_OnGrid_use"],TrialJudgement(-9.20% => improvement))
(Any["bsplines","quadratic","1d_Quadratic(Reflect())_OnGrid_construct"],TrialJudgement(-3.19% => invariant))
(Any["bsplines","quadratic","1d_Quadratic(Flat())_OnGrid_construct"],TrialJudgement(-2.93% => invariant))
(Any["bsplines","quadratic","1d_Quadratic(Periodic())_OnCell_construct"],TrialJudgement(-2.79% => invariant))
(Any["bsplines","quadratic","1d_Quadratic(Line())_OnCell_use"],TrialJudgement(-9.12% => improvement))
(Any["bsplines","quadratic","3d_Quadratic(Flat())_OnCell_construct"],TrialJudgement(-6.30% => improvement))
(Any["bsplines","quadratic","2d_Quadratic(Flat())_OnCell_construct"],TrialJudgement(-1.43% => invariant))
(Any["bsplines","quadratic","3d_Quadratic(InPlaceQ())_OnCell_use"],TrialJudgement(-8.45% => improvement))
(Any["bsplines","quadratic","1d_Quadratic(InPlace())_OnCell_construct"],TrialJudgement(-10.25% => improvement))
(Any["bsplines","quadratic","1d_Quadratic(Free())_OnCell_construct"],TrialJudgement(-3.67% => invariant))
(Any["bsplines","quadratic","2d_Quadratic(InPlaceQ())_OnCell_use"],TrialJudgement(-10.74% => improvement))
(Any["bsplines","quadratic","3d_Quadratic(Reflect())_OnCell_construct"],TrialJudgement(-2.25% => invariant))
(Any["bsplines","quadratic","1d_Quadratic(Free())_OnGrid_construct"],TrialJudgement(-5.14% => improvement))
(Any["bsplines","quadratic","1d_Quadratic(Periodic())_OnGrid_construct"],TrialJudgement(-3.01% => invariant))
(Any["bsplines","quadratic","3d_Quadratic(Flat())_OnCell_use"],TrialJudgement(-8.46% => improvement))
(Any["bsplines","quadratic","3d_Quadratic(Reflect())_OnCell_use"],TrialJudgement(-7.23% => improvement))
(Any["bsplines","quadratic","1d_Quadratic(Reflect())_OnCell_construct"],TrialJudgement(-3.35% => invariant)) Since on average this is faster than master, I'll merge it once it gets through CI. |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
The master branch of Interpolations assumes that the input arrays start indexing at 1. Starting with julia-0.5, it's possible to support arrays that index over arbitrary
AbstractUnitRanges. This PR appears to support such arrays. While not tested yet, I presume the first consumer of this functionality will be ImageTransformations.jl, where we use non-1 based indices to keep track of where pixels "really" are.In theory, internally the most thorough way to do this would have been to rework how we apply padding: it might be somewhat more natural to pad a
1:nrange so that the internal padded version has indices0:n+1rather than1:n+2. (Among other things, that would ensure thatx = 2.2would mean the same "physical" location no matter whether one was referring to the original or padded array.) However, that would make Interpolations strictly dependent on a package like OffsetArrays, and I wasn't sure whether people wanted that. So I decided to keep our current schemes and modify Interpolations so that OffsetArrays would be supported but not required.One consequence of this decision is that any user of this functionality needs to extend
A_ldiv_B_md!for whatever non1-array types (s)he wants to use. This seems necessary because our linear algebra infrastructure all assumes 1-based indices.The point of having added benchmarks was to prepare for this PR, so...here are the regressions and improvements:
Most of this is probably noise, but I suspect the regressions on
Constant(esp. for 3d) are real. If anyone is concerned about that, I suspect the next step would be to turnmap_repeatinto a@generatedfunction. I'd be interested in knowing how folks feel about the "cleanliness" vs "performance" tradeoff in this case.