IntervalArithmetic benchmark

mwe2_hull.jl

# MWE 2 — `hull`: the widest per-operation gap between the two backends
# (450 ns vs 9.0 us for a 999-step reduction over 1000 intervals, 20x).
#
# Both backends load in one session; no TaylorRemainders preference involved:
#
#   julia --project=env_ifm mwe2_hull.jl
#
# Each variant is timed two ways, which separate two distinct effects:
#   * scalar    — one `hull` call on loop-invariant arguments: per-call cost.
#   * reduce    — `reduce(hull, v)` over 1000 intervals: per-call cost *plus*
#                 whatever the compiler can or cannot vectorize.
# A variant that is only modestly slower per call but far slower under `reduce`
# is losing to vectorization, not to arithmetic.
#
# IA's `hull(::Interval, ::Interval)` layers three things on top of the bare
# min/max: NaI checks, decoration propagation (`min(decoration(x),
# decoration(y))`), and a `dec` keyword routed through `_set_decoration` ->
# `setdecoration`, which re-examines the result for emptiness and
# unboundedness. `hull_manual` performs the same computation with the
# decoration fixed at the value the `:default` keyword selects.
#
# Profiling entry points: loop_hull_ifm(v), loop_hull_ia(v),
# loop_hull_bare(v), loop_hull_manual(v), each taking a vector.

using BenchmarkTools
import IntervalFastMath as IFM
import IntervalArithmetic as IA
using ThickNumbers: hull
using Random

const IFMI = IFM.Interval{Float64}
const IAI = IA.Interval{Float64}

# --- hull variants -----------------------------------------------------------

hull_ifm(x::IFMI, y::IFMI) = hull(x, y)
hull_ia(x::IAI, y::IAI) = IA.hull(x, y)                       # what `hull` maps to
hull_bare(x::IA.BareInterval{Float64}, y::IA.BareInterval{Float64}) = IA.hull(x, y)

# Endpoint min/max plus decoration propagation, without the NaI test or the
# `_set_decoration`/`setdecoration` round trip.
# IA encodes the empty `Float64` interval as `(NaN, NaN)`, so plain `min`/`max`
# would propagate the NaN rather than ignore it; `hull(::BareInterval,
# ::BareInterval)` therefore opens with an emptiness test, and that branch is
# what stops the reduction from vectorizing. `ifelse` gives the same NaN-aware
# result with no branch. (ThickNumbers' generic `hull` needs no such test:
# `emptyset` there is `(typemax, typemin)`, which min/max handle directly.)
@inline _minz(a, b) = ifelse(isnan(a), b, ifelse(isnan(b), a, min(a, b)))
@inline _maxz(a, b) = ifelse(isnan(a), b, ifelse(isnan(b), a, max(a, b)))

hull_bare_branchless(x::IA.BareInterval{Float64}, y::IA.BareInterval{Float64}) =
    IA._unsafe_bareinterval(Float64, _minz(IA.inf(x), IA.inf(y)), _maxz(IA.sup(x), IA.sup(y)))

# `inf`/`sup` are not field accesses: each normalizes NaN (empty) and, for
# `inf`, signed zero, on every read. `_inf`/`_sup` are the raw fields, so this
# variant measures IA's floor for its own representation --- everything that
# remains is the `_normalisezero` pair inside `_unsafe_bareinterval`.
hull_raw(x::IA.BareInterval{Float64}, y::IA.BareInterval{Float64}) =
    IA._unsafe_bareinterval(Float64, _minz(IA._inf(x), IA._inf(y)), _maxz(IA._sup(x), IA._sup(y)))

function hull_manual(x::IAI, y::IAI)
    bx, by = IA.bareinterval(x), IA.bareinterval(y)
    r = IA._unsafe_bareinterval(Float64, min(IA.inf(bx), IA.inf(by)), max(IA.sup(bx), IA.sup(by)))
    d = min(IA.trv, IA.decoration(x), IA.decoration(y))
    return IA._unsafe_interval(r, d, IA.isguaranteed(x) & IA.isguaranteed(y))
end

# --- reduction kernels for the profiler --------------------------------------

loop_hull_ifm(v) = reduce(hull_ifm, v)
loop_hull_ia(v) = reduce(hull_ia, v)
loop_hull_bare(v) = reduce(hull_bare, v)
loop_hull_manual(v) = reduce(hull_manual, v)
loop_hull_branchless(v) = reduce(hull_bare_branchless, v)
loop_hull_raw(v) = reduce(hull_raw, v)

# `hull_bare_branchless` must agree with `IA.hull` on every empty/nonempty
# combination, not just on the random data used for timing.
function check_branchless()
    e = IA.emptyinterval(IA.BareInterval{Float64})
    a = IA.bareinterval(Float64, 1.0, 2.0)
    b = IA.bareinterval(Float64, -3.0, 0.5)
    for x in (a, b, e), y in (a, b, e)
        r1, r2 = IA.hull(x, y), hull_bare_branchless(x, y)
        IA.isequal_interval(r1, r2) ||
            error("hull_bare_branchless disagrees on ($x, $y): $r1 vs $r2")
    end
    return true
end

# --- data --------------------------------------------------------------------

function makevec(::Type{I}, n, rng) where {I}
    out = Vector{I}(undef, n)
    for i in eachindex(out)
        a = 10 * rand(rng) - 5
        w = 3 * rand(rng)
        out[i] = I === IFMI ? IFMI(a, a + w) : IA.interval(Float64, a, a + w)
    end
    return out
end

function main(n = 1000)
    rng = Xoshiro(0xABCD1234)
    vifm = makevec(IFMI, n, rng)
    rng = Xoshiro(0xABCD1234)
    via = makevec(IAI, n, rng)
    vbare = IA.bareinterval.(via)

    x1, x2 = vifm[1], vifm[2]
    y1, y2 = via[1], via[2]
    b1, b2 = vbare[1], vbare[2]

    check_branchless()
    println("--- scalar hull (one call) ---")
    ri, ra, rb = Ref(x1), Ref(y1), Ref(b1)
    for (name, t) in (("IFM  hull(x,y)", @belapsed $ri[] = hull_ifm($x1, $x2)),
                      ("IA   hull(x,y)          full", @belapsed $ra[] = hull_ia($y1, $y2)),
                      ("IA   hull(bare,bare)", @belapsed $rb[] = hull_bare($b1, $b2)),
                      ("IA   hull_manual", @belapsed $ra[] = hull_manual($y1, $y2)),
                      ("IA   hull bare branchless", @belapsed $rb[] = hull_bare_branchless($b1, $b2)),
                      ("IA   hull raw fields", @belapsed $rb[] = hull_raw($b1, $b2)))
        println(rpad(name, 34), lpad(round(t * 1e9; digits = 2), 8), " ns")
    end

    println("\n--- reduce(hull, v) over $n intervals ---")
    for (name, t) in (("IFM  reduce", @belapsed loop_hull_ifm($vifm)),
                      ("IA   reduce             full", @belapsed loop_hull_ia($via)),
                      ("IA   reduce bare", @belapsed loop_hull_bare($vbare)),
                      ("IA   reduce manual", @belapsed loop_hull_manual($via)),
                      ("IA   reduce bare branchless", @belapsed loop_hull_branchless($vbare)),
                      ("IA   reduce raw fields", @belapsed loop_hull_raw($vbare)))
        println(rpad(name, 34), lpad(round(t * 1e9; digits = 1), 9), " ns  (",
                lpad(round(t / (n - 1) * 1e9; digits = 2), 6), " ns/element)")
    end

    # Same reduction expressed over the raw endpoints: the floor a fully
    # vectorized hull can reach on this machine.
    los = [IA.inf(x) for x in via]
    his = [IA.sup(x) for x in via]
    t = @belapsed (minimum($los), maximum($his))
    println(rpad("raw min/max over endpoints", 34), lpad(round(t * 1e9; digits = 1), 9), " ns  (",
            lpad(round(t / (n - 1) * 1e9; digits = 2), 6), " ns/element)")
    return nothing
end

main()