Johnson Shortest Path for Sparse Graphs (#884)

* Johnson Shortest Path for Sparse Graphs

Johnson Shortest Path for Sparse Graphs

* Improved memory efficiency if distmx is mutable

* Improved memory efficiency for parallel implementation

src/LightGraphs.jl

@@ -97,9 +97,9 @@ MaximumAdjacency, AbstractMASVisitor, mincut, maximum_adjacency_visit,
 
 # a-star, dijkstra, bellman-ford, floyd-warshall
 a_star, dijkstra_shortest_paths, bellman_ford_shortest_paths,
-has_negative_edge_cycle, enumerate_paths, floyd_warshall_shortest_paths,
-transitiveclosure!, transitiveclosure, transitivereduction, yen_k_shortest_paths,
-parallel_multisource_dijkstra_shortest_paths,
+has_negative_edge_cycle, enumerate_paths, johnson_shortest_paths,
+floyd_warshall_shortest_paths, transitiveclosure!, transitiveclosure, transitivereduction, 
+yen_k_shortest_paths, parallel_multisource_dijkstra_shortest_paths,
 
 # centrality
 betweenness_centrality, closeness_centrality, degree_centrality,
@@ -217,6 +217,7 @@ include("edit_distance.jl")
 include("shortestpaths/astar.jl")
 include("shortestpaths/bellman-ford.jl")
 include("shortestpaths/dijkstra.jl")
+include("shortestpaths/johnson.jl")
 include("shortestpaths/floyd-warshall.jl")
 include("shortestpaths/yen.jl")
 include("linalg/LinAlg.jl")

src/shortestpaths/johnson.jl

@@ -0,0 +1,104 @@
+
+"""
+    struct JohnsonState{T, U}
+An [`AbstractPathState`](@ref) designed for Johnson shortest-paths calculations.
+"""
+struct JohnsonState{T<:Real, U<:Integer} <: AbstractPathState
+    dists::Matrix{T}
+    parents::Matrix{U}
+end
+
+@doc_str """
+    johnson_shortest_paths(g, distmx=weights(g); parallel=false)
+
+### Implementation Notes
+Use the [Johnson algorithm](https://en.wikipedia.org/wiki/Johnson%27s_algorithm)
+to compute the shortest paths between all pairs of vertices in graph `g` using an
+optional distance matrix `distmx`.
+If the parameter parallel is set true, dijkstra_shortest_paths will run in parallel.
+Parallel bellman_ford_shortest_paths is currently unavailable
+Return a [`LightGraphs.JohnsonState`](@ref) with relevant
+traversal information.
+Behaviour in case of negative cycle depends on bellman_ford_shortest_paths.
+Throws NegativeCycleError() if a negative cycle is present.
+### Performance
+Complexity: O(|V|*|E|)
+If distmx is not mutable or of type, DefaultDistance than a sparse matrix will be produced using distmx.
+In the case that distmx is immutable, to reduce memory overhead,  
+if edge (a, b) does not exist in g then distmx[a, b] should be set to 0.
+### Dependencies from LightGraphs
+bellman_ford_shortest_paths
+parallel_multisource_dijkstra_shortest_paths
+dijkstra_shortest_paths
+"""
+function johnson_shortest_paths(
+    g::AbstractGraph{U},
+    distmx::AbstractMatrix{T} = weights(g);
+    parallel::Bool = false
+) where T<:Real where U<:Integer
+
+    nvg = nv(g)
+    type_distmx = typeof(distmx)
+    #Change when parallel implementation of Bellman Ford available
+    wt_transform = bellman_ford_shortest_paths(g, vertices(g), distmx).dists
+
+    if !type_distmx.mutable && type_distmx !=  LightGraphs.DefaultDistance
+        distmx = sparse(distmx) #Change reference, not value
+    end
+
+    #Weight transform not needed if all weights are positive.
+    if type_distmx !=  LightGraphs.DefaultDistance
+        for e in edges(g)
+        	distmx[src(e), dst(e)] += wt_transform[src(e)] - wt_transform[dst(e)] 
+        end
+    end
+
+    if !parallel
+        dists = Matrix{T}(undef, nvg, nvg)
+        parents = Matrix{U}(undef, nvg, nvg)
+    	for v in vertices(g)
+    		dijk_state = dijkstra_shortest_paths(g, v, distmx)
+    		dists[v, :] = dijk_state.dists
+      		parents[v, :] = dijk_state.parents
+    	end
+    else
+    	dijk_state = parallel_multisource_dijkstra_shortest_paths(g, vertices(g), distmx)
+    	dists = dijk_state.dists
+    	parents = dijk_state.parents
+    end
+
+    broadcast!(-, dists, dists, wt_transform)
+    for v in vertices(g)
+        dists[:, v] .+= wt_transform[v] #Vertical traversal prefered
+    end
+
+    if type_distmx.mutable
+        for e in edges(g)
+            distmx[src(e), dst(e)] += wt_transform[dst(e)] - wt_transform[src(e)]
+        end
+    end
+
+    return JohnsonState(dists, parents)
+end
+
+function enumerate_paths(s::JohnsonState{T, U}, v::Integer) where T<:Real where U<:Integer
+    pathinfo = s.parents[v, :]
+    paths = Vector{Vector{U}}()
+    for i in 1:length(pathinfo)
+        if (i == v) || (s.dists[v, i] == typemax(T))
+            push!(paths, Vector{U}())
+        else
+            path = Vector{U}()
+            currpathindex = i
+            while currpathindex != 0
+                push!(path, currpathindex)
+                currpathindex = pathinfo[currpathindex]
+            end
+            push!(paths, reverse(path))
+        end
+    end
+    return paths
+end
+
+enumerate_paths(s::JohnsonState) = [enumerate_paths(s, v) for v in 1:size(s.parents, 1)]
+enumerate_paths(st::JohnsonState, s::Integer, d::Integer) = enumerate_paths(st, s)[d]

test/runtests.jl

@@ -39,6 +39,7 @@ tests = [
     "shortestpaths/astar",
     "shortestpaths/bellman-ford",
     "shortestpaths/dijkstra",
+    "shortestpaths/johnson",
     "shortestpaths/floyd-warshall",
     "shortestpaths/yen",
     "traversals/bfs",

test/shortestpaths/johnson.jl

@@ -0,0 +1,44 @@
+@testset "Johnson" begin
+    g3 = PathGraph(5)
+    d = LinearAlgebra.Symmetric([0 1 2 3 4; 1 0 6 7 8; 2 6 0 11 12; 3 7 11 0 16; 4 8 12 16 0])
+    for g in testgraphs(g3)
+      z = @inferred(johnson_shortest_paths(g, d))
+      @test z.dists[3, :][:] == [7, 6, 0, 11, 27]
+      @test z.parents[3, :][:] == [2, 3, 0, 3, 4]
+
+      @test @inferred(enumerate_paths(z))[2][2] == []
+      @test @inferred(enumerate_paths(z))[2][4] == enumerate_paths(z, 2)[4] == enumerate_paths(z, 2, 4) == [2, 3, 4]
+
+      z = @inferred(johnson_shortest_paths(g, d, parallel=true))
+      @test z.dists[3, :][:] == [7, 6, 0, 11, 27]
+      @test z.parents[3, :][:] == [2, 3, 0, 3, 4]
+
+      @test @inferred(enumerate_paths(z))[2][2] == []
+      @test @inferred(enumerate_paths(z))[2][4] == enumerate_paths(z, 2)[4] == enumerate_paths(z, 2, 4) == [2, 3, 4]
+
+    end
+
+    g4 = PathDiGraph(4)
+    for g in testdigraphs(g4)
+      z = @inferred(johnson_shortest_paths(g))
+      @test length(enumerate_paths(z, 4, 3)) == 0
+      @test length(enumerate_paths(z, 4, 1)) == 0
+      @test length(enumerate_paths(z, 2, 3)) == 2
+
+      z = @inferred(johnson_shortest_paths(g, parallel=true))
+      @test length(enumerate_paths(z, 4, 3)) == 0
+      @test length(enumerate_paths(z, 4, 1)) == 0
+      @test length(enumerate_paths(z, 2, 3)) == 2
+    end 
+
+    g5 = DiGraph([1 1 1 0 1; 0 1 0 1 1; 0 1 1 0 0; 1 0 1 1 0; 0 0 0 1 1])
+    d = [0 3 8 0 -4; 0 0 0 1 7; 0 4 0 0 0; 2 0 -5 0 0; 0 0 0 6 0]
+    for g in testdigraphs(g5)
+      z = @inferred(johnson_shortest_paths(g, d))
+      @test z.dists == [0 1 -3 2 -4; 3 0 -4 1 -1; 7 4 0 5 3; 2 -1 -5 0 -2; 8 5 1 6 0]
+
+      z = @inferred(johnson_shortest_paths(g, d, parallel=true))
+      @test z.dists == [0 1 -3 2 -4; 3 0 -4 1 -1; 7 4 0 5 3; 2 -1 -5 0 -2; 8 5 1 6 0]
+    end 
+
+end