never ending max flow algorithm

[etiennedeg]

Hi,

The default max-flow algorithm is stuck in an endless loop with this graph:

The capacities of the first layer are all 0.1 The capacities of the second layer are all 1 The capacities of the third layer are all 0

Here is how to reproduce it :

using LightGraphsFlows
import LightGraphs
const lg = LightGraphs

flow_graph = lg.DiGraph(8)
capacity_matrix = zeros(8,8)
flow_edges = [
(1,3,0.1),(1,4,0.1),(1,5,0.1),(3,6,1),(3,7,1),(3,8,1),
(4,6,1),(4,7,1),(4,8,1),
(5,6,1),(5,7,1),(5,8,1),(6,2,0),(7,2,0),
(8,2,0)
]
for e in flow_edges
    u, v, f = e
    lg.add_edge!(flow_graph, u, v)
    capacity_matrix[u,v] = f
end
println(collect(lg.edges(flow_graph))) 
flow,flow_matrix = maximum_flow(flow_graph, 1, 2, capacity_matrix)

[etiennedeg]

After debugging, it appears it is due to a floating point error...

[matbesancon]

Thanks for opening this up and sorry for the delay (on your other PR that is). If this is due to the simplex method, we don't have the hand on the actual algorithms, Float64 is often hardcoded and numerical precision can be an issue.

In the pure-Julia algorithms maybe we could do something, but the c = 0 case is tricky, because you have both inequalities, theoretically yielding f ==0 but as you mention, floating point ¯_(ツ)_/¯

[etiennedeg]

Thanks for your response. The floating point error is in the push relabel algorithm, when it check if the excess is 0 to deactivate a node. One way to workaround would be to convert the float numbers to rationals, one other way would be to add a tolerance (with an epsilon depending on the precision of the float)

[matbesancon]

Yes there is a ≈ operator in Julia which could be used exactly for that. If you can make a PR that would be great

[sbromberger]

Chiming in from vacation: would it be better to use isapprox explicitly and have the tolerance be a default argument that can be overridden, or does that make things too complicated?

[etiennedeg]

I thought to use isapprox( n, 0, eps(typeof(n)) I think this should do the work for almost any case, as we only compare against 0. There is no multiplications in the algorithm, so errors can't be amplified to much Having a fixed parameter would not consider the type of the capacity matrix.

[etiennedeg]

Here is a PR: #31

[etiennedeg]

I think floating points are less harmful on the other algorithms, but it might be useful to handle it anyway