Approximate equality

[stevengj]

@jlperla

So in that case I think the warning should also be to call it with a ≈ b for vectors if the infinity norm is good enough

Actually, it defaults to the Euclidean norm (vecnorm) for arrays, and you can pass a different norm if desired via the norm keyword. It only falls back to elementwise isapprox if the norm is Inf or NaN. This is all documented. You don’t have to guess at the behavior.

We shouldn’t warn people against using this for arrays — if anything, we should warn people against implementing their own approximate equality for arrays, because experience suggests they will probably get it wrong. (e.g. NumPy got it wrong: their isclose function is only elementwise, has a nonzero default atol, and even has the property that isclose(x,y) is inequivalent to isclose(y,x)!)

I don’t think it should throw a DomainError for x ≈ 0. This is perfectly well defined, it just means x == 0 since the definition is norm(x-y) ≤ rtol*max(norm(x),norm(y)) for atol==0. It would be very weird for 0 ≈ 0 to throw an error rather than returning the correct answer (true), and this would be an especially unfortunate behavior for elementwise isapprox comparisons.

[Stephen_Vavasis]

When I proposed earlier that isapprox should be dropped from the language, I wasn’t thinking about CI testing. I agree with all the respondents who said that isapprox is useful for testing, so I’ll change my proposal: isapprox should be moved out of Base and into Base.Test. Consider for example Stefan Karpinski’s snippet above showing how FFTW testing uses isapprox. Suppose that I am an undergrad at a lesser school than MIT, and my professor tells me to check in Julia whether ifft(fft(x)) is equal to x, but to be careful of roundoff error. So I look at how FFTW’s testing checks whether ifft(fft(x)==x holds in the presence of roundoff, and I crib that code. See below. My professor is wrong: ifft(fft(x)) is not equal to x even accounting for roundoff error.

After reading my example, you might say that bone-heads deserve what they get, but I would say that the language shouldn’t attempt to put training wheels on a Ferrari.

julia> x = exp.(4*(0:63)*2*pi*1im/64);
julia> y = exp.(8*(0:63)*2*pi*1im/64);
julia> z = ifft(fft(x-y));
julia> for i = 1 : length(x)
       print(isapprox(x[i]-y[i], z[i]), " ")
       end
false true true true true true true true true true true true true true true true false true true true true true true true true true true true true true true true false true true true true true true true true true true true true true true true false true true true true true true true true true true true true true true true

[jlperla]

I completely agree on the suggested move to Base.Test. This sort of thing shouldn’t be used in normal algorithms (and, hence, can have a lot of potentially inefficient checks and warnings).

IMy suggested implementation should have already had the 0=0 check before any warnings, so 0 ≈ 0 would be true. Also realized that it should check to make sure that atol==0 before any warning or domain error as well.

The norm on vectors is great, and that would explain why I was having trouble finding the code (and still can’t) I assume it finds the relative tolerance of a ≈ b though taking norm(a) and norm(b)?

Then things would work well. Basically, instead of abs(x-y) <= atol + rtol*max(abs(x), abs(y))) being the test, it would be something like
norm(x-y) <= atol + rtol*max(norm(x), norm(y)))
If so, then it fits great into my suggestion as well. If norm(a) or norm(b) is identically zero, then it would spit out a similar warning as the scalar one.

Whether it is a domain error or not, I think a warning is necessary. I have heard many great explanations for why every suggested default tolerance for a ≈ 0 doesn’t work, but I haven’t heard a single example of where it makes any sense, or represents an intent of the user… except in the knife-edge 0 ≈ 0 case which is equivalently 0 == 0.

But just to defend the function: it is extremely useful, especially for the less informed writing tests who haven’t thought through relative tolerances (I consider myself in that group).

Edit: it might also be useful to export the calculation of the tolerance itself, as it is simple but extremely useful for newer users. Something along the lines of (copying code from floatfuncs)

#Directly from the code.
rtoldefault(::Type{T}) where {T<:AbstractFloat} = sqrt(eps(T))
rtoldefault(::Type{<:Real}) = 0
rtoldefault(x::Union{T,Type{T}}, y::Union{S,Type{S}}) where {T<:Number,S<:Number} = max(rtoldefault(real(T)),rtoldefault(real(S)))

#Copying the inside of isapprox:
relativetol(x, y, ; rtol::Real=rtoldefault(x,y), atol::Real=0) = atol + rtol*max(abs(x), abs(y)))

[ChrisRackauckas]

@jlperla

but I haven’t heard a single example of where it makes any sense, or represents an intent of the user

It’s strict not lenient. It will accidentally lead to test failures instead of test passing. So if you put it on a test and it fails, then that informs you that you should set tolerances to something you think is reasonable. An error when zero is quite arbitrary, whereas for pretty much any test you need to take a second look at your test and set tolerances if it doesn’t pass, which is basically what needs to be done in any case anyways.

[jlperla]

What are you opposed to a warning to tell (likely less technical) people that they are doing something invalid almost everywhere? Consider the undergrads I have planned to write up tests for DiffEqOperators.jl

My gut tells me the solution is NOT to set tolerances for most cases, but rather to reorganize the abs(a-b) ≈ 0 to a ≈ b and rely on those carefully designed tolerances. But I could be wrong.

[ChrisRackauckas]

@jlperla

My gut tells me the solution is NOT to set tolerances for most cases, but rather to reorganize the abs(a-b) ≈ 0 to a ≈ b and rely on those carefully designed tolerances. But I could be wrong.

Yes, that is correct. That’s the only thing that makes sense with floating point unless you have reason to believe abs(a-b) should be sufficiently close by some absolute tolerance, in which case you can place that tolerance. The issue isn’t the testing setup, it’s the unintuitiveness of floating point numbers. The relative differences between floating point numbers tends to stay the same, which means that the absolute differences between floating point numbers doesn’t stay the same.

julia> eps()
2.220446049250313e-16

That’s the smallest value for which 1+eps() != 1. But note that

julia> eps(0.1)
1.3877787807814457e-17

When dealing with numbers like 0.1, the relative spacing of the floating point numbers is the same (1e-16 less) but the absolute difference between numbers (the size of floating point errors) is smaller.

So the phrase “the same up to floating point error” is ill-defined: it depends on the size of the values you’re talking about. That’s why this phrase has a non-zero definition for a ≈ b but, since you loose information about relative size, abs(a-b) ≈ 0 has a definition of “must be zero” since the absolute differences between floating point numbers goes to zero near zero

julia> eps(0.0)
5.0e-324

Another way to understand this is to realize that half of all floating point numbers reside in [-1,1]. This stuff is just a reality of floating point and does need to be directly thought about when writing tests.

@jlperla

Consider the undergrads I have planned to write up tests for DiffEqOperators.jl

They will likely need guidance. Floating point errors is the least of the problems. Understanding how to properly test that algorithms are converging with a correct order of accuracy, and regression testing on nonlinear problems without analytical solutions, is an art that requires good knowledge of the problem you’re solving. Finding the correct tests is as much of a problem as implementing the algorithms. For example, for DiffEqOperators the idea will have to exploit the local linearity of the problem: testing Fokker-Plank solutions directly is a bad idea (except for integration tests) but you can account for the behaviors against analytical solutions of the heat equation and the advection equation and piece those together to give complete unit tests. This is probably not something that would be obvious to a student who thinks the project is to write something that solves the Fokker-Plank equation.

[Tamas_Papp]

@jlperla

My suggested implementation should have already had the 0=0 check before any warnings, so 0 ≈ 0 would be true. Also realized that it should check to make sure that atol==0 before any warning or domain error as well.

You would replace a very clear and simple definition with a complicated one that has at least two branches to exceptions that the programmer needs to think and reason about. Also, consider how this would affect comparison of containers, eg arrays.

[jlperla]

@Tamas_Papp

consider how this would affect comparison of containers, eg arrays.

If it is just a warning, there is only a single branch? As for containers, I think it fits perfectly. Warn if
max(norm(x), norm(y))) == 0 && abstol == 0 for whatever norm is passed in. Or could do max(abs(x),abs(y)) == 0 && abstol == 0 for the Number case.

If there is disagreement with a DomainError then fine. My order of preferences is:

1. Warning
2. Explicit and clear statement in the docs that you should almost NEVER do x\approx0 and should attempt to rewrite tests in a x\approxy form.
3. Nothing, where people figure there is something wrong when x\approx0 tests fail, but they then pick some tolerance out of the sky and don’t understand the issue.

Keep in mind that I am probably an above average programmer for my field, and I have never understood the whole default tolerance thing until this thread. A warning can avoid you guys having to answer a whole bunch of these sorts of questions.

[Tamas_Papp]

@jlperla

I have never understood the whole default tolerance thing until this thread

I find the documentation of isapprox very clear, but if you don’t think it is, you could make a PR extending it.

[JeffreySarnoff]

imo moving isapprox to Test is inappropriate. It is a proper and time-tested way to establish convergence of an iteratively converging process.

[stevengj]

I just pushed a doc PR to clarify isapprox(x,0) and the default atol.

[mvhulten]

Now that I understand isapprox() much better than before, I believe that it is fine as it is. It is clean and can be useful even outside code testing, so in my opinion it can stay in Base.

The additional ?approx documentation of @stevengj clarifies the original issue of x ≈ 0.

[jlperla]

The new docs look great. Thanks for educating me.

[billtubbs]

@stevengj

NumPy got it wrong: their isclose function is only elementwise, has a nonzero default atol , and even has the property that isclose(x,y) is inequivalent to isclose(y,x) !

I might be biased but I find the Numpy isclose function intuitive. As many people have pointed out testing, is often an engineering task rather than a scientific task (if you need a sound test you should be writing a custom bound yourself).

np.isclose(0, 1e-10)  # True
np.isclose(10000000, 10000099)  # True
np.isclose(0, 0.0001)  # False
np.isclose(0, 0.0001, atol=0.001)  # True
np.isclose(1e-10, 2e-10, atol=0)  # False
np.isclose(10000000, 10000099, rtol=0, atol=1)  # False

I agree the default settings of rtol=1e-05, atol=1e-08 might be arbitrary but at least they reflect ‘common’ requirements.

[draftman9]

I think the result is decided by the ratio of subtraction of numbers where the equal signs are between and the number used to compare. Maybe following codes could help you:

julia> exp(-pi*im) + 1
0.0 - 1.2246467991473532e-16im

julia> imag(exp(-pi*im) + 1)
-1.2246467991473532e-16

julia> imag(exp(-pi*im) + 1) ≈ -1
false

julia> imag(exp(-pi*im) + 1) ≈ -1.22e-16
false

julia> imag(exp(-pi*im) + 1) +1 ≈ 1
true

julia> abs(exp(-pi*im) + 1) +1 ≈ 1
true