fix floating point (#3)

Author: etiennedeg

docs/make.jl

@@ -12,6 +12,7 @@ makedocs(
         "Multiroute flows"   => "multiroute.md",
         "Min-cost flows"     => "mincost.md",
         "Min-cut"            => "mincut.md",
+        "Internals"          => "internals.md",
     ]
 )
 

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.
 

docs/src/internals.md

@@ -0,0 +1,8 @@
+# Internals
+
+```@docs
+GraphsFlows.DefaultCapacity
+GraphsFlows.AbstractFlowAlgorithm
+GraphsFlows.residual
+GraphsFlows.is_zero
+```

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"]

src/dinic.jl

@@ -11,17 +11,19 @@ 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
 
     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)

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)

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

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