Add benchmarks

[Kolaru]

This draft adds some basic benchmarks to the package.

MPFI.jl is benchmarked alongside the package.

There is a script that analyze and plot the results of the benchmark. I consider the minimum execution time of over a trial (recommend to reduce trial noise), and normalize it so that it is 1 for the Float64 BareInterval. For a selection of basic functions it gives the following:

It seems that the decoration are relatively not too expensive (except for very simple functions) and that we are not doing too bad compared to MPFI.jl (which doesn't seem to have support for Float64).

Currently there is only quite basics functions in the benchmark, let me know if you think something should be added.

PS: I do realize that the benchmark suite should not depend on Makie and DataFrames, I just need to think how to deal with that

[OlivierHnt]

This is great! The figures don't seem too bad indeed 🙂

[Kolaru]

I did some additional micro-benchmarks and noticed that our constructor for BareInterval{BigFloat} is slightly faster than MPFI.jl.

At least for abs the problem is our algorithm: because we focused on simplicity and readability, some checks are performed multiple times.

Currently we have (using abs as an example)

function Base.abs(x::BareInterval{T}) where {T<:NumTypes}
    isempty_interval(x) && return x
    return _unsafe_bareinterval(T, mig(x), mag(x))
end

Which looks fine, but mig and mag both check if the interval is empty again.

By changing it to

function Base.abs(x::BareInterval{T}) where {T<:NumTypes}
    isempty_interval(x) && return x
    inf(x) > 0 && return x
    sup(x) < 0 && return _unsafe_bareinterval(T, -sup(x), -inf(x))
    return _unsafe_bareinterval(T, zero(T), max(-inf(x), sup(x)))
end

I get between 1.5x and 4x speedup, depending on which branch is taken. In the simplest case (no-op), we even get faster than MPFI.

It feels a bit wrong, and it may be: while intervals are immutable, BigFloats are mutable in julia, so returning the input is technically not safe, as the value of the BigFloat may be changed from outside (using e.g. MutableArithmetics.jl). This is a concern more generally, every time we reuse a BigFloat.

Doing deep copies would again kill performances, so I believe that it is a tradeoff that we must take (and that maybe MPFI does as well).

Overall, my current conclusion is that we can probably make a lot of basic operations faster by hunting repeated checks and computations, at the possible cost of duplicating code.

Considering yesterday's discussion, I would assume that this is a step we would like to take @OlivierHnt @dpsanders @lbenet

[lbenet]

Very very interesting!

Would there be a gain if we @inline some "short" functions, e.g., isempty_interval?

I think I would agree to hunt down repeated checks and avoid them, but write somewhere the possible dangers when dealing with BigFloats, which I guess, we currently have.

[OlivierHnt]

With mig and mag we still have this problem with the mutability, right?

Indeed we were careful to have readable code, easier to debug and maintain over time. Not sure if @inbounds on isempty_interval will help (it may already be inlined) but it's worth a try.

We should probably first understand why (with BigFloat) +, -, *, / and inv are also a bit slower than MPFI before changing the code too much.

[Kolaru]

With mig and mag we still have this problem with the mutability, right?

Yes, it is everywhere. Even in the line

return _unsafe_bareinterval(T, zero(T), max(-inf(x), sup(x)))

There sup(x) return a reference to the upper bound, which is then passed to the newly constructed interval as is. Therefore if x.hi is mutated somewhere, it will affect the new interval.

We should probably first understand why (with BigFloat) +, -, *, / and inv are also a bit slower than MPFI before changing the code too much.

I agree, I will look into those next.

[Kolaru]

The next main offender was inf

inf(x::BareInterval{T}) where {T<:AbstractFloat} = ifelse(isnan(x.lo), typemax(T), ifelse(iszero(x.lo), copysign(x.lo, -1), x.lo))

By making it branching (and thus short-circuiting), we get much better performance for the simple operations:

We may very well be able to become faster than MPFI.jl if we keep digging a bit :D (unclear if we would beat native MPFI though, or intlab, so getting comparison with them is probably more important at this point).

[OlivierHnt]

Nice! I don't know if it's worth "digging more" at this stage. I think we can merge this PR, after a few tweaks.
Unless you are planning on adding more things before remove the "draft status".

[OlivierHnt]

Specialise to T<:BigFloat, and keep the branch free version

sup(x::BareInterval{T}) where {T<:AbstractFloat} = ifelse(isnan(x.hi), typemin(T), x.hi)

[Kolaru]

Are there reason to believe that in a larger code the branching would make a significant difference?

Looking at only microbenchmarks, both versions seem to be equivalent for Float64:

julia> s1(x) = ifelse(isnan(x.hi), typemin(x.hi), x.hi)
s1 (generic function with 1 method)

julia> function s2(x)
           isnan(x.hi) && return typemin(x.hi)
           return x.hi
       end
s2 (generic function with 1 method)

julia> @benchmark s1(X) setup=(X = bareinterval(rand(), rand() + 1))
BenchmarkTools.Trial: 10000 samples with 1000 evaluations per sample.
 Range (min … max):  2.654 ns … 24.137 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     2.805 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   2.822 ns ±  0.711 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  █    ▃                                                      
  █▁▁▁▁█▄▂▅▂▁▂▂▂▂▂▂▁▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂ ▂
  2.65 ns        Histogram: frequency by time         4.2 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

julia> @benchmark s2(X) setup=(X = bareinterval(rand(), rand() + 1))
BenchmarkTools.Trial: 10000 samples with 1000 evaluations per sample.
 Range (min … max):  2.669 ns … 22.299 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     2.806 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   2.829 ns ±  0.550 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▅    ▇█                                                     
  █▂▁▁▁██▂▂▂▁▁▁▁▂▂▂▂▁▁▂▁▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂ ▂
  2.67 ns        Histogram: frequency by time           4 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

[Kolaru]

That being said in this benchmark, all samples are taking the same branch, it may bias the benchmark.

[OlivierHnt]

Yes; though I am not sure I'd say "significant". From my understanding it is a "free" micro-optimisation (which the compiler might end up doing anyways) when the ifelse version bears no extra cost. This is the case except for BigFloat.

[OlivierHnt]

Specialise to T<:BigFloat and keep the branch free version

inf(x::BareInterval{T}) where {T<:AbstractFloat} = ifelse(isnan(x.lo), typemax(T), ifelse(iszero(x.lo), copysign(x.lo, -1), x.lo))

[OlivierHnt]

We still have a gap for sqrt. If there is an easy fix that's great, otherwise let's not get hung up on this and move forward.

EDIT: could the $y$-axis indicate the time unit (nano sec, micro sec)?

[OlivierHnt]

It's possible the extra cost in sqrt is due to the domain check:

domain = _unsafe_bareinterval(T, zero(T), typemax(T))
x = intersect_interval(x, domain)
isempty_interval(x) && return x

Can you check if changes these lines with

!issubset_interval(x, domain) && return emptyinterval(BareInterval{T})

fixes it?