Sbromberger/new shortest paths (#1242)

* first cut at aggregation

* abstractified paths

* bellman-ford

* A*, and removing getproperty

* floyd-warshall

* D'Esopo-Pape

* submodule strategy

* export ShortestPaths

* Update src/Experimental/ShortestPaths/ShortestPaths.jl

Co-Authored-By: Mathieu Besançon <[email protected]>

* Johnson

* spfa and bfs

* Update johnson.jl (#1241)

clean up docs for johnson

* tests

* failing tests

* fixed tests

* dijksta/desopo and more tests

* Result

* tests

* docs

* distmx is not untyped.

* exports and doc fix

* doc updates

* DFS algorithm

* refactor BFS shortest paths for NOOPSort

src/Experimental/Experimental.jl

@@ -8,11 +8,13 @@ export description,
     VF2, vf2, IsomorphismProblem, SubgraphIsomorphismProblem, InducedSubgraphIsomorphismProblem,
     could_have_isomorph, has_isomorph, all_isomorph, count_isomorph,
     has_induced_subgraphisomorph, count_induced_subgraphisomorph, all_induced_subgraphisomorph,
-    has_subgraphisomorph, count_subgraphisomorph, all_subgraphisomorph
+    has_subgraphisomorph, count_subgraphisomorph, all_subgraphisomorph,
 
+    ShortestPaths
 description() = "This module contains experimental graph functions."
 
 include("isomorphism.jl")
 include("vf2.jl") # Julian implementation of VF2 algorithm
+include("ShortestPaths/ShortestPaths.jl")
 
 end

src/Experimental/ShortestPaths/ShortestPaths.jl

@@ -0,0 +1,203 @@
+module ShortestPaths
+using SparseArrays: sparse
+using LightGraphs
+using LightGraphs: AbstractGraph, AbstractEdge
+using LightGraphs.SimpleGraphs: AbstractSimpleGraph
+
+using DataStructures:PriorityQueue, enqueue!, dequeue!
+
+
+# TODO: figure out how we keep environmental params.
+# struct LGEnvironment
+#     threaded::Bool
+#     parallel::Bool
+#     LGEnvironment() = new(false, false)
+# end
+
+abstract type AbstractGraphResult end
+abstract type AbstractGraphAlgorithm end
+
+"""
+    ShortestPathResult <: AbstractGraphResult
+
+Concrete subtypes of `ShortestPathResult` contain the
+results of a shortest-path calculation using a specific
+[`ShortestPathAlgorithm`](@ref).
+
+In general, the fields in these structs should not be
+accessed directly; use the [`dists`](@ref) and
+[`paths`](@ref) functions to obtain the results of the
+calculation.
+"""
+abstract type ShortestPathResult <: AbstractGraphResult end
+
+
+"""
+    ShortestPathAlgorithm <: AbstractGraphAlgorithm
+
+Concrete subtypes of `ShortestPathAlgorithm` are used to specify
+the type of shortest path calculation used by [`shortest_paths`](@ref).
+Some concrete subtypes (most notably [`Dijkstra`](@ref) have fields
+that specify algorithm parameters.
+
+See [`AStar`](@ref), [`BellmanFord`](@ref), [`BFS`](@ref),
+[`DEspopoPape`](@ref), [`Dijkstra`](@ref), [`FloydWarshall`](@ref),
+[`Johnson`](@ref), and [`SPFA`](@ref) for specific requirements and
+usage details.
+"""
+abstract type ShortestPathAlgorithm <: AbstractGraphAlgorithm end
+
+include("astar.jl")
+include("bellman-ford.jl")
+include("bfs.jl")
+include("desopo-pape.jl")
+include("dijkstra.jl")
+include("floyd-warshall.jl")
+include("johnson.jl")
+include("spfa.jl")
+
+
+################################
+# Shortest Paths via algorithm #
+################################
+# if we don't pass in distances but specify an algorithm, use weights.
+"""
+    shortest_paths(g, s, distmx, alg)
+    shortest_paths(g, s, t, alg)
+    shortest_paths(g, s, alg)
+    shortest_paths(g, s)
+    shortest_paths(g)
+
+Return a `ShortestPathResult` that allows construction of the shortest path
+between sets of vertices in graph `g`. Depending on the algorithm specified,
+other information may be required: (e.g., a distance matrix `distmx`, and/or
+a target vertex `t`). Some algorithms will accept multiple source vertices
+`s`; algorithms that do not accept any source vertex `s` produce all-pairs
+shortest paths. 
+
+See `ShortestPathAlgorithm` for more details on the algorithm specifications.
+
+### Implementation Notes
+The elements of `distmx` may be of any type that has a [Total Ordering](https://en.m.wikipedia.org/wiki/Total_order)
+and valid comparator, `zero` and `typemax` functions. Concretely, this means that
+distance matrices containing complex numbers are invalid.
+
+### Examples
+```
+g = path_graph(4)
+w = zeros(4, 4)
+for i in 1:3
+    w[i, i+1] = 1.0
+    w[i+1, i] = 1.0
+end
+
+s1 = shortest_paths(g)                   # `alg` defaults to `FloydWarshall`
+s2 = shortest_paths(g, 1)                # `alg` defaults to `BFS`
+s3 = shortest_paths(g, 1, w)             # `alg` defaults to `Dijkstra`
+s4 = shortest_paths(g, 1, BellmanFord())
+s5 = shortest_paths(g, 1, w, DEsopoPape())
+```
+"""
+shortest_paths(g::AbstractGraph, s, alg::ShortestPathAlgorithm) =
+    shortest_paths(g, s, weights(g), alg)
+
+# If we don't specify an algorithm AND there are no dists, use BFS.
+shortest_paths(g::AbstractGraph{T}, s::Integer) where {T<:Integer} = shortest_paths(g, s, BFS())
+shortest_paths(g::AbstractGraph{T}, ss::AbstractVector) where {T<:Integer} = shortest_paths(g, ss, BFS())
+
+# Full-formed methods.
+"""
+    paths(state[, v])
+    paths(state[, vs])
+    paths(state[, v, d]))
+
+Given the output of a [`shortest_paths`](@ref) calculation of type
+[`ShortestPathResult`](@ref), return a vector (indexed by vertex)
+of the paths between the source vertex used to compute the shortest path
+and a single destination vertex `v`, a vector of destination vertices `vs`,
+or the entire graph.
+
+For multiple destination vertices, each path is represented by a vector of
+vertices on the path between the source and the destination. Nonexistent
+paths will be indicated by an empty vector. For single destinations, the
+path is represented by a single vector of vertices, and will be length 0
+if the path does not exist.
+
+For [`ShortestPathAlgorithm`](@ref)s that compute all shortest paths for all
+pairs of vertices, `paths(state)` will return a vector (indexed by source
+vertex) of vectors (indexed by destination vertex) of paths. `paths(state, v)`
+will return a vector (indexed by destination vertex) of paths from source `v`
+to all other vertices. In addition, `paths(state, v, d)` will return a vector
+representing the path from vertex `v` to vertex `d`.
+"""
+function paths(state::ShortestPathResult, vs::AbstractVector{<:Integer})
+    parents = state.parents
+    T = eltype(parents)
+
+    num_vs = length(vs)
+    all_paths = Vector{Vector{T}}(undef, num_vs)
+    for i = 1:num_vs
+        all_paths[i] = Vector{T}()
+        index = T(vs[i])
+        if parents[index] != 0 || parents[index] == index
+            while parents[index] != 0
+                pushfirst!(all_paths[i], index)
+                index = parents[index]
+            end
+            pushfirst!(all_paths[i], index)
+        end
+    end
+    return all_paths
+end
+
+paths(state::ShortestPathResult, v::Integer) = paths(state, [v])[1]
+paths(state::ShortestPathResult) = paths(state, [1:length(state.parents);])
+
+"""
+    dists(state[, v])
+
+Given the output of a [`shortest_paths`](@ref) calculation of type
+[`ShortestPathResult`](@ref), return a vector (indexed by vertex)
+of the distances between the source vertex used to compute the
+shortest path and a single destination vertex `v` or the entire graph.
+
+For [`ShortestPathAlgorithm`](@ref)s that compute all-pairs shortest
+paths, `dists(state)` will return a matrix (indexed by source and destination
+vertices) of distances.
+"""
+dists(state::ShortestPathResult, v::Integer) = state.dists[v]
+dists(state::ShortestPathResult) = state.dists
+
+"""
+    has_negative_weight_cycle(g[, distmx=weights(g), alg=BellmanFord()])
+
+Given a graph `g`, an optional distance matrix `distmx`, and an optional
+algorithm `alg` (one of [`BellmanFord`](@ref) or [`SPFA`](@ref)), return
+`true` if any cycle detected in the graph has a negative weight.
+# Examples
+
+```jldoctest
+julia> g = complete_graph(3);
+
+julia> d = [1 -3 1; -3 1 1; 1 1 1];
+
+julia> has_negative_weight_cycle(g, d)
+true
+
+julia> g = complete_graph(4);
+
+julia> d = [1 1 -1 1; 1 1 -1 1; 1 1 1 1; 1 1 1 1];
+
+julia> has_negative_weight_cycle(g, d, SPFA())
+false
+```
+"""
+has_negative_weight_cycle(g::AbstractGraph, distmx::AbstractMatrix=weights(g)) = has_negative_weight_cycle(g, distmx, BellmanFord())
+has_negative_weight_cycle(g::AbstractSimpleGraph) = false
+
+export ShortestPathAlgorithm
+export paths, dists, shortest_paths, has_negative_weight_cycle
+export Dijkstra, AStar, BellmanFord, FloydWarshall, DEsopoPape, Johnson, SPFA, BFS
+export NegativeCycleError
+
+end  # module

src/Experimental/ShortestPaths/astar.jl

@@ -0,0 +1,125 @@
+# Parts of this code were taken / derived from Graphs.jl. See LICENSE for
+# licensing details.
+
+# A* shortest-path algorithm
+
+"""
+    struct AStar <: ShortestPathAlgorithm
+
+The structure used to configure and specify that [`shortest_paths`](@ref)
+should use the [A* search algorithm](http://en.wikipedia.org/wiki/A%2A_search_algorithm).
+An optional `heuristic` function may be supplied. If missing, the heuristic is set to
+`n -> 0`.
+
+### Implementation Notes
+`AStar` supports the following shortest-path functionality:
+- non-negative distance matrices / weights
+- single destination
+"""
+struct AStar{F<:Function} <: ShortestPathAlgorithm
+    heuristic::F
+end
+
+AStar(T::Type=Float64) = AStar(n -> zero(T))
+
+struct AStarResult{T, U<:Integer} <: ShortestPathResult
+    path::Vector{U}
+    dist::T
+end
+
+function reconstruct_path!(total_path, # a vector to be filled with the shortest path
+    came_from, # a vector holding the parent of each node in the A* exploration
+    end_idx, # the end vertex
+    g) # the graph
+
+    curr_idx = end_idx
+    while came_from[curr_idx] != curr_idx
+        pushfirst!(total_path, came_from[curr_idx])
+        curr_idx = came_from[curr_idx]
+    end
+    push!(total_path, end_idx)
+end
+
+function a_star_impl!(g, # the graph
+    goal, # the end vertex
+    open_set, # an initialized heap containing the active vertices
+    closed_set, # an (initialized) color-map to indicate status of vertices
+    g_score, # a vector holding g scores for each node
+    f_score, # a vector holding f scores for each node
+    came_from, # a vector holding the parent of each node in the A* exploration
+    distmx,
+    heuristic)
+
+    T = eltype(g)
+    total_path = Vector{T}()
+
+    @inbounds while !isempty(open_set)
+        current = dequeue!(open_set)
+
+        if current == goal
+            reconstruct_path!(total_path, came_from, current, g)
+            return total_path
+        end
+
+        closed_set[current] = true
+
+        for neighbor in outneighbors(g, current)
+            closed_set[neighbor] && continue
+
+            tentative_g_score = g_score[current] + distmx[current, neighbor]
+
+            if tentative_g_score < g_score[neighbor]
+                g_score[neighbor] = tentative_g_score
+                priority = tentative_g_score + heuristic(neighbor)
+                enqueue!(open_set, neighbor, priority)
+                came_from[neighbor] = current
+            end
+        end
+    end
+    return total_path
+end
+
+function calc_dist(path, distmx)
+    T = eltype(distmx)
+    isempty(path) && return typemax(T)
+    rest_of_path = copy(path)
+    dist = zero(T)
+    u = pop!(rest_of_path)
+    while !isempty(rest_of_path)
+        v = pop!(rest_of_path)
+        dist += distmx[u, v]
+        u = v
+    end
+    return dist
+end
+
+function shortest_paths(g::AbstractGraph, s::Integer, t::Integer, distmx::AbstractMatrix, alg::AStar)
+    T = eltype(distmx)
+
+    # if we do checkbounds here, we can use @inbounds in a_star_impl!
+    checkbounds(distmx, Base.OneTo(nv(g)), Base.OneTo(nv(g)))
+
+    open_set = PriorityQueue{Integer, T}()
+    enqueue!(open_set, s, 0)
+
+    closed_set = zeros(Bool, nv(g))
+
+    g_score = fill(Inf, nv(g))
+    g_score[s] = 0
+
+    f_score = fill(Inf, nv(g))
+    f_score[s] = alg.heuristic(s)
+
+    came_from = -ones(Integer, nv(g))
+    came_from[s] = s
+
+    path = a_star_impl!(g, t, open_set, closed_set, g_score, f_score, came_from, distmx, alg.heuristic)
+    return AStarResult(path, calc_dist(path, distmx))
+end
+
+shortest_paths(g::AbstractGraph, s::Integer, t::Integer, alg::AStar) = shortest_paths(g, s, t, weights(g), alg)
+paths(s::AStarResult) = [s.path]
+paths(s::AStarResult, v::Integer) = throw(ArgumentError("AStar produces at most one path."))
+dists(s::AStarResult) = [[s.dist]]
+dists(s::AStarResult, v::Integer) = throw(ArgumentError("AStar produces at most one path."))
+