Constructing intervals

[dpsanders]

I suggest that we should have the following mechanisms for constructing intervals:

• interval(a, b) (with a small i): the current Interval(a, b) -- includes the checks at the moment of construction; constructs an interval from exactly the given value of a to the given value of b; fast

• a..b: Gives a guaranteed enclosure of the interval, including both a and b; see Speed up a..b #34; fast.
  Does not try to do anything clever with floats; just uses prevfloat and nextfloat

• convert(Interval, a): Gives the smallest possible enclosure of the number a as an interval; does fancy things (rationalize or parse(string(x))) to interpret e.g. 0.1 as the real number 0.1; slow (same as currently)

• @interval(a, b): The most general method, but slowest.

• a ± b: Only for numbers a and b (not intervals)

cc @lbenet

[dpsanders]

Thinking about it, I think convert(Interval{Float64}, x::Float64) also needs to be fast, e.g. for

julia> x = 1..2
[1, 2]

julia> x * 0.3

Here, 0.3 gets converted to the type of x in the promotion.

[dpsanders]

However, we still want some way of making the smallest interval around 0.3.
This is currently @interval(0.3); however, that uses convert(Interval, 0.3) internally.

[lbenet]

Thanks to bring this up again; perfect timing!

The idea of this comment is to exploit dispatch rather than conversion. I think that such an approach may be better and "more julian".

Just to recall it, we should have two types, one which is "safe" (in the sense of checks and correct rounding), and the other that is "fast". Both should be subtypes of AbstractInterval{T}. Naming such types could be Interval or SafeInterval, and FastInterval or UnsafeInterval, or whatever you prefer. In terms of constructors, I think we would essentially need one internal constructor for SafeInterval which checks the correct order of the bounds, and could also handle rounding.

From what you propose above, I think I like the idea of using different short-forms for each type, for example, a..b could build a FastInterval.

[dpsanders]

Refactor of .. in PR #37

[lbenet]

I'm playing with the idea above, or having a fast constructor for intervals (FastInterval). Do you have some benchmark case?

[dpsanders]

On master:

julia> using BenchmarkTools

julia> @btime Interval(1, 2)
  7.994 ns (0 allocations: 0 bytes)

On simplify_constructor branch:

julia> @btime Interval(1, 2)  # without checks
  1.393 ns (0 allocations: 0 bytes)

julia> @btime interval(1, 2)   # with checks
  8.908 ns (0 allocations: 0 bytes)

We should actually make a proper benchmark suite.

[dpsanders]

(Edited above to add interval)

[dpsanders]

Note that any benchmark that includes arithmetic will currently be dominated by those arithmetic operations, since they are currently really slow! That is what the fast_rounding branch solves (which is included in simplify_constructor, possibly unfortunately).

[dpsanders]

Having said that, it is no longer clear to me that removing the checks from the constructor really has that much influence on the overall time for a real calculation, which will indeed be dominated by other costs.

[dpsanders]

For example,

julia> f1(N) = sum(Interval(1) / (i^2) for i in 1:N)
f1 (generic function with 1 method)

julia> f2(N) = sum(interval(1) / (i^2) for i in 1:N)
f2 (generic function with 1 method)

julia> f1(10)
[1.54976, 1.54977]

julia> f2(10)
[1.54976, 1.54977]

julia> @btime f1(100)
  3.627 μs (3 allocations: 96 bytes)
[1.63498, 1.63499]

julia> @btime f2(100)
  4.399 μs (3 allocations: 96 bytes)
[1.63498, 1.63499]

There is a 20% increase in time due to the extra checks in the constructor.
(Note that the constructor is called, for example, to create the result of the division.)

[lbenet]

This is what I get:

julia> @btime Interval(pi, 7pi)
  11.307 ns (0 allocations: 0 bytes)
[3.14159, 21.9912]

julia> @btime @interval(pi, 7pi)
  1.122 μs (7 allocations: 192 bytes)
[3.14159, 21.9912]

julia> @btime FastInterval(pi, 7pi)
  8.488 ns (0 allocations: 0 bytes)
IntervalArithmetic.FastInterval{Float64}(3.141592653589793, 21.991148575128552)

FastInterval does not check the proper bounds of the interval, so it is the simplest possible constructor, Interval does check the bounds but performs no rounding, and @interval does both. This is consistent with your last remark, that checking the bounds does influence the timing, but not too crucially, while the rounding business is what costs allocations and time.

Let me advance further and then I push it to GitHub.

[lbenet]

Among others, I still have to adapt the arithmetic definitions...

[dpsanders]

There is some time being spent in 7*pi. Try with just numbers.

The fact that you have to adapt the arithmetic definitions is one reason why this doesn't seem to me to be an attractive direction to take!

[lbenet]

The pi and 7*pi was on purpose. The following is another example:

julia> @btime Interval(1, 2)
  11.307 ns (0 allocations: 0 bytes)
[1, 2]

julia> @btime FastInterval(1, 2)
  8.487 ns (0 allocations: 0 bytes)
IntervalArithmetic.FastInterval{Float64}(1.0, 2.0)

julia> @btime @interval(1, 2)
  756.105 ns (5 allocations: 128 bytes)
[1, 2]

[lbenet]

Regarding adapting the operations, the idea is to exploit that both Interval and FastInterval are subtypes of AbstractInterval.

[lbenet]

Not related to FastInterval, but still interesting:

julia> @btime convert(Interval{Float64}, 0.1)
  303.565 ns (3 allocations: 128 bytes)
[0.0999999, 0.100001]

julia> @btime @interval 0.1
  675.000 ns (5 allocations: 192 bytes)
[0.0999999, 0.100001]

[dpsanders]

Yes, that is interesting -- it seems to be due to this:

julia> @macroexpand @interval 0.1
:((IntervalArithmetic.Interval)((IntervalArithmetic.convert)((IntervalArithmetic.Interval){IntervalArithmetic.parameters.precision_type}, 0.1)))

i.e. it must do a dynamic dispatch (i.e. at runtime) to convert to the correct type.

[dpsanders]

Note the following (with simplify_constructor), that seems to be due to instruction pipelining:

julia> using AdjacentFloats

julia> @btime Interval(prev_float(1.0), next_float(2.0))
  1.386 ns (0 allocations: 0 bytes)
[0.999999, 2.00001]

[lbenet]

I just pushed the fastinterval branch, where the idea of having Interval or FastInterval is basically implemented. I didn't implement any tests yet. Both (as well as DecoratedInterval) are subtypes of AbstractInterval, which is used whenever possible to extend the methods.

Here are some benchmarks using Julia v0.6.0-rc2, as we outlined above:

julia> using IntervalArithmetic, BenchmarkTools

julia> @btime Interval(1, 2)
  11.307 ns (0 allocations: 0 bytes)
[1, 2]

julia> @btime FastInterval(1, 2) # no checks
  8.173 ns (0 allocations: 0 bytes)
IntervalArithmetic.FastInterval{Float64}(1.0, 2.0)

julia> f1(N) = sum(Interval(1) / (i^2) for i in 1:N)
f1 (generic function with 1 method)

julia> f2(N) = sum(FastInterval(1) / (i^2) for i in 1:N)
f2 (generic function with 1 method)

julia> setformat(:full)
6

julia> f1(10)
Interval(1.5497677311665397, 1.549767731166541)

julia> f2(10)
IntervalArithmetic.FastInterval{Float64}(1.5497677311665408, 1.5497677311665408)

julia> @btime f1(100) # uses Interval; slow
  17.795 μs (3 allocations: 96 bytes)
Interval(1.634983900184882, 1.6349839001849027)

julia> @btime f2(100) # uses FastInterval
  2.876 μs (3 allocations: 96 bytes)
IntervalArithmetic.FastInterval{Float64}(1.6349839001848923, 1.6349839001848923)

[dpsanders]

I was just about to write that I didn't understand why FastInterval was so much faster, but it is not doing any rounding?

[dpsanders]

The amazing thing about the fastrounding branch is that doing correct rounding is almost as fast as no rounding at all. There is one bug left to fix before it is ready.

[lbenet]

You are right, it is doing no rounding. The subtlety is related to @round that returns an Interval. I'll try to implement it in such a way that it rounds correctly.

[dpsanders]

You could change @round so that it just returns the tuple of rounded values, and then wrap that in the correct type inside each function. Painful... I really don't think that this is the right way to go.

[lbenet]

I thought that changing @round(expr1, expr2) to @round(T, expr1, expr2), with T being either Interval or FastInterval would work, but I am having problems...

[dpsanders]

I rewrote @round to get rid of putting explicit types in...

[lbenet]

I don't like my solution to include rounding, but it's up there. Currently, I get the following

julia> using IntervalArithmetic, BenchmarkTools

julia> @btime Interval(1, 2)
  11.307 ns (0 allocations: 0 bytes)
[1, 2]

julia> @btime FastInterval(1, 2)
  8.173 ns (0 allocations: 0 bytes)
IntervalArithmetic.FastInterval{Float64}(1.0, 2.0)

julia> f1(N) = sum(Interval(1) / (i^2) for i in 1:N)
f1 (generic function with 1 method)

julia> f2(N) = sum(FastInterval(1) / (i^2) for i in 1:N)
f2 (generic function with 1 method)

julia> setformat(:full)
6

julia> f1(10), f2(10)
(Interval(1.5497677311665397, 1.549767731166541), IntervalArithmetic.FastInterval{Float64}(1.5497677311665397, 1.549767731166541))

julia> @btime f1(100) # uses Interval; slow
  17.632 μs (3 allocations: 96 bytes)
Interval(1.634983900184882, 1.6349839001849027)

julia> @btime f2(100) # uses FastInterval
  16.094 μs (3 allocations: 96 bytes)
IntervalArithmetic.FastInterval{Float64}(1.634983900184882, 1.6349839001849027)

From this, it seems to me that the whole issue is not so much to check the bounds, but to round faster. And currently we don't do it.

[dpsanders]

Yes, I completely agree -- this is all pointless with the current slooow rounding.
That's why it's urgent to get the fastrounding branch ready!

(Actually that branch is not the problem; there are two remaining bugs to fix in the packages that it will depend on.)

[lbenet]

I partially agree with you: I think there are two points here:

• Are the checks (in Interval) relevant? (The point of this issue.) The answer is "not quite"
• Make an efficient rounding scheme (independent of checking or not the end-points).

I'll take a look on #25.

[dpsanders]

I'm not sure what you mean by "not quite". (Or by "revelant".)

[lbenet]

What I mean by "not quite" is that using Interval (that checks the correctness of the bounds and rounds) and FastInterval (that does not check the bounds but rounds) yields timings that differ in ~10%. On the other hand, the differences of the naive benchmarks when FastInterval includes rounding and when it doesn't yields a factor 8x in speed-up.

[dpsanders]

You can try mixing your branch and fast_rounding. Something like

git checkout fast_rounding
git checkout -b fast_rounding_fast_interval
git merge fast_interval

[lbenet]

While I managed to merge both branches (only needed to correct two conflicts) I am having problems with FastRounding (actually on Nemo)...

[dpsanders]

I'm not sure what you mean by "on nemo". I guess youve done Pkg.add("FastRounding")?

[lbenet]

I'm not sure what was I doing wrong...

Anyway, here are the new benchmarks:

julia> using IntervalArithmetic, BenchmarkTools

julia> setrounding(Interval, :errorfree)

julia> @btime Interval(1, 2)
  10.992 ns (0 allocations: 0 bytes)
[1, 2]

julia> @btime FastInterval(1, 2)
  8.173 ns (0 allocations: 0 bytes)
IntervalArithmetic.FastInterval{Float64}(1.0, 2.0)

julia> f1(N) = sum(Interval(1) / (i^2) for i in 1:N)
f1 (generic function with 1 method)

julia> f2(N) = sum(FastInterval(1) / (i^2) for i in 1:N)
f2 (generic function with 1 method)

julia> setformat(:full)
6

julia> f1(10), f2(10)
(Interval(1.5497677311665397, 1.549767731166541), IntervalArithmetic.FastInterval{Float64}(1.5497677311665397, 1.549767731166541))

julia> @btime f1(100)
  8.307 μs (3 allocations: 96 bytes)
Interval(1.634983900184882, 1.6349839001849027)

julia> @btime f2(100)
  5.963 μs (3 allocations: 96 bytes)
IntervalArithmetic.FastInterval{Float64}(1.634983900184882, 1.6349839001849027)

So, f1 is 2.1x faster, and f2 is 2.7x faster. Pretty impressive.

All our tests pass; ITF1788 tests also pass, except for two tests related with cancelPlus and cancelMinus. I'll check this (and push to github the fast_rounding_fastinterval branch.

[dpsanders]

That's a bit disappointing -- I was expecting a better speedup.

Those two tests are from a (very) subtle bug that is on my radar: JeffreySarnoff/FastRounding.jl#11

[dpsanders]

Although the code for cancelplus and cancelminus seems to be surprisingly complicated.
(And I still don't really understand the need for those two functions, to be honest.)

[lbenet]

These function provide something like the inverse for addition and subtraction; I included them because they are in the standard. The complicated programing was the way I found to get the tests passing...

[lbenet]

I just pushed another commit where cancelminus is rewritten. This current version does not plays with the precision, and seems to work for both setrounding(Interval, :correct) and setrounding(Interval, :errorfree). All tests should pass.

[lbenet]

I can make a new PR with the sole correction of cancelminus, which I think is worth to have in master.

[dpsanders]

That would be great, yes! And cancelplus is in terms of cancelminus?

[lbenet]

Yes: cancelplus(a,b) = cancelminus(a,-b).

I will prepare a separate PR to handle this.

[dpsanders]

Is there anything left to do here @lbenet ?

[lbenet]

I must admit that I completely forgot about this discussion. I think the discussion is quite outdated and a lot has changed since we were talking about this: For instance, we have now atomic, and convert behaves quite differently.

Maybe your question was related to cancelminus and cancelplus. If so, I can look up again into the fastinterval branch and try to rescue the relevant PRs related to those functions.

[dpsanders]

The current behaviour is now roughly what was sketched at the start of the thread and is documented in the docs.