Add benchmarks (#738)

* Initial basic benchmarks

* Clean up a bit and add analysis

* WIP Make BigFloat interval faster

* Fix inf, sup, and abs

* Optimize sqrt for BigFloat

* Add the benchmark GitHub action from AirspeedVelocity.jl

* Clean up benchmarks

Author: Kolaru

.github/workflows/benchmark.yml

@@ -0,0 +1,14 @@
+name: Benchmark this PR
+on:
+  pull_request_target:
+    branches: [ master ]  # change to your default branch
+permissions:
+  pull-requests: write    # needed to post comments
+
+jobs:
+  bench:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: MilesCranmer/AirspeedVelocity.jl@action-v1
+        with:
+          julia-version: '1'

benchmark/Project.toml

@@ -0,0 +1,8 @@
+[deps]
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
+MPFI = "e50a78b8-6e2f-482e-b7c3-776aa2235d1a"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+
+[sources]
+MPFI = {url = "https://gitlab.inria.fr/ckatsama/mpfi.jl.git"}

benchmark/analyze.jl

@@ -0,0 +1,69 @@
+using BenchmarkTools
+using CairoMakie
+using DataFrames
+
+include("benchmarks.jl")
+
+results = run(SUITE ; verbose = true)
+
+df = DataFrame(; constructor = String[], suite = String[], f = String[], trial = BenchmarkTools.Trial[])
+
+for (name, T) in interval_constructors
+    # for suite in suites
+    begin
+        suite = "basics"
+        suite_df = DataFrame(results[name][suite], [:f, :trial])
+        suite_df[:, :constructor] .= name
+        suite_df[:, :suite] .= suite
+        df = vcat(df, suite_df)
+    end
+end
+
+transform!(df,
+    :trial => ByRow(trial -> minimum(trial.times)) => :minimum,
+    :trial => ByRow(trial -> mean(trial.times)) => :mean,
+    :trial => ByRow(trial -> median(trial.times)) => :median
+)
+
+df[:, :relative] .= 0.0
+
+for group in groupby(df, :f)
+    group[:, :relative] .= group[:, :median] ./ only(group[group.constructor .== "bareinterval", :median])
+end
+
+begin
+    fig = Figure(size = (800, 500))
+    
+    data = df[df.suite .== "basics", :]
+
+    fs = vcat(string.(basic_arithmetic), string.(basic_functions))
+    to_x = Dict(f => k for (k, f) in enumerate(fs))
+    xx = [to_x[f] for f in data.f]
+    
+    constructors = ["bareinterval", "interval", "BigFloat bareinterval", "BigFloat interval", "BigFloat MPFI"]
+    to_dodge = Dict(constructor => k for (k, constructor) in enumerate(constructors))
+    dodge = [to_dodge[constructor] for constructor in data.constructor]
+
+    ax = Axis(fig[1, 1] ;
+        xlabel = "Benchmarked function",
+        xticks = (1:length(fs), fs),
+        ylabel = "Relative execution time (log scale)",
+        yscale = log10,
+    )
+
+    barplot!(ax, xx, data.relative ;
+        dodge,
+        color = dodge,
+        colormap = :mpetroff_10,
+        colorrange = (1, 10),
+    )
+
+    Legend(fig[0, 1],
+        [LineElement(linewidth = 20, color = to_colormap(:mpetroff_10)[k]) for k in eachindex(constructors)],
+        constructors ;
+        tellwidth = false,
+        orientation = :horizontal
+    )
+
+    fig
+end

benchmark/benchmarks.jl

@@ -0,0 +1,49 @@
+using BenchmarkTools
+using IntervalArithmetic
+
+import MPFI
+import Random
+
+MPFI.BigInterval(x::ExactReal) = MPFI.BigInterval(x.value)
+Base.:(^)(x::MPFI.BigInterval, y::ExactReal) = x^y.value
+
+Random.seed!(0)
+
+basic_arithmetic = [+, *, -, /]
+basic_functions = [exp, cosh, sinh, tanh, inv, sqrt, abs, log, sin, cos, tan, acos, asin, atan]
+
+interval_constructors = Dict(
+    "bareinterval" => bareinterval,
+    "interval" => interval,
+    "BigFloat bareinterval" => (x, y) -> bareinterval(BigFloat, x, y),
+    "BigFloat interval" => (x, y) -> interval(BigFloat, x, y),
+    "BigFloat MPFI" => MPFI.BigInterval
+)
+
+suites = ["basics", "Dietmar-Ratz"]
+
+bounds = map(range(0.01, 10 ; length = 100)) do i
+    x = i * randn()
+    y = i * randn()
+    x > y && return (y, x)
+    return (x, y)
+end
+
+SUITE = BenchmarkGroup()
+
+for (name, T) in interval_constructors
+    xx = [T(x, y) for (x, y) in bounds]
+    yy = reverse(xx)
+
+    SUITE[name] = BenchmarkGroup(split(name))
+    SUITE[name]["basics"] = BenchmarkGroup(["arithmetic"])
+    SUITE[name]["Dietmar-Ratz"] = BenchmarkGroup(["arithmetic"])
+
+    for f in basic_functions
+        SUITE[name]["basics"][string(f)] = @benchmarkable ($f).($xx)
+    end
+
+    for op in basic_arithmetic
+        SUITE[name]["basics"][string(op)] = @benchmarkable ($op).($xx, $yy)
+    end
+end

src/intervals/arithmetic/absmax.jl

@@ -10,7 +10,9 @@ Implement the `abs` function of the IEEE Standard 1788-2015 (Table 9.1).
 """
 function Base.abs(x::BareInterval{T}) where {T<:NumTypes}
     isempty_interval(x) && return x
-    return _unsafe_bareinterval(T, mig(x), mag(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
 
 Base.abs(x::Interval) = _unsafe_interval(abs(bareinterval(x)), decoration(x), isguaranteed(x))

src/intervals/arithmetic/basic.jl

@@ -238,10 +238,9 @@ end
 Implement the `sqrt` function of the IEEE Standard 1788-2015 (Table 9.1).
 """
 function Base.sqrt(x::BareInterval{T}) where {T<:AbstractFloat}
-    domain = _unsafe_bareinterval(T, zero(T), typemax(T))
-    x = intersect_interval(x, domain)
-    isempty_interval(x) && return x
-    return @round(T, sqrt(inf(x)), sqrt(sup(x)))
+    sup(x) < 0 && return emptyinterval(BareInterval{T})
+    lo = inf(x) < 0 ? zero(T) : inf(x)
+    return @round(T, sqrt(lo), sqrt(sup(x)))
 end
 
 Base.sqrt(x::BareInterval{<:Rational}) = sqrt(float(x))

src/intervals/construction.jl

@@ -78,7 +78,8 @@ Internal constructor which assumes that `is_valid_interval(lo, hi) == true`.
 _unsafe_bareinterval(::Type{T}, a, b) where {T<:NumTypes} =
     _unsafe_bareinterval(T, _round(T, a, RoundDown), _round(T, b, RoundUp))
 
-_normalisezero(a) = ifelse(iszero(a), zero(a), a)
+_normalisezero(a) = ifelse(iszero(a), zero(a), a)  # Avoid branch in general
+_normalisezero(a::BigFloat) = iszero(a) ? zero(a) : a  # For BigFloat we avoid the allocation as much as possible
 # used only to construct intervals; needed to avoid `inf` and `sup` normalization
 _inf(x::BareInterval) = x.lo
 _sup(x::BareInterval) = x.hi

src/intervals/interval_operations/numeric.jl

@@ -14,7 +14,12 @@ Implement the `inf` function of the IEEE Standard 1788-2015 (Table 9.2).
 See also: [`sup`](@ref), [`bounds`](@ref), [`mid`](@ref), [`diam`](@ref),
 [`radius`](@ref) and [`midradius`](@ref).
 """
-inf(x::BareInterval{T}) where {T<:AbstractFloat} = ifelse(isnan(x.lo), typemax(T), ifelse(iszero(x.lo), copysign(x.lo, -1), x.lo))
+function inf(x::BareInterval{T}) where {T<:AbstractFloat}
+    isnan(x.lo) && return typemax(T)
+    iszero(x.lo) && return copysign(x.lo, -1)
+    return x.lo
+end
+
 inf(x::BareInterval{<:Rational}) = x.lo
 
 function inf(x::Interval{T}) where {T<:AbstractFloat}
@@ -38,7 +43,11 @@ Implement the `sup` function of the IEEE Standard 1788-2015 (Table 9.2).
 See also: [`inf`](@ref), [`bounds`](@ref), [`mid`](@ref), [`diam`](@ref),
 [`radius`](@ref) and [`midradius`](@ref).
 """
-sup(x::BareInterval{T}) where {T<:AbstractFloat} = ifelse(isnan(x.hi), typemin(T), x.hi)
+function sup(x::BareInterval{T}) where {T<:AbstractFloat}
+    isnan(x.hi) && return typemin(T)
+    return x.hi
+end
+
 sup(x::BareInterval{<:Rational}) = x.hi
 
 function sup(x::Interval{T}) where {T<:AbstractFloat}