New macro for guaranteed operations

[Kolaru]

I discussed with Luis and I put up a quick prototype on how we could extend the current @interval macro.

The strategy drafted here is to convert every literals into ExactNumber with the @exact macro.

An ExactNumber promote to real when mixed with reals, and to a guaranteed interval when interacting with intervals.

julia> x = ExactNumber(2.2)
ExactNumber{Float64}(2.2)

julia> x + 1
3.2

julia> x + interval(1, 1)
[3.19999, 3.20001]_com

The @exact macro works on both expression and function definition (thanks to MacroTools.jl)

julia> y = interval(1, 2)
[1.0, 2.0]_com

julia> @exact 1.2 * y + 3
[4.19999, 5.40001]_com

julia> @exact f(x) = 1.2 * x + 3
f (generic function with 1 method)

julia> f(y)
[4.19999, 5.40001]_com

julia> f(mid(y))
4.8

I haven't yet checked how it would play out with IntervalRootFinding or ForwardDiff.

[codecov-commenter]

Codecov Report

Attention: Patch coverage is 72.00000% with 7 lines in your changes are missing coverage. Please review.

Project coverage is 82.86%. Comparing base (c55dd23) to head (b5b584c).
Report is 1 commits behind head on master.

❗ Current head b5b584c differs from pull request most recent head d0842c3. Consider uploading reports for the commit d0842c3 to get more accurate results

Files	Patch %	Lines
src/intervals/exact_literals.jl	72.00%	7 Missing ⚠️

❗ Your organization needs to install the Codecov GitHub app to enable full functionality.

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #627      +/-   ##
==========================================
- Coverage   82.98%   82.86%   -0.13%     
==========================================
  Files          26       27       +1     
  Lines        2181     2206      +25     
==========================================
+ Hits         1810     1828      +18     
- Misses        371      378       +7     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[Kolaru]

ForwardDiff seems to not complain

julia> using IntervalArithmetic

julia> using ForwardDiff: derivative

julia> @exact f(x) = 3*sin(10x) - 2.2/x
f (generic function with 1 method)

julia> derivative(f, 1.2)
26.843396539752543

julia> derivative(f, interval(1.2, 1.2))
[26.8433, 26.8435]_com

Except when it uses comparisons of course

julia> @exact g(x) = x^2
g (generic function with 1 method)

julia> derivative(g, 1.2)
2.4

julia> derivative(g, interval(1.2, 1.2))
ERROR: ArgumentError: `==` is purposely not supported for intervals. See instead `isequal_interval`

[Kolaru]

Had a bit of fun and found a way to make the floating point input problem apparent:

julia> x = interval(0, 1)
[0.0, 1.0]_com

julia> @exact 0.5x + 0.1
┌ Warning: Float64 displayed as 0.1 is equal to 0.1000000000000000055511151231257827021181583404541015625
└ @ IntervalArithmetic ~/.julia/dev/IntervalArithmetic/src/intervals/exact_numbers.jl:33
[0.0999999, 0.600001]_com

The same trick should probably be used with @interval to warn the user.

[lbenet]

This is very very nice! 🚀

Let me begin with the following (side) remark: we can use the same trick to find the actual rounding (up or down) when converting a real number $x\in \mathbb{R}$ into a Float64. This could be used, e.g., to create a tight interval from the user provided $x$:

julia> x = 0.1
0.1

julia> ExactNumber(x)
ExactNumber{Float64}(0.1000000000000000055511151231257827021181583404541015625)

julia> string(x) < string(ExactNumber(x)) # so 0.1::Float64 is rounded up!
true

# A tight interval representation of the real 0.1 is:
julia> interval(prevfloat(x), x)
Interval{Float64}(0.09999999999999999, 0.1, com)

julia> x = 0.3
0.3

julia> ExactNumber(x)
ExactNumber{Float64}(0.299999999999999988897769753748434595763683319091796875)

julia> string(x) > string(ExactNumber(x)) # 0.3::Float64 is rounded down!
true

# A tight interval representation of the real 0.3 is:
julia> interval(x, nextfloat(x))
Interval{Float64}(0.3, 0.30000000000000004, com)

(I guess certain changes are required for negative numbers, but the idea is similar.)

The above raises the question: Do we want to "operate" with interval(0.1) or with the tight interval associated to $x = 0.1\in \mathbb{R}$?

[lbenet]

Aside from the comment above, and knowing this is still a prototype, what about other numeric formats (Rationals, BigFloats, ``Intervals? While (correctly) @exact` throws no warning for `Rational`s, there is a `MethodError` involving `//` (and also `/`) and `ExactNumber{Int64}`. Similar errors occur for `BigFloat's...

[Kolaru]

Let me begin with the following (side) remark: we can use the same trick to find the actual rounding (up or down)

Unfortunately, not quite, because we never know what the user did input. Anything that is parsed as 0.1 is equal for us.

julia> @exact 0.1
┌ Warning: Float64 displayed as 0.1 is equal to 0.1000000000000000055511151231257827021181583404541015625
└ @ IntervalArithmetic ~/.julia/dev/IntervalArithmetic/src/intervals/exact_numbers.jl:33
ExactNumber{Float64}(0.1000000000000000055511151231257827021181583404541015625)

julia> @exact 0.09999999999999999999782
┌ Warning: Float64 displayed as 0.1 is equal to 0.1000000000000000055511151231257827021181583404541015625
└ @ IntervalArithmetic ~/.julia/dev/IntervalArithmetic/src/intervals/exact_numbers.jl:33
ExactNumber{Float64}(0.1000000000000000055511151231257827021181583404541015625)

julia> @exact 0.1000000000000000055511151231257827021181583404541015625
┌ Warning: Float64 displayed as 0.1 is equal to 0.1000000000000000055511151231257827021181583404541015625
└ @ IntervalArithmetic ~/.julia/dev/IntervalArithmetic/src/intervals/exact_numbers.jl:33
ExactNumber{Float64}(0.1000000000000000055511151231257827021181583404541015625)

Aside from the comment above, and knowing this is still a prototype, what about other numeric formats

If we agree this is the way to go I'll add the required methods for Rational, BigFloat and the basic arithmetic operations.

[OlivierHnt]

So the purpose of the macro @exact is to notify that every literal numbers are "guaranteed", but not only, it would also make something like x= 1; @exact interval(1) + x be also guaranteed, correct?
I think this is nice.

EDIT: sorry I did not follow the discussion at the time. My understanding of the macro is that "whatever is typed in, is guaranteed". In this regard, we should not use an atomic interval for floats; for instance, @exact 0.1x should generate interval(Float64, 0.1)*x and not atomic(Float64, 0.1) * x.

[OlivierHnt]

Looking back at the current implementation, I find the goal of @exact a bit confusing.

I thought the purpose was to allow us to write general code (code used for both floating-points arithmetic and interval arithmetic), and avoid propagating the "NG" flag.
Of course, it is the responsibility of the user to know that their code complies with the specifications of both arithmetics.

To be concrete, suppose I have a function f(x) = 1 + 2x. Then this functions is completely valid for both floats and interval inputs. But currently the result is marked by "NG" if I use an interval; typically I would duplicate the code, and write a new method f(x::Interval) = @interval Float64 1+2x.

Now if I define g(x) = @exact 1+2x, I do not want to always get an interval output (as @interval would do). I would want g to also be used on floats (non-rigorous arithmetic). So g(1) == 3, while g(interval(1)) is interval(3) (with no "NG" flag).

This means that an operation involving ExactReal is only rigorous if all the other operands involved are intervals.
Thus, @exact sqrt(3) * x + 1 reads sqrt(ExactReal(3)) * x + ExactReal(1) with sqrt(ExactReal(3)) == sqrt(3).

I think this would be a very powerful feature. One situation I have in mind is rigorously computing a specific zero of a function f using Newton's method. Typically we want to perform as many iterations in floats, and only use intervals once we get close enough to the zero.

[Kolaru]

I think this PR is roughly where we want it to be.

The formulation of the doc would probably benefit from the review.

Once think that this PR is not doing is to deal with conversion of ExactReal{T} to Interval{S} when T and S are different (and T is not an Integer). In principle with different precision floating points number some care is required. I think it is niche enough to leave it unimplemented for now.

[OlivierHnt]

The CI failure seems to be due to broken tests introduced in #619, which now pass on nightly.

One tiny note: I am trying to think of a single word for the name of the macro. How about @guarantee, @literals? Or just @exact, since the purpose is to create an ExactReal?

[OlivierHnt]

@Kolaru
I made very minor changes to the PR. I reorganised the docstrings a bit.
I also removed the constructor ExactReal{T} since this is a fallback used everywhere for Real in Julia.

Lastly I ran into issues when the expression contained im, so I made an exception for this symbol in the macro. Can you double check this?

Also, I do not understand how the display works. Is the string being prompted truly the exact mathematical number?

[Kolaru]

One tiny note: I am trying to think of a single word for the name of the macro. How about @guarantee, @literals? Or just @exact, since the purpose is to create an ExactReal?

The most descriptive name would probably be @user_guaranteed_literals... I think the most reasonable is @exact.

Also, I do not understand how the display works. Is the string being prompted truly the exact mathematical number?

Yes. I ask the Ryu library (I have no idea what it is, but that's what Base uses) to write the number with 2000 decimals and truncate the extra 0s. According to my tests and estimation, 2000 decimals is more than enough to have the exact decimal representation of any Float64.

Lastly I ran into issues when the expression contained im, so I made an exception for this symbol in the macro. Can you double check this?

I fixed the constructor for complex such that @exact 1.2 + 3.4im now works. It may still be a bit brittle, if you have any additional example feel free to add them to the tests and I can have a look.

PS: @exact is a github user, so we need to be a bit careful to not randomly tag her ^^'

[OlivierHnt]

Looks good to merge, once the tests are successful. I glanced at it but I am not sure what went wrong.

[Kolaru]

I forgot to rename the macro.

[Kolaru]

I don't know what is going on with nightly or the doctests. Do the doctests use nightly julia ? It may be the same problem.

[OlivierHnt]

It seems that the problem with the doctests is due to the printing of the intervals. Somehow it expects setdisplay(:full).