Introducing iteration on vertices

[ghyatzo]

Hi,

following #1127 and #1131 I'm starting this draft PR.
I've added the main iteration logic for single geometry and multi geometry types.

This si incomplete as I would require some input on the following points

1. check any wrong assumption i've made so far
2. relevant types in the package
3. code placement (where in the source tree) and general code organization

in particular, I am not very familiar with all the nomenclature and type system structure used in this package, so I would appreciate an indication of which types we are interested in and/or for which onese we want to be able to call vertices.

WRT to point 2, I've gone for stuffing all iteration related code into a single file for now, let me know if this has to change.

cheers!

[juliohm]

Thank you @ghyatzo for the contribution. To facilitate the usage of this iterator, we just need to implement a new eachvertex function for Polytope Multi, and Mesh subtypes. Can you please refactor the PR to add this function instead?

For the Multi case you can use Iterators.flatten for example.

[ghyatzo]

Hi, yes absolutely can do.

Regarding flatten, I've done some tries, but the flatten iterator implementations I could come up with are both slower and worse allocating than the custom one here:

julia> @btime begin                                                                                            
           s = 0                                                                                               
           for _ in flat_multiitr                                                                              
               s+=1                                                                                            
           end                                                                                                 
           s                                                                                                   
       end                                                                                                     
  2.472 μs (45 allocations: 3.38 KiB)
9

julia> @btime begin                                                                                            
           s = 0                                                                                               
           for _ in custom_multiitr                                                                            
               s+=1                                                                                            
           end                                                                                                 
           s                                                                                                   
       end                                                                                                     
  1.167 μs (36 allocations: 1.27 KiB)
9

julia> @btime Base.iterate($flat_multiitr)                                                                     
  448.391 ns (5 allocations: 304 bytes)                                                                        
(Point(x: 0.011963944663157644 m, y: 0.13991304491415013 m), (2, VertexItr{Triangle{𝔼{2}, CoordRefSystems.Cartesian2D{CoordRefSystems.NoDatum, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}}}(Triangle((x: 0.0119639 m, y: 0.139913 m), (x: 0.347587 m, y: 0.958666 m), (x: 0.145339 m, y: 0.109886 m))), 2))          
                                                                                                               
julia> @btime Base.iterate($custom_multiitr)                                                                   
  53.077 ns (1 allocation: 32 bytes)
(Point(x: 0.011963944663157644 m, y: 0.13991304491415013 m), (1, 2))

flatten requires an iterator of iterators, which means that we have to build an iterator for each subelement in the multigeometry.
While the custom one can cheat and implement a pseudo iterator since we can manually call and it doesn't have to go through Base.iterate.
We could (i think) circumvent this by defining Base.iterate directly for the single elements, i.e. Base.iterate(p::Polytope), but that would get in the way of one day maybe be able to iterate over the edges or faces of a geometry instead of the vertices.
I think making a specific VertexIterator that acts on geometries is a more solid approach.

[juliohm]

@ghyatzo I meant that we simply need to define the function eachvertex(::Polytope) without additional constructs. It can be as simple as:

eachvertex(p::Polytope) = (v for v in vertices(p))

The geometry that will benefit the most is the Multi, which will avoid the allocations.

[ghyatzo]

The reason why I didn't use simply a generator (v for v in vertices(p)) is to make the Multi iterator (which is the one we care about) more performant. if we simply rely on generators it means that we have to instantiate a generator for each geometry.
while with the approach I went with you have a single iterator that returns the correct element at each step.
It is similar to how the iterator for an AbstractArray is the array itself, here we can use the geometry itself to be the iterator object, without the need to instantiate extra objetcs.

It should still be generic enough since underneat we are calling vertex(el, idx)

Please let me know if there is anything I overlooked or wrong with my reasoning.
I'll try to refactor it more to your preference in the mean time.

[ghyatzo]

I've improved the flatten approach using a generator of generators and not an array as I was naively doing before.
But I still see much more allocations and time than the custom iterator approach.

julia> eachvertex_custom(e::Union{Mesh,Polytope,Multi}) = VertexItr(e)
julia> eachvertex_gen(p::Polytope) = (v for v in vertices(p))
julia> eachvertex_flat(e::Multi) = Iterators.flatten(eachvertex_gen(g) for g in e.geoms)                       
eachvertex_flat (generic function with 1 method)
julia> m                                                                                                       
MultiNgon
├─ Triangle((x: 0.0119639 m, y: 0.139913 m), (x: 0.347587 m, y: 0.958666 m), (x: 0.145339 m, y: 0.109886 m))
└─ Hexagon((x: 0.25018 m, y: 0.85677 m), ..., (x: 0.57332 m, y: 0.3497 m))

julia> @btime begin                                                                                            
           s = 0                                                                                               
           for _ in eachvertex_flat($m)                                                                                                                                                                                       
               s+=1                                                                                            
           end                                                                                                 
           s                                                                                                   
       end                                                                                                     
  1.104 μs (39 allocations: 2.95 KiB)
9

julia> @btime begin                                                                                            
           s = 0                                                                                               
           for _ in eachvertex_custom($m)                                                                                                                                                                                     
               s+=1                                                                                            
           end                                                                                                 
           s                                                                                                   
       end                                                                                                     
  497.851 ns (9 allocations: 288 bytes)
9

those 9 allocations are from an Int64 allocation in the iterate function that If ironed out can become completely non allocating.

Anyway, I completely understand prefering the simpler solution.

Also: is it enough to just have the eachvertex function just for Polytopes, Multi and Mesh? let me know if I've missed something here as well.

[juliohm]

Isn't ok to define

eachvertex(m::Multi) = (v for g in m.geoms for v in vertices(g))

in the Multi geometry case? Why do we need nested generators?

[ghyatzo]

(v for g in m.geoms for v in vertices(g))

is shorthand for

Iterators.flatten((v for v in vertices(g)) for g in e.geoms)

from lowering the code:

julia> eachvertex_flatter(m::Multi) = (v for g in m.geoms for v in vertices(g))

julia> @code_lowered eachvertex_flatter(m)                                                                     
CodeInfo(
1 ─ %1 = Main.:(var"#13#14")
│        #13 = %new(%1)
│   %3 = #13
│   %4 = Base.getproperty(e, :geoms)
│   %5 = Base.Generator(%3, %4)
│   %6 = Base.Flatten(%5) <---- we construct a flatten iterator anyway
└──      return %6
)

and has similar performances

@btime begin                                                                                            
           s = 0                                                                                               
           for _ in eachvertex_flatter($m)                                                                                                                                                                                    
               s+=1                                                                                            
           end                                                                                                 
           s                                                                                                   
       end                                                                                                     
  1.092 μs (39 allocations: 2.95 KiB)
9

[juliohm]

is it enough to just have the eachvertex function just for Polytopes, Multi and Mesh?

Yes, but please restrict to MultiPolytope instead of general Multi geometries.

(v for g in m.geoms for v in vertices(g))

is shorthand for

Iterators.flatten((v for v in vertices(g)) for g in e.geoms)

Interesting. Thanks for clarifying this.

[juliohm]

Appreciate if you can keep it simple for now without an auxiliary iterator struct. Please define the eachvertex function in the respective geometry files, and add new tests with @allocated.

We can reconsider the auxiliary iterator struct in the future. I understand it allocates a fixed amount no matter the size of the MultiPolytope.

[ghyatzo]

Alright, will do. shall I make a separate PR and keep this as draft with the struct approach?

[juliohm]

Feel free to update the PR here directly. You can copy the VertexIter struct in the comments to save it for future reference.

[ghyatzo]

For future reference:

struct VertexItr{T}
    el::T
end

_v_iterate(el::Polytope, i) =
    (@inline; (i - 1) % UInt < nvertices(el) % UInt ? (@inbounds vertex(el, i), i + 1) : nothing)

_v_iterate(el::Mesh, i) =
    (@inline; (i - 1) % UInt < nvertices(el) % UInt ? (@inbounds vertex(el, i), i + 1) : nothing)

Base.iterate(itr::VertexItr{<:Mesh}, i=1) = _v_iterate(itr.el, i)

Base.iterate(itr::VertexItr{<:Polytope}, i=1) = _v_iterate(itr.el, i)

Base.iterate(itr::VertexItr{<:MultiPolytope}, state=(1, 1)) = begin
    ig, ivg = state
    ig > length(itr.el.geoms) && return nothing

    is = _v_iterate(itr.el.geoms[ig], ivg)
    is === nothing && return Base.iterate(itr, (ig + 1, 1))

    v, ivg = is
    return (v, (ig, ivg))
end

Base.length(itr::VertexItr{<:Mesh}) = nvertices(T)
Base.length(itr::VertexItr{<:Polytope}) = nvertices(T)
Base.length(itr::VertexItr{<:MultiPolytope}) = sum(nvertices, itr.el.geoms)

Base.IteratorSize(itr::VertexItr) = Base.HasLength()
Base.eltype(::VertexItr) = Point

eachvertex(e::T) where {T} = error("Vertex iterator not implemented for type $T")
eachvertex(e::Union{Mesh,Polytope,MultiPolytope}) = VertexItr(e)

[ghyatzo]

I would like some pointers WRT the new tests, I'm a bit unsure where best to put them
and what you mean exactly when you meant with @allocated.

I've also took the liberty of rewriting some vertices() function to instead collecting the iterator instead of repeating the same thing. Let me know if you agree with the change.

[juliohm]

I would like some pointers WRT the new tests, I'm a bit unsure where best to put them

You can add tests in the corresponding test files. For instance, test/polytopes.jl, test/multigeoms.jl and test/meshes.jl.

and what you mean exactly when you meant with @allocated.

I shared a comment above about that. Just a test that uses the eachvertex in a for loop, making sure that the for loop doesn't allocate. You could even nest for loops to make sure that everything is working as expected, motivated by the original discussion on Discourse.

I've also took the liberty of rewriting some vertices() function to instead collecting the iterator instead of repeating the same thing.

Does it lead to the same type inference? I know that list comprehensions are very efficient, and that the compiler can do magic with them. I am not sure if this is the case with the collect of generator.

[ghyatzo]

Does it lead to the same type inference? I know that list comprehensions are very efficient, and that the compiler can do magic with them. I am not sure if this is the case with the collect of generator.

I'll do some tests, but I expect it does, since I have the sneaking suspicion that doing [.. for ... in ...] is just syntax sugar for collect((... for ... in ...))

Edit:
Indeed

julia> Meta.@lower [x for x in 1:10]
:($(Expr(:thunk, CodeInfo(
    @ none within `top-level scope`
1 ─ %1 = 1:10
│   %2 = Base.Generator(Base.identity, %1)
│   %3 = Base.collect(%2)
└──      return %3
))))

[ghyatzo]

ok accepted everything, can we close this up now?

[eliascarv]

Not yet, I will add more tests before merging.

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 87.71%. Comparing base (243d3f0) to head (6c5e330).
Report is 4 commits behind head on master.

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1133      +/-   ##
==========================================
+ Coverage   87.69%   87.71%   +0.02%     
==========================================
  Files         191      191              
  Lines        6005     6007       +2     
==========================================
+ Hits         5266     5269       +3     
+ Misses        739      738       -1     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

---

🚨 Try these New Features:

• Flaky Tests Detection - Detect and resolve failed and flaky tests

[juliohm]

Thank you @ghyatzo for initiating this work, and @eliascarv for finishing it up!

We will work on the final detail now with MultiPolyArea in a separate PR to make all cases allocation-free.