Performance improvements and bugfix to Tarjan's algorithm

[saolof]

This is the same as this PR in the old lightgraphs repository.

From substantial benchmarking, this is a significant speedup in most cases and an asymptotic improvement from quadratic to linear for star graphs, with a performance regression of a few percent only in the limit of extremely sparse graphs. Furthermore, it opens up the possibility for future speedups in other functions provided by this package, since it also computes data that can be directly used in those other functions.

This PR does not change the API for the function itself, but it does add a couple of performance and type hint functions that other library types that inherit from AbstractGraphs can optionally overload to improve/tune performance, for example if they store adjacency lists in a linked list and so have different performance characteristics than array representations.

[codecov]

Codecov Report

Merging #32 (b1018bf) into master (11f54ad) will increase coverage by 2.07%.
The diff coverage is 88.33%.

❗ Current head b1018bf differs from pull request most recent head 83cda6f. Consider uploading reports for the commit 83cda6f to get more accurate results

@@            Coverage Diff             @@
##           master      #32      +/-   ##
==========================================
+ Coverage   97.26%   99.33%   +2.07%     
==========================================
  Files         114      106       -8     
  Lines        6579     5571    -1008     
==========================================
- Hits         6399     5534     -865     
+ Misses        180       37     -143     

[ViralBShah]

@saolof While I can't review this in great detail, I am wondering what the level of testing we have is. If it is reasonable, then passing tests would make it easier to quickly merge.

[saolof]

@saolof While I can't review this in great detail, I am wondering what the level of testing we have is. If it is reasonable, then passing tests would make it easier to quickly merge.

Tests are reasonably all-encompassing. One gotcha is that there are a few probabilistic tests in the codebase that will fail in ~1% of test runs, and can cause tests to fail for unrelated code changes. So in those cases it is worth checking which test failed.

For this particular PR, it passes the CI tests, and in addition to that the original PR thread also has a ton of benchmarks which act as additional tests on a wide range of different random graphs to check that the output is the same as on stable.

[ViralBShah]

@jpfairbanks Thoughts?

[jpfairbanks]

I also haven't reviewed in detail, but I don't see anything crazy. @simonschoelly, you were reviewing this on the old repo. Can you take a look before we merge?

[ViralBShah]

Bump. Is this good to merge? Do we have enough tests to give us confidence in merging?

[oscardssmith]

@saolof can you add a few more tests to this so the new codepaths are fully tested? It looks like this PR currently decreases test coverage. Other than that, this is looks like a really nice improvement.

[saolof]

Sure. I'll add tests for infer_nb_iterstate_type by creating a mock graph type for it to infer the node iterator type of.

[gdalle]

Bump, has anyone been able to review this in detail?

[oscardssmith]

not yet.

[gdalle]

Presumably fails because of a bug in Aqua: JuliaTesting/Aqua.jl#139

[gdalle]

Aqua complains because this is ill-defined: Union{A,B} is the same as Union{B,A}

[gdalle]

This seems needlessly convoluted, is it really needed?

[saolof]

The thing that is actually needed is the infer_nb_iterstate_type function below. You could just make that return Any for the non-AbstractSimpleGraph case, and overload the is_large_vertex function to a much higher threshold to avoid regressions for user-defined graph types, and leave improved type inference for a later commit.

That could still mess with type stability though (relying on whomever defines a third party graph type to overload the functions to get the full performance benefit), while this solution seems to always be type stable even for arbitrary user-defined types. Another option is to rewrite this into a recursive function that uses the call stack for large nodes, but it's tricky to do that with no regressions.

[saolof]

The thing that the function actually does is fairly simple though, it just unpacks the type of the objects returned by an iterator. Referring to that type would be straightforward to do in any statically typed language where iterators have a typed interface instead of being a duck typed protocol, but is slightly tricky in Julia

[etiennedeg]

outneighbors should return vertices of the graph, which should be of type eltype(g), why do you need to use infer_nb_iterstate_type method, where you can just return T?

[saolof]

eltype is the first tuple element, I need the second (usually an Int, but could be a linked list node for example), so that when outneighbours is suspended I can restart from the same location instead of restarting from the beginning and making the algorithm quadratic in the number of edges.

See https://julialang.org/blog/2018/07/iterators-in-julia-0.7/ on the julia iteration protocol where the type I need is the type of the state object in the tuple

[gdalle]

Do we have correctness and performance tests? It seems this is the only file that changed in the PR

[saolof]

The original commit this was based on had fairly extensive benchmarking and testing in the thread, on a fairly wide range of graph types. I make no guarentees as to whether those still hold two years later after possible rebases/merges though

[gdalle]

Would it be possible to put those tests in the Graphs.jl test suite?

[saolof]

The tests relied on checking that the output was exactly the same as the old function:

println("Constructing large random graphs...")
example_graphs_labels =  ["path_digraph(10000)", "random_regular_digraph(50000,3)", "random_regular_digraph(50000,8)", "random_regular_digraph(50000,200)"," random_orientation_dag(random_regular_graph(50000,200))", "random_regular_digraph(50000,2000)"," random_tournament_digraph(10000)"," random_orientation_dag(complete_graph(10000))"," star_digraph(100000)"]
example_graphs =         [ path_digraph(10000),   random_regular_digraph(50000,3),   random_regular_digraph(50000,8),   random_regular_digraph(50000,200),   random_orientation_dag(random_regular_graph(50000,200)) ,   random_regular_digraph(50000,2000),   random_tournament_digraph(10000),   random_orientation_dag(complete_graph(10000)),   star_digraph(100000)]
println("Benchmarking:")

for f in (strongly_connected_components,strongly_connected_components_2)
    println(" ------------------------------")
    println("testing function $(string(f)): ")

    for i in eachindex(example_graphs) 
        gg = example_graphs[i]
        println(example_graphs_labels[i])

        @assert sort.(f(gg)) == sort.(strongly_connected_components(gg)) "incorrect implementation"
        @btime ($f($(gg));)
    end
end

If you still have a separate Kosaraju implementation, one possible test could be checking that the two give the same output.

[etiennedeg]

I just fixed a similar quadratic hit in #266. Rather than pushing neighbors on a separate stack (here, depending on number of neighbors), I did it by always pushing the neighbors directly on the dfs stack, and using a more refined coloring scheme than just visited / unvisited. Maybe it is worth doing the same here ?

[saolof]

I just fixed a similar quadratic hit in #266. Rather than pushing neighbors on a separate stack (here, depending on number of neighbors), I did it by always pushing the neighbors directly on the dfs stack, and using a more refined coloring scheme than just visited / unvisited. Maybe it is worth doing the same here ?

I don't push neighbours onto a stack at all. That would be quadratic memory. The algorithm only has a dfs stack for visited stack, and the iteration index for neighbours when iterating over nodes that have a lot of neighbours.

[etiennedeg]

For me, this is pushing on a stack ?

    if is_large_vertex(g, s)
        push!(largev_iterstate_stack, iterate(outneighbors(g, s)))
    end

This is not quadratic memory, this is linear in the number of edges.
Edit: Oh, you just push the iteration state

[saolof]

For me, this is pushing on a stack ?

    if is_large_vertex(g, s)
        push!(largev_iterstate_stack, iterate(outneighbors(g, s)))
    end

This is not quadratic memory, this is linear in the number of edges.

I am pushing the current node on the stack, not the neighbours. In the Julia iteration protocol, iterate returns a tuple of the current iterated object and the iteration index.

In a worst case memory scenario like a tournament graph where all nodes have many neighbours (most of which are already visited) that pushes all nodes but not all edges, and in that case this PR gives a 4x improvement in runtime compared to the previous implementation.

Tarjan's algorithm is implicitly worst case linear memory in the number of nodes, and among other things this PR improves the constant factor compared to the previous implementation.

[etiennedeg]

Ok, I think I have the solution: We need to use Stateful iterators. This boils down to the same, but at least, this is part of the julia API, and it provides a cleaner implementation.