Skip to content

Iteration of Interpolation Knots #397

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

Merged
merged 12 commits into from
Jan 12, 2021
Merged

Conversation

awadell1
Copy link
Contributor

@awadell1 awadell1 commented Jan 10, 2021

Overview of Change

This pull request addresses Issue #346 by adding knots(itp) method to get a knots iterator for interpolations and extrapolations.

Support for the following has been implemented / tests (Please let me know if I missed something)

  • Directional Extrapolation Boundary Conditions
  • Nd Interpolations / Extrapolations (Iteration is handled on each axis separately and combined with Iterators.product)
  • Non-equally spaced grid (i.e., GriddedInterpolation)

Request for Comments

This would be my first public pull request, so here's hoping it's decent :)

Currently, my test suite is passing (And hooked into the rest of the test suite), but if I missed a corner case or support for an interpolant, please let me know.

Examples

The behavior for 1 and Nd interpolations essentially just spits out the results of getknots but flattened

julia> using Interpolations;
julia> itp = interpolate(1:5, BSpline(Linear()));
julia> knots(itp) |> collect
5-element Array{Int64,1}:
 1
 2
 3
 4
 5

julia> itp = interpolate(rand(3,3), BSpline(Linear()));
julia> knots(itp) |> collect
3×3 Array{Tuple{Int64,Int64},2}:
 (1, 1)  (1, 2)  (1, 3)
 (2, 1)  (2, 2)  (2, 3)
 (3, 1)  (3, 2)  (3, 3)

Iteration over extrapolations will yield an infinite sequence of knots (For Periodic / Reflect) boundary conditions

julia> x = [1.0, 2.3, 2.4, 2.8, 3.0];
julia> etp = LinearInterpolation(x, x.^2, extrapolation_bc=Periodic());
julia> knots(etp) |> x -> Iterators.take(x, 8) |> collect
8-element Array{Float64,1}:
 1.0
 2.3
 2.4
 2.8
 3.0
 4.3
 4.4
 4.8

Next Steps (?)

I'm interested in a knot iterator like this as a stepping stone to implementing Gaussian Quadrature integration. The next step in that direction would be to add a knotsbetween(itp; start, stop) method to iterate over all knots between start and stop inclusively.

Is that a direction worth heading, or would that be out of scope for this package?

Best,

Alex

Test shows support for 1 and 2+ dimensions and extrapolation boundary
conditions (Periodic, Reflect)

Currently lacking support for the following:
 - Direction based boundary conditions
 - Reverse iteration
Prior version would return Int(Inf) for inifinite seqences of knots
    - Which errors as Int(::Float64) is inexact
    - Really should just not exist as a method
New version dispatches to _knot_length(::KnotIterator, ::IteratorSize)
so only Knots with a length will have a method

Also replaced Iterators.ProductInterator -> Iterators.product
    - Exact same (Later dispatches to the former)
    - Former is undocumented (ie. it's internal)
Ie. extrapolation_bc = ((Line(),Flat()), Flat()) now works

Reworked Construction, IteratorSize and IteratorEltype as a result
    - Fixed construction to handle Directional BC (And not keep
      Dispatching, due to the nested Tuple
    - Length is now undefined iff knots are infinite (As expected)

Fixed IteratorSize so that KnotIterator, and KnotIterator{T} (ie.
missing ET) were Base.Unknown not Base.HasLength (As could be HasLength
or IsInfinite)

Added more tests to check directional BC and iteration interface methods
Also fixed type errors (Knots should be Ints not Floats) in
test/iterate.jl (No impact on test results)
Split Periodic / Reflective Base.iterate methods, and dropped nextstate
    - Differences between Periodic / Reflect occur for both the current
      item and computing the next state -> Common method wasn't cutting
      it
    - New methods pass all existing test + new test looking at an
      unequal grid

Added FwdExtrapSpec macro to generate Union{ET,Tuple{BC,ET}} type specs
    - Could write it out, but it was getting cluttered
    - Happens pre-compile so no runtime penalty
@codecov
Copy link

codecov bot commented Jan 10, 2021

Codecov Report

Merging #397 (a88aa45) into master (318ebc8) will increase coverage by 0.58%.
The diff coverage is 97.91%.

Impacted file tree graph

@@            Coverage Diff             @@
##           master     #397      +/-   ##
==========================================
+ Coverage   81.26%   81.84%   +0.58%     
==========================================
  Files          22       23       +1     
  Lines        1329     1377      +48     
==========================================
+ Hits         1080     1127      +47     
- Misses        249      250       +1     
Impacted Files Coverage Δ
src/Interpolations.jl 75.20% <ø> (ø)
src/iterate.jl 97.91% <97.91%> (ø)

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 318ebc8...a88aa45. Read the comment docs.

Removed isnothing and only as not available in Julia 1.0.5 (Current LTS)
src/iterate.jl Outdated
KnotIterator(k::AbstractArray{T}, bc::ET) where {T,ET <: ExtrapDimSpec} = KnotIterator{T,ET}(k, bc)

const RepeatKnots = Union{Periodic,Reflect}
Base.IteratorSize(::Type{KnotIterator}) = Base.SizeUnknown()
Copy link
Collaborator

Choose a reason for hiding this comment

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

Rather than prefixing Base like this, perhaps we should explicitly import these in Interpolations.jl ?

src/iterate.jl Outdated
Comment on lines 19 to 21
KnotIterator(k::NTuple{N,AbstractArray}, bc::NTuple{N, ExtrapSpec}) where {N} = map(KnotIterator, k, bc)
KnotIterator(k::NTuple{N,AbstractArray}, bc::BoundaryCondition) where {N} = map(x -> KnotIterator(x, bc), k)
KnotIterator(k::AbstractArray{T}, bc::ET) where {T,ET <: ExtrapDimSpec} = KnotIterator{T,ET}(k, bc)
Copy link
Collaborator

Choose a reason for hiding this comment

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

Add tests

src/iterate.jl Outdated
Comment on lines 28 to 30
_knot_iter_size(::Type{<:BoundaryCondition}) = Base.HasLength()
_knot_iter_size(::Type{<:RepeatKnots}) = Base.IsInfinite()
_knot_iter_size(::Type{Tuple{RevBC, FwdBC}}) where {RevBC,FwdBC} = _knot_iter_size(FwdBC)
Copy link
Collaborator

Choose a reason for hiding this comment

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

Add tests

src/iterate.jl Outdated
Comment on lines 32 to 34
Base.length(iter::KnotIterator) = _knot_length(iter, Base.IteratorSize(iter))
_knot_length(iter::KnotIterator, ::Base.HasLength) = iter.nknots
Base.size(iter::KnotIterator) = (length(iter),)
Copy link
Collaborator

Choose a reason for hiding this comment

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

Add tests

src/iterate.jl Outdated
Comment on lines 36 to 37
Base.IteratorEltype(::Type{KnotIterator{T,ET}}) where {T,ET} = Base.HasEltype()
Base.eltype(::Type{KnotIterator{T,ET}}) where {T,ET} = T
Copy link
Collaborator

Choose a reason for hiding this comment

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

Add tests

src/iterate.jl Outdated
@@ -0,0 +1,138 @@

export knots
Copy link
Collaborator

Choose a reason for hiding this comment

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

Please move to Interpolations.jl file so the imports and exports are all together.

@mkitti
Copy link
Collaborator

mkitti commented Jan 10, 2021

Thank you for the pull request. The idea looks good. Ping me when the tests are passing please. Thanks.

Regarding knotsbetween(itp; start, stop) why should the arguments be keywords?

Also got rid of using Iterators.flatten for interpolations -> Now
everything use KnotIterators

Coverage just under 100%, but it's showing no hits for:

iterate(iter::KnotIterator{T,ET}) where {T,ET <: @FwdExtrapSpec RepeatKnots}

Which, if true, would mean the state was getting set to (1, zero(T))
somewhere else...

Added more comments to code and unit tests
Copy link
Collaborator

@mkitti mkitti left a comment

Choose a reason for hiding this comment

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

Great. This looks like it is almost mergeable. The only thing left is to shore up the documentation for the KnotIterator struct type and making the example code into a jldoctest per Documenter.jl.

A markdown file in the docs directory as well into the Documenter.jl docs would be appreciated.

src/iterate.jl Outdated
Comment on lines 64 to 73
```julia-repl
julia> etp = LinearInterpolation([1.0, 1.2, 2.3, 3.0], rand(4); extrapolation_bc=Periodic())
julia> Iterators.take(knots(etp), 5) |> collect
5-element Array{Float64,1}:
1.0
1.2
2.3
3.0
3.2
```
Copy link
Collaborator

Choose a reason for hiding this comment

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

macro FwdExtrapSpec(bc)
:( Union{$bc,Tuple{BoundaryCondition,$bc}} )
end
struct KnotIterator{T,ET <: ExtrapSpec}
Copy link
Collaborator

Choose a reason for hiding this comment

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

Add proper doc string for KnotIterator

src/iterate.jl Outdated
iterate(iter::KnotIterator, idx::Integer) = idx <= iter.nknots ? (iter.knots[idx], idx+1) : nothing

# For repeating knots state is the knot index + offset value
function iterate(iter::KnotIterator{T,ET}) where {T,ET <: @FwdExtrapSpec RepeatKnots}
Copy link
Collaborator

Choose a reason for hiding this comment

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

Codecov is saying this the only line left not covered by tests. I'm not sure why though.

@awadell1
Copy link
Contributor Author

Matt,

Yeah, not sure what's going on with the lack of coverage for that line. It has to get hit otherwise, state doesn't get initialized correctly. Plus line 99 gets hit, so 98 must get hit too?

Regarding knotsbetween(itp; start, stop) why should the arguments be keywords?

So my thought here was to mimic the syntax of range: range(start; length::Union{Integer,Nothing}=nothing, stop=nothing, step=nothing).

So the syntax would be:

  • knotsbetween(itp; start=x) Iterates over all knots greater than x
  • knotsbetween(itp; stop=x) Iterates over all knots less than x
  • knotsbetween(itp; start=x, stop=y) Iterates over all knots between x and y

Although maybe that should be rolled into knots?

Thank you for the fast review, I'm working through your last round of comments, and will ping you when it's done.

Best,

Alex

Plus added doctests to check that the examples are working
Julia seems to have difficult tracking coverage for macros -> Switch to
using a type alias to avoid this

JuliaLang/julia#36283

Still not having coverage issues, but they go away if we disable
inlining -> julia --inline=no ...
I thought it wasn't an issue, because running docs/make.jl would error
on the first time but not the 2nd.

But this just doesn't error
@awadell1
Copy link
Contributor Author

Locally I can show full coverage for iterate.jl if I run julia with inlining disabled:

julia --inline=no --project -e '\
		using Pkg, Coverage; \
		Pkg.test("Interpolations", coverage=true); \
		LCOV.writefile("lcov.info", process_folder());'

Doctests wise, it looks like everything is passing now. Do you have a script you use to build the docs locally? I think I have everything working correctly locally, aside from the following warnings:

julia --project make.jl
[ Info: SetupBuildDirectory: setting up build directory.
[ Info: Doctest: running doctests.
[ Info: ExpandTemplates: expanding markdown templates.
[ Info: CrossReferences: building cross-references.
[ Info: CheckDocument: running document checks.
┌ Warning: 1 docstring potentially missing:
│ 
│     Interpolations.l4_2d_cs
└ @ Documenter.DocChecks ~/.julia/packages/Documenter/PLD7m/src/DocChecks.jl:63
[ Info: Populate: populating indices.
[ Info: RenderDocument: rendering document.
[ Info: HTMLWriter: rendering HTML pages.
┌ Warning: invalid local image: unresolved path in control.md
│   link = "doc/images/parametric_spline.png"
└ @ Documenter.Writers.HTMLWriter ~/.julia/packages/Documenter/PLD7m/src/Writers/HTMLWriter.jl:1718
┌ Warning: Documenter could not auto-detect the building environment Skipping deployment.
└ @ Documenter ~/.julia/packages/Documenter/PLD7m/src/deployconfig.jl:41

@mkitti
Copy link
Collaborator

mkitti commented Jan 12, 2021

There are some other issues we need to fix with the documentation (not related with this PR). This looks good. If you're done, I'll merge this. Thanks. -Mark

@awadell1
Copy link
Contributor Author

Yeah, I'm done with the knots work. I'll open a separate pull request for adding knotsbetween when that's ready.

Best,
Alex

@mkitti mkitti merged commit b735c08 into JuliaMath:master Jan 12, 2021
@awadell1 awadell1 deleted the feature/iterate branch January 13, 2021 21:31
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants