Modify constructor to avoid Interval(Inf,Inf) or Interval(-Inf,-Inf) #51

lbenet · 2017-06-07T18:27:15Z

Solves #45.

dpsanders · 2017-06-07T22:46:25Z

I think Interval(Inf, Inf) should actually be an error and should be dealt with in the Interval constructor itself (or whatever replaces it for checks when constructing an interval).

coveralls · 2017-06-07T22:55:06Z

Coverage decreased (-0.5%) to 91.8% when pulling 5142ee8 on near_infinity into b54e103 on master.

lbenet · 2017-06-07T23:03:25Z

There are some examples in the ITF1788 test suite that go in a slightly different direction. One involves cancelminus (with an equivalent for cancelplus) and the other involves atanh. Both are already considered, though contrary to what you propose, they do not lead to an error.

The first one is equivalent to #45. In that case it is correct to have Interval(realmax(),Inf).

The second one corresponds to atanh(1.0, Inf). The domain for atanh is the interval [-1,1], and atanh(1.0) is Inf. So, by intersecting with the natural domain atanh([1.0,Inf]) becomes atanh([1.0,1.0]) which is the interval [Inf,Inf], hence the requirement that it yields the emptyinterval. This is precisely considered here.

dpsanders · 2017-06-07T23:06:35Z

Ah, I see: Interval(Inf, Inf) contains no real numbers, and is hence the empty interval.
In that case, that can also be done right in the constructor of Interval, and not by special casing those particular functions.

lbenet · 2017-06-07T23:10:07Z

The case of atanh is the only (so far) treated as a special case, since in that case the answer should not be Interval(prevfloat(Inf),Inf) as would be done by the constructor, but the empty interval.

lbenet · 2017-06-07T23:12:49Z

I just noticed some inconsistency; hold on a little bit...

lbenet · 2017-06-07T23:34:02Z

Actually, the problem of #45 is not the constructor, but rationalize, since rationalize(1e300) yields 1//0.

codecov-io · 2017-06-07T23:52:50Z

Codecov Report

Merging #51 into master will increase coverage by 0.03%.
The diff coverage is 100%.

@@            Coverage Diff             @@
##           master      #51      +/-   ##
==========================================
+ Coverage   91.64%   91.67%   +0.03%     
==========================================
  Files          21       21              
  Lines         945      961      +16     
==========================================
+ Hits          866      881      +15     
- Misses         79       80       +1

Impacted Files	Coverage Δ
src/intervals/intervals.jl	`95.65% <100%> (-4.35%)`	⬇️
src/intervals/conversion.jl	`61.53% <100%> (+5.01%)`	⬆️
src/intervals/hyperbolic.jl	`100% <100%> (ø)`	⬆️
src/intervals/functions.jl	`96.07% <100%> (+0.33%)`	⬆️
src/intervals/special.jl	`100% <100%> (ø)`	⬆️

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 2ecefe5...a176424. Read the comment docs.

lbenet · 2017-06-08T00:23:30Z

This PR has a problem: @interval(1e300) yields Interval(1.7976931348623157e308, Inf), which is obviously wrong.

The problem lies in this line of a convert method, since rationalize for large numbers returns 1//0, and thus the wrong result.

dpsanders · 2017-06-08T01:05:34Z

test/interval_tests/construction.jl

@@ -35,6 +35,10 @@ using Base.Test
    @test Interval{BigFloat}(1) == Interval{BigFloat}(big(1.0), big(1.0))
    @test Interval{BigFloat}(pi) == Interval{BigFloat}(big(pi), big(pi))

+    # a < Inf and b > -Inf
+    @test Interval(Inf, Inf) == Interval(prevfloat(Inf), Inf)


I don't understand this line.
Interval(Inf, Inf) should either return the empty interval or error.

It is an error in the octave interval package.

I agree it should probably yield an error or an empty interval; however, due to rationalize in the convert method which is used, we get Interval(Inf,Inf) when it shouldn't. So for the time being I think we should focus in getting around rationalize.

Three tests are skipped since they fail; the rest is passing.

lbenet · 2017-06-08T13:45:39Z

I just pushed another commit (after rebasing to current master).

Except for three of our tests, all tests pass! I am skipping checking those cases (a test comparing float intervals with the corresponding rationals, a test involving pi, and a test involving large values and tan) for the moment, so I can analyze carefully what is going on there.

dpsanders · 2017-06-08T13:55:58Z

src/intervals/conversion.jl

-    end
-    II
-end
+convert{T<:AbstractFloat}(::Type{Interval{T}}, x::Float64) =


This is just far too slow for such a basic operation that is used all over the place, e.g. in Interval(0.1) * 0.3:

julia> @btime Interval(0.1) * 0.3 700.447 ns (10 allocations: 416 bytes) Interval(0.03, 0.03000000000000001) julia> @btime Interval(0.1) 7.612 ns (0 allocations: 0 bytes) Interval(0.1, 0.1) julia> @btime convert(Interval{Float64}, 0.3) 610.783 ns (10 allocations: 416 bytes) Interval(0.3, 0.30000000000000004) julia> @btime Interval(prevfloat(0.3), nextfloat(0.3)) 12.497 ns (0 allocations: 0 bytes)

(A factor of 50 difference!)

I suggest we use parse for something like convert(Interval, "0.3") (when a string is given), and simply use convert(::Type{Interval{Float64}, x) = Interval(prevfloat(x), nextfloat(x)) for floats.

The only thing I did is to eliminate rationalize and the conditional, so that method directly converts to a string which is parsed and rounded.

Unfortunately the conversion to and from a string is very slow, that's what I'm pointing out. (And probably rationalize was too.)

Hmm, it seems maybe I'm wrong and this is not the function that is actually used?

No, I'm right.

Rationalize has two problems: sometimes it yields a zero for too small floating points, and sometimes it yields Inf for too large values. I'll see if I can improve the speed...

lbenet · 2017-06-08T14:24:30Z

I think there are two issues in this discussion we should separate: the construction of an interval (which contains whatever checks we wish), and the use of convert.

As for the construction, i.e., Interval(a,b) there is essentially no difference.

As you write, we should improve the convert method, which is used whenever we round, i.e., always. Your suggestion is not viable because it does not respect the tightest character that the interval must have.

dpsanders · 2017-06-08T14:34:53Z

Well, part of the problem as usual is that it is not clear what one means when writing 0.3 -- do you mean the decimal number "0.3", or the floating point number that is printed as 0.3.

As far as I can tell, the standard only talks about what to do with strings, in which case certainly we must give the tightest interval, however long that takes: that is, if the user uses a string, it is a sign that (s)he wants the tightest possible interval.

With floats there is, as far as I can tell, basically no way of making it performant to get the tightest interval if you interpret 0.3 as "the exact decimal number 0.3". If you interpret 0.3 as "some
real number that is closest to the floating-point number 0.3" (i.e. an uncertain unum-like object) then it actually represents a small interval of width eps(0.3)/2 (roughly), and so it makes total sense to do outward rounding with prevfloat and nextfloat.

lbenet · 2017-06-08T15:02:28Z

I disagree: the purpose of the standard is to set-up definitions and rules to get the same; returning a tightest interval is part of it.

Quoting from the standard (italics are mine):

"The aim of the standard is to improve the availability of reliable computing in modern hardware and software environments by defining the basic building blocks needed for performing interval arithmetic. There are presently many systems for interval arithmetic in use; lack of a standard inhibits development, portability, and ability to verify correctness of codes."

" This standard specifies:
...
Constructors for intervals from numeric and character sequence data"

coveralls · 2017-06-08T15:22:26Z

Coverage increased (+0.009%) to 91.649% when pulling 772cac4 on near_infinity into 2ecefe5 on master.

dpsanders · 2017-06-08T16:14:27Z

The relevant section of the standard is 10.5.8 (pg. 38): there are only two required constructors, namely (in our notation):

Interval(a, b) for floats a and b [this is numsToInterval in the Standard]
textToInterval(s), for a text string s

(We certainly need to work on the second one for it to be compliant with the standard.)

dpsanders · 2017-06-08T16:25:22Z

Here's another relevant piece (pg. 30):
"An implementation may support an extended form of literals, e.g., using number literals in the syntax of the host language of the implementation. It may restrict the support of literals at Level 2, by relaxing conversion accuracy of hard cases: rational number literals, long strings, etc. What extensions and restrictions of this kind are permitted is flavor-defined."

(Italics are mine.)

lbenet · 2017-06-08T16:28:58Z

All clause 10 is related to Level 1 description, i.e., to the mathematical theory; the level 2 description (finite precision) corresponds to the clause 12. The measures of accuracy is treated in section 12.10. The mode we are discussing is called accurate and is the intermediate one with respect to the tightness.

lbenet · 2017-06-08T16:29:30Z

Perhaps we should use traits also for the constructors...

dpsanders · 2017-06-08T16:32:03Z

Yes, that is certainly a possibility. What did you have in mind? Something like

Interval(:tight, 0.3)

vs

Interval(0.3)

which may not be tight, could work.

Or

Interval(Tight(), 0.3)

to use dispatch.

dpsanders · 2017-06-08T16:38:53Z

Also pg. 60:
"An implementation may provide more than one version of some operations for a given type. For instance, it may provide an “accurate” version in addition to a required “tightest” one, to offer a trade-off of accuracy versus speed or code size. How such a facility is provided is language- or implementation-defined."

lbenet · 2017-06-08T16:51:01Z

Regarding traits, I was thinking in something like Interval{T, RM}, where RM is a rounding mode (e.g., :tight), but you know all these things better than I do.

Regarding the last quotation you provided, the key part for me is "... in addition to a required “tightest” one,...". That's the reason I think we should focus on this one, and the allow for others.

dpsanders · 2017-06-08T22:40:07Z

Shall we merge this for now, @lbenet?

dpsanders · 2017-06-09T14:05:21Z

I'm not sure why the result for tan changed (do you understand why?), but I see no problem with returning [-Inf, Inf] for x = Interval(2)^1000. (But apparently that causes other issues?)

dpsanders · 2017-06-09T14:07:51Z

Could you give the code sequence that gives Interval(Inf, Inf)? I don't understand where that is.

lbenet · 2017-06-09T14:41:28Z

I'm not sure why the result for tan changed (do you understand why?), but I see no problem with returning [-Inf, Inf] for x = Interval(2)^1000. (But apparently that causes other issues?)

I agree that [-Inf,Inf] contains the solution, but in master the test case returns a tighter interval which contains the solution. That is what I wanted to reproduce.

Could you give the code sequence that gives Interval(Inf, Inf)? I don't understand where that is.

This is my understanding: find_quadrants creates the interval temp and rounds to the previous (smaller) integer the bounds of the interval. In master, the line temp = x / half_pi(x) returns Interval(Inf,Inf), which is wrong(!). Then, lo_quadrant_mod and hi_quadrant_mod which are Inf succeed skipping all if's and then the returned interval is evaluated as @round(tan(a.lo), tan(a.hi)), which happens to be a tighter interval (than the entire one), seemingly containing the true answer.

Yet, this PR is excluding precisely intervals of the form Interval(Inf,Inf). Then, temp becomes Interval(6.821435655968254e300, 6.821435655968257e300), which is such a wide interval that the interval returned is the entire interval.

The current result [-Inf,Inf] is correct, but maybe we can do better...

dpsanders · 2017-06-09T14:44:00Z

The line with temp only ever acts on floats, not on intervals?

lbenet · 2017-06-09T14:44:20Z

One of the tests that is not passing (due to this PR) is this one. This, I think, is related to the changes in convert, and more specifically, to the fact that it doesn't use anymore rationalize.

What troubles me is the following:

julia> x = 1/3
0.3333333333333333

julia> @interval(x)
Interval(0.33333333333333326, 0.3333333333333333)

julia> big(1)/3 ∈ ans  # AUCH!!!
false

On the other hand:

julia> big(1)/3
3.333333333333333333333333333333333333333333333333333333333333333333333333333348e-01

julia> @interval(1/3)
Interval(0.3333333333333333, 0.33333333333333337)

julia> big(1)/3 ∈ ans
true

So rationalize was somehow helping. Any idea?

To go back to rationalize means to have another conditional when rationalize yields 1//0 which still can be handled as a Float64 number, aside from the performance issue.

lbenet · 2017-06-09T14:45:53Z

The line with temp only ever acts on floats, not on intervals?

I agree, but it creates an interval; in master we get:

julia> x = 2.0^1000
1.0715086071862673e301

julia> x/IntervalArithmetic.half_pi(x)
Interval(Inf, Inf)

dpsanders · 2017-06-09T15:21:09Z

Not at my computer currently but I thought I tried that and did not get that result. That would definitely be a bug, in division apparently?

lbenet · 2017-06-09T15:24:29Z

i haven't checked the division, but I think the problem is related to #45; somewhere we must rely on convert(Interval{Float64}, x) with x a large enough number.

dpsanders · 2017-06-09T15:27:56Z

Ah yes, good catch: since x is divided by an Interval, the first thing that happens is that it gets converted to an Interval.

dpsanders

Interval(Inf,Inf) is defined to be an error by the standard. Apart from that, this looks good, thanks.

dpsanders · 2017-06-09T15:29:38Z

src/intervals/conversion.jl

-                        parse(T, string(x), RoundUp))
-    end
-    II
+    isinf(x) && return Interval{T}(x)


If x is infinite, then it should be an error to convert it to an Interval (though this can happen in the constructor itself).

Good point. Currently this was returning one of the edge intervals, e.g. Interval(realmax(), Inf), which was handled at the constructor. The idea to include that here was to skip the rest of the method. (I will change the constructor as you suggest below, and the adapt this.)

dpsanders · 2017-06-09T15:32:01Z

src/intervals/intervals.jl

+        # The IEEE Std 1788-2015 states that a < Inf and b > -Inf;
+        # see page 6, "bounds".
+        if a == Inf
+            return new(prevfloat(a), b)


a == Inf should give an error.

There is a case (cancelminus I believe) where the subtlety is important.

In any case, I think you are right. I will change this to throw an error (ArgumentError?) and treat the special case separately.

lbenet · 2017-06-09T15:46:22Z

@dpsanders Thanks for the review. I will incorporate the changes. Should I also clean-up (cherry-pick and push -f) this PR?

… method Currently `@interval(Inf)` returns `Interval(realmax(),Inf)`

dpsanders · 2017-06-09T21:30:17Z

Sure, do whatever needs doing. Is rebase not enough?

lbenet · 2017-06-09T22:55:56Z

Pushing what I have now; it's messy (with respect to the end result) but I think it works fine. I want to let travis do its job.

coveralls · 2017-06-09T23:38:57Z

Coverage increased (+0.04%) to 91.675% when pulling a176424 on near_infinity into 2ecefe5 on master.

lbenet · 2017-06-10T00:25:27Z

Tests are passing, including few additions :)

So this is a summary of what I did:

I modified the constructor as we discussed, so Interval(Inf,Inf) throws an error;
rationalize is back, so h=1/3; @interval(h * h) yields an interval containing big(1)/9;
I defined the (unexported) function wideinterval(x) which returns the interval constructed with prevfloat(x) and nextfloat(x);
I included conditionals of the form (a.lo == Inf || a.hi ==-Inf) && return wideinterval(a.lo), which is sometimes needed (e.g. ^).

The following is a feature: @interval(Inf) == wideinterval(Inf).

lbenet · 2017-06-10T00:27:50Z

Incidentally, @round(a, b) may lead to headaches in case a and b may become both Inf...

lbenet · 2017-06-10T00:44:23Z

Only the test introduced in #22 is skipped; the result of tan( Interval(2.0)^1000 ) is the entire interval...

dpsanders · 2017-06-10T14:38:30Z

This is looking very nice, thanks!

What was the problem with the power with infinite limits?

In my opinion this is ready to merge. If you agree, then please go ahead!

lbenet · 2017-06-10T15:33:44Z

The problem is that, sometimes, @round(a,b) may lead to Interval(Inf,Inf); this was throwing an error (triggered by the constructor) despite that it could be to a computable interval. That's the reason I had to split the operations and use @round_down and @round_up. (I tried to do it directly in the macros, but was unable.) We should keep this in mind.

Another comment: as far as I can see, rationalize is really required by @interval, but seemingly not when using @round. The only tests that were failing without it were things like @interval( h*h ) for h=1/3. Maybe we can skip using it in the convert method where it is now by adding a specialize convert method which @interval will use. What do you think about this?

Aside from the last comment, which may be left for another PR, I agree that this should be merged.

lbenet · 2017-06-11T02:13:24Z

Merging!

dpsanders · 2017-06-11T02:32:30Z

Great, thanks!

…

On 10 Jun 2017 9:13 p.m., "Luis Benet" ***@***.***> wrote: Merged #51 <#51>. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#51 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AALtTrVCFnxWJC-BzVBNcnqKBlpC9E6Vks5sC01KgaJpZM4NzGqE> .

dpsanders mentioned this pull request Jun 7, 2017

Make Interval(Inf, Inf) an error #52

Closed

lbenet mentioned this pull request Jun 8, 2017

Avoid use of rationalize in conversion? #53

Closed

dpsanders reviewed Jun 8, 2017

View reviewed changes

Luis Benet added 2 commits June 8, 2017 08:36

Modify constructor to avoid Interval(Inf,Inf) or Interval(-Inf,-Inf)

ee01bb4

Remove rationalize from convert method

772cac4

Three tests are skipped since they fail; the rest is passing.

lbenet force-pushed the near_infinity branch from 5142ee8 to 772cac4 Compare June 8, 2017 13:41

dpsanders reviewed Jun 8, 2017

View reviewed changes

Revert cancelminus to a former version (tests were not passing)

e237a2c

dpsanders requested changes Jun 9, 2017

View reviewed changes

Throw an error for Interval(Inf,Inf) and reuse rationalize in convert…

09ceda6

… method Currently `@interval(Inf)` returns `Interval(realmax(),Inf)`

Add wideinterval, retouch convert and fix ^

a176424

dpsanders approved these changes Jun 10, 2017

View reviewed changes

lbenet merged commit 2ed8164 into master Jun 11, 2017

lbenet mentioned this pull request Jun 11, 2017

Bug creating intervals from large floats #45

Closed

This was referenced Jun 11, 2017

Bug in tan #21

Closed

convert(Interval, x) is wrong for large x JuliaIntervals/ValidatedNumerics.jl#259

Closed

dpsanders deleted the near_infinity branch August 2, 2017 19:40

cbilz mentioned this pull request Sep 28, 2022

convert can fail to return a valid interval #539

Closed

Modify constructor to avoid Interval(Inf,Inf) or Interval(-Inf,-Inf) #51

Modify constructor to avoid Interval(Inf,Inf) or Interval(-Inf,-Inf) #51

Conversation

lbenet commented Jun 7, 2017

dpsanders commented Jun 7, 2017

coveralls commented Jun 7, 2017 • edited Loading

lbenet commented Jun 7, 2017

dpsanders commented Jun 7, 2017

lbenet commented Jun 7, 2017

lbenet commented Jun 7, 2017

lbenet commented Jun 7, 2017

codecov-io commented Jun 7, 2017 • edited Loading

Codecov Report

lbenet commented Jun 8, 2017

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

lbenet commented Jun 8, 2017

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

lbenet commented Jun 8, 2017

dpsanders commented Jun 8, 2017

lbenet commented Jun 8, 2017

coveralls commented Jun 8, 2017 • edited Loading

dpsanders commented Jun 8, 2017 • edited Loading

dpsanders commented Jun 8, 2017

lbenet commented Jun 8, 2017

lbenet commented Jun 8, 2017

dpsanders commented Jun 8, 2017

dpsanders commented Jun 8, 2017

lbenet commented Jun 8, 2017

dpsanders commented Jun 8, 2017

dpsanders commented Jun 9, 2017

dpsanders commented Jun 9, 2017

lbenet commented Jun 9, 2017

dpsanders commented Jun 9, 2017

lbenet commented Jun 9, 2017

lbenet commented Jun 9, 2017

dpsanders commented Jun 9, 2017

lbenet commented Jun 9, 2017

dpsanders commented Jun 9, 2017

dpsanders left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

lbenet commented Jun 9, 2017

dpsanders commented Jun 9, 2017 • edited Loading

lbenet commented Jun 9, 2017

coveralls commented Jun 9, 2017 • edited Loading

lbenet commented Jun 10, 2017

lbenet commented Jun 10, 2017

lbenet commented Jun 10, 2017

dpsanders commented Jun 10, 2017

lbenet commented Jun 10, 2017

lbenet commented Jun 11, 2017

dpsanders commented Jun 11, 2017 via email

coveralls commented Jun 7, 2017 •

edited

Loading

codecov-io commented Jun 7, 2017 •

edited

Loading

coveralls commented Jun 8, 2017 •

edited

Loading

dpsanders commented Jun 8, 2017 •

edited

Loading

dpsanders commented Jun 9, 2017 •

edited

Loading

coveralls commented Jun 9, 2017 •

edited

Loading