handling floating point errors #31
Conversation
Codecov Report
@@ Coverage Diff @@
## master #31 +/- ##
===========================================
+ Coverage 65.65% 98.26% +32.61%
===========================================
Files 10 10
Lines 428 289 -139
===========================================
+ Hits 281 284 +3
+ Misses 147 5 -142
Continue to review full report at Codecov.
|
src/push_relabel.jl
Outdated
|
|
||
| Test if the value is equal to zero. It handles floating point errors. | ||
| """ | ||
| function is_zero end |
There was a problem hiding this comment.
This function is not necessary, you can use isapprox(x, 0) or x ≈ 0 in the code
There was a problem hiding this comment.
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
There was a problem hiding this comment.
Would set_zero_subnormals() help here?
There was a problem hiding this comment.
I don't think so, 0.1+0.2-0.3 is much greater than subnormal numbers
There was a problem hiding this comment.
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)
|
I also run into an infinite loop with Dinic's algorithm. |
Did you fix this in this PR or is it still in progress? |
Yes, this is fixed with the new commit |
|
I forgot the sqrt for the epsilon tolerance. |
| @@ -178,3 +178,19 @@ function maximum_flow( | |||
| end | |||
| return maximum_flow(flow_graph, source, target, capacity_matrix, algorithm) | |||
| end | |||
|
|
|||
| """ | |||
| is_zero(value) | |||
There was a problem hiding this comment.
You could replace iszero with isapprox directly, which handles both integers and abstract floats correctly
There was a problem hiding this comment.
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 returningtrue,isapprox(0.2+0.1-0.3, 0.0)is returningfalse. In fact,x ≈ 0is equivalent tox == 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 endSee Base.isapprox
and floatfuncs lines 285-291 for the definition of rtol
|
@etienneINSA what's the status on this? |
|
Sorry for the delay, I think this can be merged, if you agree on the fact that we need a separate function to compare the excess against 0 |
I don't know how to use Trait.
This code should pass the tests for both integers and floats