Skip to content
This repository was 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

handling floating point errors #31

Changes from 1 commit
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
1fdf02a
handling floating point errors
EtiDg Sep 18, 2019
cc87014
handles floating point on dinic's algorithm
etiennedeg Nov 25, 2019
fe2aeff
wider tolerance for floating point
etiennedeg Nov 26, 2019
60f2f8a
Update src/push_relabel.jl
matbesancon Oct 12, 2020
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
19 changes: 17 additions & 2 deletions src/push_relabel.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -169,6 +169,21 @@ function relabel! end
nothing
end

"""
is_zero(value)

Test if the value is equal to zero. It handles floating point errors.
"""
function is_zero end
Copy link
Member

Choose a reason for hiding this comment

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

This function is not necessary, you can use isapprox(x, 0) or x ≈ 0 in the code

Copy link
Member Author

Choose a reason for hiding this comment

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

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

Copy link

Choose a reason for hiding this comment

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

Would set_zero_subnormals() help here?

Copy link
Member Author

Choose a reason for hiding this comment

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

I don't think so, 0.1+0.2-0.3 is much greater than subnormal numbers

Copy link
Member Author

Choose a reason for hiding this comment

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

After experimentation, comparing to the epsilon is not sufficient for large float.
It seems sqrt(epsilon) is a common choice for tolerance, and it is justified by mathematical error analysis.
Maybe an even better choice would be to choose the tolerance depending on the capacity matrix values. For example, tol = sqrt(eps(maximum(capacity_matrix) or tol = sqrt(eps(sum(capacity_matrix)

function is_zero(
value::T
) where {T}
if isa(value,AbstractFloat)
return isapprox(value, 0, atol = eps(T))
else
return value == 0
end
end

"""
discharge!(residual_graph, v, capacity_matrix, flow_matrix, excess, height, active, count, Q)
Expand All @@ -189,11 +204,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 ! is_zero(excess[v])
if count[height[v] + 1] == 1
gap!(residual_graph, height[v], excess, height, active, count, Q)
else
Expand Down