diff --git a/docs/make.jl b/docs/make.jl index 5ab4fba..53c162d 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -12,6 +12,7 @@ makedocs( "Multiroute flows" => "multiroute.md", "Min-cost flows" => "mincost.md", "Min-cut" => "mincut.md", + "Internals" => "internals.md", ] ) diff --git a/docs/src/index.md b/docs/src/index.md index 7953f65..83722b7 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -1,4 +1,4 @@ - # GraphsFlows.jl: flow algorithms for Graphs.jl +# GraphsFlows.jl: flow algorithms for Graphs.jl ```@meta CurrentModule = GraphsFlows @@ -8,12 +8,6 @@ DocTestSetup = quote end ``` -```@autodocs -Modules = [GraphsFlows] -Pages = ["GraphsFlows.jl"] -Order = [:function, :type] -``` - This is the documentation page for `GraphsFlows`. In all pages, we assume GraphsFlows has been imported into scope and that Graphs.jl has been imported. diff --git a/docs/src/internals.md b/docs/src/internals.md new file mode 100644 index 0000000..be13aa9 --- /dev/null +++ b/docs/src/internals.md @@ -0,0 +1,8 @@ +# Internals + +```@docs +GraphsFlows.DefaultCapacity +GraphsFlows.AbstractFlowAlgorithm +GraphsFlows.residual +GraphsFlows.is_zero +``` diff --git a/docs/src/maxflow.md b/docs/src/maxflow.md index 2dbaa97..5ae30a0 100644 --- a/docs/src/maxflow.md +++ b/docs/src/maxflow.md @@ -1,5 +1,13 @@ ## Max flow algorithms +```@docs +maximum_flow +EdmondsKarpAlgorithm +DinicAlgorithm +PushRelabelAlgorithm +BoykovKolmogorovAlgorithm +``` + ```@autodocs Modules = [GraphsFlows] Pages = ["maxflow.jl", "boykov_kolmogorov.jl", "push_relabel.jl", "dinic.jl", "edmonds_karp.jl"] diff --git a/src/dinic.jl b/src/dinic.jl index 1f000d4..249fa6c 100644 --- a/src/dinic.jl +++ b/src/dinic.jl @@ -11,8 +11,10 @@ function dinic_impl end residual_graph::::Graphs.IsDirected, # the input graph source::Integer, # the source vertex target::Integer, # the target vertex - capacity_matrix::AbstractMatrix{T} # edge flow capacities + capacity_matrix::AbstractMatrix{T}; # edge flow capacities + tolerance::T = (T <: AbstractFloat) ? sqrt(eps(T)) : zero(T) ) where {T} + n = Graphs.nv(residual_graph) # number of vertexes flow_matrix = zeros(T, n, n) # initialize flow matrix P = zeros(Int, n) # Sharable parent vector @@ -20,8 +22,8 @@ function dinic_impl end flow = 0 while true - augment = blocking_flow!(residual_graph, source, target, capacity_matrix, flow_matrix, P) - augment == 0 && break + augment = blocking_flow!(residual_graph, source, target, capacity_matrix, flow_matrix, P; tolerance = tolerance) + is_zero(augment; atol = tolerance) && break flow += augment end return flow, flow_matrix @@ -42,7 +44,8 @@ function blocking_flow! end target::Integer, # the target vertex capacity_matrix::AbstractMatrix{T}, # edge flow capacities flow_matrix::AbstractMatrix, # the current flow matrix - P::AbstractVector{Int} # Parent vector to store Level Graph + P::AbstractVector{Int}; # Parent vector to store Level Graph + tolerance = (T <: AbstractFloat) ? sqrt(eps(T)) : zero(T) ) where {T} n = Graphs.nv(residual_graph) # number of vertexes fill!(P, -1) @@ -80,7 +83,7 @@ function blocking_flow! end end end - flow == 0 && continue # Flow cannot be augmented along path + is_zero(flow; atol = tolerance) && continue # Flow cannot be augmented along path v = target u = bv @@ -108,11 +111,13 @@ blocking_flow( source::Integer, # the source vertex target::Integer, # the target vertex capacity_matrix::AbstractMatrix, # edge flow capacities - flow_matrix::AbstractMatrix, # the current flow matrix - ) = blocking_flow!( + flow_matrix::AbstractMatrix{T}; # the current flow matrix + tolerance = (T <: AbstractFloat) ? sqrt(eps(T)) : zero(T) + ) where {T}= blocking_flow!( residual_graph, source, target, capacity_matrix, flow_matrix, - zeros(Int, Graphs.nv(residual_graph))) + zeros(Int, Graphs.nv(residual_graph)); + tolerance = tolerance) diff --git a/src/maximum_flow.jl b/src/maximum_flow.jl index d8d5467..bcbe226 100644 --- a/src/maximum_flow.jl +++ b/src/maximum_flow.jl @@ -178,3 +178,11 @@ function maximum_flow( end return maximum_flow(flow_graph, source, target, capacity_matrix, algorithm) end + +""" + is_zero(value; tolerance) + +Test if the value is equal to zero. It handles floating point errors. +""" +is_zero(value::T; atol = sqrt(eps(T))) where {T<:AbstractFloat} = isapprox(value, zero(T), atol = atol) +is_zero(value; atol) = iszero(value) diff --git a/src/push_relabel.jl b/src/push_relabel.jl index 6aacc10..cf96e00 100644 --- a/src/push_relabel.jl +++ b/src/push_relabel.jl @@ -12,7 +12,8 @@ function push_relabel end residual_graph::::Graphs.IsDirected, # the input graph source::Integer, # the source vertex target::Integer, # the target vertex - capacity_matrix::AbstractMatrix{T} # edge flow capacities + capacity_matrix::AbstractMatrix{T}; # edge flow capacities + tolerance::T = (T <: AbstractFloat) ? sqrt(eps(T)) : zero(T) ) where {T} n = Graphs.nv(residual_graph) @@ -43,7 +44,7 @@ function push_relabel end while length(Q) > 0 v = pop!(Q) active[v] = false - discharge!(residual_graph, v, capacity_matrix, flow_matrix, excess, height, active, count, Q) + discharge!(residual_graph, v, capacity_matrix, flow_matrix, excess, height, active, count, Q; tolerance = tolerance) end return sum([flow_matrix[v, target] for v in Graphs.inneighbors(residual_graph, target)]), flow_matrix @@ -180,20 +181,21 @@ function discharge! end @traitfn function discharge!( residual_graph::::Graphs.IsDirected, # the input graph v::Integer, # vertex to be discharged - capacity_matrix::AbstractMatrix, + capacity_matrix::AbstractMatrix{T}, flow_matrix::AbstractMatrix, excess::AbstractVector, height::AbstractVector{Int}, active::AbstractVector{Bool}, count::AbstractVector{Int}, - Q::AbstractVector # FIFO queue - ) + Q::AbstractVector; # FIFO queue + tolerance = (T <: AbstractFloat) ? sqrt(eps(T)) : zero(T) + ) where {T} for to in Graphs.outneighbors(residual_graph, v) - excess[v] == 0 && break + is_zero(excess[v]; atol = tolerance) && break push_flow!(residual_graph, v, to, capacity_matrix, flow_matrix, excess, height, active, Q) end - if excess[v] > 0 + if ! is_zero(excess[v]; atol = tolerance) if count[height[v] + 1] == 1 gap!(residual_graph, height[v], excess, height, active, count, Q) else diff --git a/test/push_relabel.jl b/test/push_relabel.jl index 9eb0e9e..7e9c75b 100644 --- a/test/push_relabel.jl +++ b/test/push_relabel.jl @@ -91,4 +91,24 @@ 0 0 0 0 1 0] g448 = Graphs.DiGraph(M448) @test maximum_flow(g448, 1, 2, M448, algorithm=PushRelabelAlgorithm())[1] == 1 + + # Non regression test for floating point error (issue #28 in LightGraphsFlows.jl) + M28 = [0 0 1 1 1 0 0 0 + 0 0 0 0 0 0 0 0 + 0 0 0 0 0 1 1 1 + 0 0 0 0 0 1 1 1 + 0 0 0 0 0 1 1 1 + 0 1 0 0 0 0 0 0 + 0 1 0 0 0 0 0 0 + 0 1 0 0 0 0 0 0] + g28 = Graphs.DiGraph(M28) + C28 = [0. 0. 0.1 0.1 0.1 0. 0. 0. + 0. 0. 0. 0. 0. 0. 0. 0. + 0. 0. 0. 0. 0. 1. 1. 1. + 0. 0. 0. 0. 0. 1. 1. 1. + 0. 0. 0. 0. 0. 1. 1. 1. + 0. 0. 0. 0. 0. 0. 0. 0. + 0. 0. 0. 0. 0. 0. 0. 0. + 0. 0. 0. 0. 0. 0. 0. 0.] + @test maximum_flow(g28, 1, 2, C28, algorithm=PushRelabelAlgorithm())[1] == 0 end