Optimization: iFUB algorithm for unweighted graph diameter (#482)

Co-authored-by: Stefan Krastanov <[email protected]>

Author: LoveLow-Global

Project.toml

@@ -1,6 +1,6 @@
 name = "Graphs"
 uuid = "86223c79-3864-5bf0-83f7-82e725a168b6"
-version = "1.13.2"
+version = "1.13.3"
 
 [deps]
 ArnoldiMethod = "ec485272-7323-5ecc-a04f-4719b315124d"

src/distance.jl

@@ -95,6 +95,9 @@ end
 Given a graph and optional distance matrix, or a vector of precomputed
 eccentricities, return the maximum eccentricity of the graph.
 
+An optimizied BFS algorithm (iFUB) is used for unweighted graphs, both in [undirected](https://www.sciencedirect.com/science/article/pii/S0304397512008687) 
+and [directed](https://link.springer.com/chapter/10.1007/978-3-642-30850-5_10) cases.
+
 # Examples
 ```jldoctest
 julia> using Graphs
@@ -106,11 +109,127 @@ julia> diameter(path_graph(5))
 4
 ```
 """
+function diameter end
+
 diameter(eccentricities::Vector) = maximum(eccentricities)
-function diameter(g::AbstractGraph, distmx::AbstractMatrix=weights(g))
+
+diameter(g::AbstractGraph) = diameter(g, weights(g))
+
+function diameter(g::AbstractGraph, ::DefaultDistance)
+    if nv(g) == 0
+        return 0
+    end
+
+    connected = is_directed(g) ? is_strongly_connected(g) : is_connected(g)
+
+    if !connected
+        return typemax(Int)
+    end
+
+    return _diameter_ifub(g)
+end
+
+function diameter(g::AbstractGraph, distmx::AbstractMatrix)
     return maximum(eccentricity(g, distmx))
 end
 
+function _diameter_ifub(g::AbstractGraph{T}) where {T<:Integer}
+    nvg = nv(g)
+    out_list = [outneighbors(g, v) for v in vertices(g)]
+
+    if is_directed(g)
+        in_list = [inneighbors(g, v) for v in vertices(g)]
+    else
+        in_list = out_list
+    end
+
+    active = trues(nvg)
+    visited = falses(nvg)
+    queue = Vector{T}(undef, nvg)
+    distbuf = fill(typemax(T), nvg)
+    diam = 0
+
+    # Sort vertices by total degree (descending) to maximize pruning potential
+    vs = collect(vertices(g))
+    sort!(vs; by=v -> -(length(out_list[v]) + length(in_list[v])))
+
+    for u in vs
+        if !active[u]
+            continue
+        end
+
+        # Forward BFS from u
+        fill!(visited, false)
+        visited[u] = true
+        queue[1] = u
+        front = 1
+        back = 2
+        level_end = 1
+        e = 0
+
+        while front < back
+            v = queue[front]
+            front += 1
+
+            @inbounds for w in out_list[v]
+                if !visited[w]
+                    visited[w] = true
+                    queue[back] = w
+                    back += 1
+                end
+            end
+
+            if front > level_end && front < back
+                e += 1
+                level_end = back - 1
+            end
+        end
+        diam = max(diam, e)
+
+        # Backward BFS (Pruning)
+        dmax = diam - e
+
+        # Only prune if we have a chance to exceed the current diameter
+        if dmax >= 0
+            fill!(distbuf, typemax(T))
+            distbuf[u] = 0
+            queue[1] = u
+            front = 1
+            back = 2
+
+            while front < back
+                v = queue[front]
+                front += 1
+
+                if distbuf[v] >= dmax
+                    continue
+                end
+
+                @inbounds for w in in_list[v]
+                    if distbuf[w] == typemax(T)
+                        distbuf[w] = distbuf[v] + 1
+                        queue[back] = w
+                        back += 1
+                    end
+                end
+            end
+
+            # Prune vertices that cannot lead to a longer diameter
+            @inbounds for v in vertices(g)
+                if active[v] && distbuf[v] != typemax(T) && (distbuf[v] + e <= diam)
+                    active[v] = false
+                end
+            end
+        end
+
+        if !any(active)
+            break
+        end
+    end
+
+    return diam
+end
+
 """
     periphery(eccentricities)
     periphery(g, distmx=weights(g))

test/distance.jl

@@ -2,8 +2,12 @@
     g4 = path_digraph(5)
     adjmx1 = [0 1 0; 1 0 1; 0 1 0] # graph
     adjmx2 = [0 1 0; 1 0 1; 1 1 0] # digraph
+    adjmx3 = [0 1 0; 0 0 0; 0 0 0]
     a1 = SimpleGraph(adjmx1)
     a2 = SimpleDiGraph(adjmx2)
+    a3 = SimpleDiGraph(adjmx3)
+    a4 = blockdiag(complete_graph(5), complete_graph(5));
+    add_edge!(a4, 1, 6)
     distmx1 = [Inf 2.0 Inf; 2.0 Inf 4.2; Inf 4.2 Inf]
     distmx2 = [Inf 2.0 Inf; 3.2 Inf 4.2; 5.5 6.1 Inf]
 
@@ -44,6 +48,74 @@
             @test @inferred(center(z)) == center(g, distmx2) == [2]
         end
     end
+
+    @testset "Disconnected graph diameter" for g in test_generic_graphs(a3)
+        @test @inferred(diameter(g)) == typemax(Int)
+    end
+
+    @testset "simplegraph diameter" for g in test_generic_graphs(a4)
+        @test @inferred(diameter(g)) == 3
+    end
+
+    @testset "Empty graph diameter" begin
+        @test @inferred(diameter(SimpleGraph(0))) == 0
+        @test @inferred(diameter(SimpleDiGraph(0))) == 0
+    end
+
+    @testset "iFUB diameter" begin
+        # 1. Comparing against large graphs with known diameters  
+        n_large = 5000
+        g_path = path_graph(n_large)
+        @test diameter(g_path) == n_large - 1
+
+        g_cycle = cycle_graph(n_large)
+        @test diameter(g_cycle) == floor(Int, n_large / 2)
+
+        g_star = star_graph(n_large)
+        @test diameter(g_star) == 2
+
+        # 2. Comparing against the original implementation for random graphs
+        function diameter_naive(g)
+            return maximum(eccentricity(g))
+        end
+
+        NUM_SAMPLES = 50 # Adjust this to change test duration
+
+        for i in 1:NUM_SAMPLES
+            # Random unweighted Graphs
+            n = rand(10:1000) # Small to Medium size graphs
+            p = rand() * 0.1 + 0.005 # Sparse to medium density
+
+            # Undirected Graphs
+            g = erdos_renyi(n, p)
+            @test diameter(g) == diameter_naive(g)
+
+            ccs = connected_components(g)
+            largest_component = ccs[argmax(length.(ccs))]
+            g_lscc, _ = induced_subgraph(g, largest_component)
+
+            if nv(g_lscc) > 1
+                d_new = @inferred diameter(g_lscc)
+                d_ref = diameter_naive(g_lscc)
+                @test d_new == d_ref
+            end
+
+            # Directed Graphs
+            g_dir = erdos_renyi(n, p, is_directed=true)
+            @test diameter(g_dir) == diameter_naive(g_dir)
+
+            sccs = strongly_connected_components(g_dir)
+            largest_component_directed = sccs[argmax(length.(sccs))]
+            g_dir_lscc, _ = induced_subgraph(g_dir, largest_component_directed)
+
+            if nv(g_dir_lscc) > 1
+                d_new_dir = @inferred diameter(g_dir_lscc)
+                d_ref_dir = diameter_naive(g_dir_lscc)
+                @test d_new_dir == d_ref_dir
+            end
+        end
+    end
+
     @testset "DefaultDistance" begin
         @test size(Graphs.DefaultDistance()) == (typemax(Int), typemax(Int))
         d = @inferred(Graphs.DefaultDistance(3))