Skip to content
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

Introducing iteration on vertices #1133

Open
wants to merge 10 commits into
base: master
Choose a base branch
from

Conversation

ghyatzo
Copy link

@ghyatzo ghyatzo commented Nov 12, 2024

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
Copy link
Member

juliohm commented Nov 13, 2024

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
Copy link
Author

ghyatzo commented Nov 13, 2024

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
Copy link
Member

juliohm commented Nov 13, 2024

@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
Copy link
Author

ghyatzo commented Nov 13, 2024

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
Copy link
Author

ghyatzo commented Nov 13, 2024

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
Copy link
Member

juliohm commented Nov 13, 2024

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
Copy link
Author

ghyatzo commented Nov 13, 2024

(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
Copy link
Member

juliohm commented Nov 13, 2024

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.

src/Meshes.jl Outdated Show resolved Hide resolved
@juliohm
Copy link
Member

juliohm commented Nov 13, 2024

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
Copy link
Author

ghyatzo commented Nov 13, 2024

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

@juliohm
Copy link
Member

juliohm commented Nov 13, 2024

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

@ghyatzo
Copy link
Author

ghyatzo commented Nov 13, 2024

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 ghyatzo marked this pull request as ready for review November 13, 2024 16:07
@ghyatzo
Copy link
Author

ghyatzo commented Nov 13, 2024

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.

src/domains/meshes.jl Outdated Show resolved Hide resolved
test/polytopes.jl Outdated Show resolved Hide resolved
@juliohm
Copy link
Member

juliohm commented Nov 13, 2024

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 and others added 2 commits November 13, 2024 17:30
Co-authored-by: Júlio Hoffimann <julio.hoffimann@gmail.com>
@ghyatzo
Copy link
Author

ghyatzo commented Nov 13, 2024

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
Copy link
Author

ghyatzo commented Nov 13, 2024

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 am not sure how to exclude the allocation of the object being iterated over, like one can do with @btime.
for example:

julia> rg = RegularGrid(point1, point2, (1, 1))
julia> @btime eachindex($rg)
  2.416 ns (0 allocations: 0 bytes)
julia> @btime iterate(eachindex($rg))
  3.041 ns (0 allocations: 0 bytes)

but

julia> @allocated for v in eachvertex(rg)
       end
227008
# this does not happen for the simple polytopes
julia> @allocated for v in eachvertex(t()) #<--- t() is a short helper that builds a random Triangle
       end
0

another example:

julia> @btime for _ in eachvertex_true($rg)
       end
  16.032 ns (0 allocations: 0 bytes)

julia> @btime for _ in eachvertex_true(rg)
       end
  60.000 μs (1478 allocations: 63.55 KiB)

@juliohm
Copy link
Member

juliohm commented Nov 13, 2024

The vertices of the RegularGrid subtypes are built on the fly upon demand. It is strange to see allocations. All vertices are stored in static memory, for instance.

@ghyatzo
Copy link
Author

ghyatzo commented Nov 13, 2024

Seems like all the allocations arrive from here:

function iterate(g::Generator, s...)
    @inline
    y = iterate(g.iter, s...)
    y === nothing && return nothing
    y = y::Tuple{Any, Any} # try to give inference some idea of what to expect about the behavior of the next line
    return (g.f(y[1]), y[2])   # <---  this line here in particular
end

which is the code for iterating generators in Base.

in particular, I've seen these allocations using @check_allock:

julia> try ttt() catch err; err.errors[1] end
Allocating runtime call to "jl_f_getfield" in ./Base.jl:49
  | getproperty(x, f::Symbol) = (@inline; getfield(x, f))

Stacktrace:
 [1] getproperty
   @ ./Base.jl:49 [inlined]
 [2] iterate
   @ ./generator.jl:48 [inlined]
 [3] var"##ttt#446"()
   @ Main ./REPL[190]:2


julia> try ttt() catch err; err.errors[2] end
Allocating runtime call to "jl_get_nth_field_checked" in ./tuple.jl:31
  | getindex(@nospecialize(t::Tuple), i::Int) = getfield(t, i, @_boundscheck)

Stacktrace:
 [1] getindex
   @ ./tuple.jl:31 [inlined]
 [2] iterate
   @ ./generator.jl:48 [inlined]
 [3] var"##ttt#446"()
   @ Main ./REPL[190]:2

julia> try ttt() catch err; err.errors[3] end
Dynamic dispatch in ./generator.jl:48
  | return (g.f(y[1]), y[2])

Stacktrace:
 [1] iterate
   @ ./generator.jl:48 [inlined]
 [2] var"##ttt#446"()
   @ Main ./REPL[190]:2


julia> try ttt() catch err; err.errors[4] end
Allocating runtime call to "jl_f_getfield" in ./Base.jl:49
  | getproperty(x, f::Symbol) = (@inline; getfield(x, f))

Stacktrace:
 [1] getproperty
   @ ./Base.jl:49 [inlined]
 [2] iterate
   @ ./generator.jl:45 [inlined]
 [3] var"##ttt#446"()
   @ Main ./REPL[190]:3


julia> try ttt() catch err; err.errors[5] end
Dynamic dispatch to function iterate in ./generator.jl:45
  | y = iterate(g.iter, s...)

Stacktrace:
 [1] iterate
   @ ./generator.jl:45 [inlined]
 [2] var"##ttt#446"()
   @ Main ./REPL[190]:3


julia> try ttt() catch err; err.errors[6] end
Allocating runtime call to "jl_f_getfield" in ./Base.jl:49
  | getproperty(x, f::Symbol) = (@inline; getfield(x, f))

Stacktrace:
 [1] getproperty
   @ ./Base.jl:49 [inlined]
 [2] iterate
   @ ./generator.jl:48 [inlined]
 [3] var"##ttt#446"()
   @ Main ./REPL[190]:3


julia> try ttt() catch err; err.errors[7] end
Allocating runtime call to "jl_get_nth_field_checked" in ./tuple.jl:31
  | getindex(@nospecialize(t::Tuple), i::Int) = getfield(t, i, @_boundscheck)

Stacktrace:
 [1] getindex
   @ ./tuple.jl:31 [inlined]
 [2] iterate
   @ ./generator.jl:48 [inlined]
 [3] var"##ttt#446"()
   @ Main ./REPL[190]:3


julia> try ttt() catch err; err.errors[8] end
Dynamic dispatch in ./generator.jl:48
  | return (g.f(y[1]), y[2])

Stacktrace:
 [1] iterate
   @ ./generator.jl:48 [inlined]
 [2] var"##ttt#446"()
   @ Main ./REPL[190]:3

not completely sure what is going on. Seems like the compiler is giving up in trying to infer stuff and fallbacks to runtime allocation. In this case f is vertex(m, ind)

@juliohm
Copy link
Member

juliohm commented Nov 14, 2024

That is very interesting. The vertex(grid, ind) should be fully static. From your investigations I understand that the issue is in the nested generator, and that your initial approach with a custom VertexIter eliminates these allocations.

If that is the case, we have a good motive to introduce the custom iterator in a new iterators.jl file, use it in the eachvertex methods defined in the corresponding geometry files, and include tests for full coverage.

Thank you for the deep analysis of memory usage, it really helps with these design choices. 💯

@ghyatzo
Copy link
Author

ghyatzo commented Nov 15, 2024

I am a doofus, I was checking the allocations using non const global variables, that is why it was struggling.
I've redid the testing locally using a function enclosure and it successfully doesn't allocate for most stuff now, even by just using generators.

julia> @allocated begin                                                                                                                                                                                                                 
           for v in eachvertex(T)   # const Triangle                                                                                                                                                                                                    
           end                                                                                                                                                                                                                          
       end                                                                                                                                                                                                                              
0                                                                                                                                                                                                                                       

julia> @allocated begin                                                                                                                                                                                                                 
           for v in eachvertex(H)  # const Hexagon                                                                                                                                                                                                     
           end                                                                                                                                                                                                                          
       end                                                                                                                                                                                                                              
0

julia> @allocated begin                                                                                                                                                                                                                 
           for v in eachvertex(M)  # const Multi([T, H])                                                                                                                                                                                                     
           end                                                                                                                                                                                                                          
       end                                                                                                                                                                                                                              
3312 # <--- this still allocates.  

julia> @allocated begin                                                                                                                                                                                                                 
           for v in eachvertex(RG)   # const RegularGrid                                                                                                                                                                                                   
           end                                                                                                                                                                                                                          
       end                                                                                                                                                                                                                              
0

julia> @allocated begin                                                                                                                                                                                                                 
           for v in eachvertex(SG)  # const StructuredGrid                                                                                                                                                                                                     
           end                                                                                                                                                                                                                          
       end                                                                                                                                                                                                                              
0

The issues still seems to be the compiler choking on the nested generator at some point.
I'll setup again the custom iterator to see if we can get rid of it!

@juliohm
Copy link
Member

juliohm commented Nov 15, 2024

But if the nested generators work well with non global variables, we are set, right? It is simpler code with the desired behavior.

@ghyatzo
Copy link
Author

ghyatzo commented Nov 15, 2024

Yeah, but the code I shared above is with constants, not global variables.
And it still allocates. Also inside a function barrier:

julia> forloops(g) = for v in eachvertex(g) end
julia> @allocated forloops(m)
3024

@juliohm
Copy link
Member

juliohm commented Nov 15, 2024 via email

@ghyatzo
Copy link
Author

ghyatzo commented Nov 16, 2024

I have been hacking away at it, but I can't seem to get rid of the type instability in the Multi iteration.
for example, using as an example m = Multi([T, H]) where T is a Triangle and H is a Hexagon.

the compiler correctly infers that m is of type Multi{M, C, <:Ngon}, but when iterating over the geometries it gives up on infering the type of each geometry. It knows that m.geoms is of type Vector{Ngon{N}} where N but it stops there.

Which incurs in dynamic dispatch calls to find the number of vertices through nvertices and therefore allocations (although much less than the generators).

I've tried multiple approaches but I fear that without adding some extra information in the structure itself to help the compiler
having complete inference of all elements inside the vector might be hard.

Couldn't find a way around it yet, let me know if you have any insights.

@juliohm
Copy link
Member

juliohm commented Nov 16, 2024 via email

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