Skip to content
This repository has been archived by the owner on Oct 22, 2021. It is now read-only.

handling floating point errors #31

Open
wants to merge 4 commits into
base: master
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions src/dinic.jl
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@ function dinic_impl end

while true
augment = blocking_flow!(residual_graph, source, target, capacity_matrix, flow_matrix, P)
augment == 0 && break
is_zero(augment) && break
flow += augment
end
return flow, flow_matrix
Expand Down Expand Up @@ -80,7 +80,7 @@ function blocking_flow! end
end
end

flow == 0 && continue # Flow cannot be augmented along path
is_zero(flow) && continue # Flow cannot be augmented along path

v = target
u = bv
Expand Down
16 changes: 16 additions & 0 deletions src/maximum_flow.jl
Original file line number Diff line number Diff line change
Expand Up @@ -178,3 +178,19 @@ function maximum_flow(
end
return maximum_flow(flow_graph, source, target, capacity_matrix, algorithm)
end

"""
is_zero(value)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You could replace iszero with isapprox directly, which handles both integers and abstract floats correctly

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We already discussed of that point, unfortunately eps is not defined for integers

The default isapprox is not well suited for our problem. If isapprox(0.2+0.1, 0.3) is returning true, isapprox(0.2+0.1-0.3, 0.0) is returning false. In fact, x ≈ 0 is equivalent to x == 0.
As we compare to 0, rtol is just useless, so we must define an atol.
If we want the tolerance to rely on the type, we must separate Float and Integer case (as the default rtol do), because eps is not defined for Integer types.
This is not an absurd way o do it as rtol is defined quite similarly:

# default tolerance arguments
rtoldefault(::Type{T}) where {T<:AbstractFloat} = sqrt(eps(T))
rtoldefault(::Type{<:Real}) = 0
function rtoldefault(x::Union{T,Type{T}}, y::Union{S,Type{S}}, atol::Real) where {T<:Number,S<:Number}
    rtol = max(rtoldefault(real(T)),rtoldefault(real(S)))
    return atol > 0 ? zero(rtol) : rtol
end

See Base.isapprox
and floatfuncs lines 285-291 for the definition of rtol


Test if the value is equal to zero. It handles floating point errors.
"""
function is_zero end
function is_zero(
value::T
) where {T}
if isa(value,AbstractFloat)
return isapprox(value, 0, atol = sqrt(eps(T)))
else
return value == 0
end
end
5 changes: 2 additions & 3 deletions src/push_relabel.jl
Original file line number Diff line number Diff line change
Expand Up @@ -169,7 +169,6 @@ function relabel! end
nothing
end


"""
discharge!(residual_graph, v, capacity_matrix, flow_matrix, excess, height, active, count, Q)

Expand All @@ -189,11 +188,11 @@ function discharge! end
Q::AbstractVector # FIFO queue
)
for to in lg.outneighbors(residual_graph, v)
excess[v] == 0 && break
is_zero(excess[v]) && break
push_flow!(residual_graph, v, to, capacity_matrix, flow_matrix, excess, height, active, Q)
end

if excess[v] > 0
if excess[v] ≉ zero(excess[v])
if count[height[v] + 1] == 1
gap!(residual_graph, height[v], excess, height, active, count, Q)
else
Expand Down