From 8ae8121ceac930f4849b457f52559ffb8a4a55c2 Mon Sep 17 00:00:00 2001 From: schillic Date: Mon, 12 Feb 2024 22:31:53 +0100 Subject: [PATCH 1/3] outsource init code for IntervalArithmetic --- src/IntervalMatrices.jl | 53 +++++++--------------------------- src/init_IntervalArithmetic.jl | 41 ++++++++++++++++++++++++++ 2 files changed, 51 insertions(+), 43 deletions(-) create mode 100644 src/init_IntervalArithmetic.jl diff --git a/src/IntervalMatrices.jl b/src/IntervalMatrices.jl index 0ccab0c..e16eab8 100644 --- a/src/IntervalMatrices.jl +++ b/src/IntervalMatrices.jl @@ -12,44 +12,11 @@ import Base: copy, get, similar, ∈, ⊆, ∩, ∪, real, imag using Reexport -@reexport using IntervalArithmetic -import IntervalArithmetic: inf, sup, mid, diam, radius, hull -@static if VERSION >= v"1.9" - vIA = pkgversion(IntervalArithmetic) -else - import PkgVersion - vIA = PkgVersion.Version(IntervalArithmetic) -end -if vIA >= v"0.21" - # IntervalArithmetic v0.21 removed `convert` - Base.convert(::Type{Interval{T}}, x::Number) where {T} = interval(T(x)) - Base.convert(::Type{Interval{T}}, x::Interval{T}) where {T} = x - function Base.convert(::Type{Interval{T1}}, x::Interval{T2}) where {T1,T2} - return interval(T1(inf(x)), T1(sup(x))) - end -else - # COV_EXCL_START - # IntervalArithmetic v0.21 requires `interval`, but prior versions did not - # offer this method for `Complex` inputs - IntervalArithmetic.interval(a::Complex) = complex(interval(real(a)), interval(imag(a))) - # COV_EXCL_STOP -end -if vIA >= v"0.22" - import Base: intersect - export ±, midpoint_radius - - function Base.:(==)(A::AbstractMatrix{<:Interval}, B::AbstractMatrix{<:Interval}) - return size(A) == size(B) && all(map((a, b) -> isequal_interval(a, b), A, B)) - end -else - import IntervalArithmetic: ±, midpoint_radius - - issubset_interval(x::Interval, y::Interval) = issubset(x, y) - - in_interval(x::Number, y::Interval) = inf(y) <= x <= sup(y) - - intersect_interval(a::Interval, b::Interval) = intersect(a, b) -end + +# ================================= +# Interface with IntervalArithmetic +# ================================= +include("init_IntervalArithmetic.jl") # ================================ # Type defining an interval matrix @@ -57,9 +24,9 @@ end include("matrix.jl") include("affine.jl") -# ================================= +# ================================ # Operations for interval matrices -# ================================= +# ================================ include("operations/arithmetic.jl") include("operations/mult.jl") include("operations/norm.jl") @@ -82,9 +49,9 @@ include("exponential.jl") # ============================================================== include("correction_matrices.jl") -# ======== +# ======= # Exports -# ======== +# ======= export AbstractIntervalMatrix, IntervalMatrix, set_multiplication_mode, get_multiplication_mode @@ -109,4 +76,4 @@ export AffineIntervalMatrix, set_multiplication_mode(config[:multiplication]) -end # module +end # module diff --git a/src/init_IntervalArithmetic.jl b/src/init_IntervalArithmetic.jl new file mode 100644 index 0000000..03062ae --- /dev/null +++ b/src/init_IntervalArithmetic.jl @@ -0,0 +1,41 @@ +@reexport using IntervalArithmetic +import IntervalArithmetic: inf, sup, mid, diam, radius, hull + +@static if VERSION >= v"1.9" + vIA = pkgversion(IntervalArithmetic) +else + import PkgVersion + vIA = PkgVersion.Version(IntervalArithmetic) +end + +@static if vIA >= v"0.22" + import Base: intersect + export ±, midpoint_radius # not exported by IntervalArithmetic anymore + + function Base.:(==)(A::AbstractMatrix{<:Interval}, B::AbstractMatrix{<:Interval}) + return size(A) == size(B) && all(map((a, b) -> isequal_interval(a, b), A, B)) + end +else # vIA < v"0.22" + # COV_EXCL_START + import IntervalArithmetic: ±, midpoint_radius + + issubset_interval(x::Interval, y::Interval) = issubset(x, y) + + in_interval(x::Number, y::Interval) = inf(y) <= x <= sup(y) + + intersect_interval(a::Interval, b::Interval) = intersect(a, b) + + if vIA >= v"0.21" + # `convert` was temporarily removed in IntervalArithmetic v0.21 until v0.22 + Base.convert(::Type{Interval{T}}, x::Number) where {T} = interval(T(x)) + Base.convert(::Type{Interval{T}}, x::Interval{T}) where {T} = x + function Base.convert(::Type{Interval{T1}}, x::Interval{T2}) where {T1,T2} + return interval(T1(inf(x)), T1(sup(x))) + end + else # vIA < v"0.21" + # IntervalArithmetic v0.21 requires `interval`, but prior versions did not + # offer this method for `Complex` inputs + IntervalArithmetic.interval(a::Complex) = complex(interval(real(a)), interval(imag(a))) + end + # COV_EXCL_STOP +end From 3be49993bf082d9c9f42b832c3df1c1f6761cd07 Mon Sep 17 00:00:00 2001 From: schillic Date: Mon, 12 Feb 2024 22:57:20 +0100 Subject: [PATCH 2/3] complete tests for IntervalArithmetic interface --- test/constructors.jl | 12 --------- test/intervals.jl | 58 ++++++++++++++++++++++++++++++++++++++++++++ test/runtests.jl | 2 ++ 3 files changed, 60 insertions(+), 12 deletions(-) create mode 100644 test/intervals.jl diff --git a/test/constructors.jl b/test/constructors.jl index 12cc0d0..41e944a 100644 --- a/test/constructors.jl +++ b/test/constructors.jl @@ -1,15 +1,3 @@ -@testset "Interval conversion" begin - x = interval(1.0) - - y = convert(Interval{Float64}, x) - # `===` is not applicable here because it just checks value equivalence - # for (immutable) `Interval`s - @test y ⩵ x && y isa Interval{Float64} - - y = convert(Interval{Float32}, x) - @test y ⩵ x && y isa Interval{Float32} -end - @testset "Interval matrix construction" begin m1 = IntervalMatrix([interval(-1.1, 0.9) interval(-4.1, -3.9); interval(3.8, 4.2) interval(0, 0.9)]) diff --git a/test/intervals.jl b/test/intervals.jl new file mode 100644 index 0000000..0985d62 --- /dev/null +++ b/test/intervals.jl @@ -0,0 +1,58 @@ +# =============================================== +# helper methods for IntervalArithmetic interface +# =============================================== +using IntervalMatrices: intersect_interval, in_interval, issubset_interval + +@testset "Interval conversion" begin + x = convert(Interval{Float64}, 1.0) + + y = convert(Interval{Float64}, x) + # `===` is not applicable here because it just checks value equivalence + # for (immutable) `Interval`s + @test y ⩵ x && y isa Interval{Float64} + + y = convert(Interval{Float32}, x) + @test y ⩵ x && y isa Interval{Float32} +end + +@testset "Interval operation `intersect_interval`" begin + x = interval(1, 3) + for y in (interval(2, 4), interval(2f0, 4f0)) + @test intersect_interval(x, y) ⩵ interval(2, 3) + end + @test intersect_interval(x, x) ⩵ x + for y in (interval(4, 5), interval(4f0, 5f0)) + @test intersect_interval(x, y) ⩵ emptyinterval(Float64) + end +end + +@testset "Interval operation `in_interval`" begin + x = interval(1, 3) + for e in (1, 2.0, 3) + @test in_interval(e, x) + end + for e in (0, 0.999, 3.001, 4) + @test !in_interval(e, x) + end +end + +@testset "Interval operation `issubset_interval`" begin + x = interval(1, 3) + for y in (interval(0, 4), x) + @test issubset_interval(x, y) + end + for y in (interval(2, 4), interval(0, 2), interval(2)) + @test !issubset_interval(x, y) + end +end + +@testset "Equality for `Interval` matrix" begin + M = [interval(1) interval(2); interval(3) interval(4)] + @test M == [1 2; 3 4] + M = [interval(1, 2) interval(2, 3); interval(3, 4) interval(4, 5)] + @test M == M +end + +@testset "Interval from a complex number" begin + @test interval(1+2im) isa Complex{Interval{Float64}} +end diff --git a/test/runtests.jl b/test/runtests.jl index 05fa594..5d6f69c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -12,10 +12,12 @@ using IntervalMatrices: _truncated_exponential_series, ⩵(x::Interval, y::Interval) = isequal_interval(x, y) else ⩵(x::Interval, y::Interval) = ==(x, y) + using IntervalMatrices: issubset_interval end include("models.jl") +include("intervals.jl") include("constructors.jl") include("arithmetic.jl") include("setops.jl") From 97030e07c82b2827088b2fc3bdffdd8eff8c0ab3 Mon Sep 17 00:00:00 2001 From: schillic Date: Mon, 12 Feb 2024 23:08:18 +0100 Subject: [PATCH 3/3] do not support :fast multiplication with IA v0.22 --- src/operations/mult.jl | 120 ++++++++++++++++++++++------------------- test/arithmetic.jl | 2 + 2 files changed, 68 insertions(+), 54 deletions(-) diff --git a/src/operations/mult.jl b/src/operations/mult.jl index a711789..227aa0f 100644 --- a/src/operations/mult.jl +++ b/src/operations/mult.jl @@ -16,15 +16,22 @@ Sets the algorithm used to perform matrix multiplication with interval matrices. - `:fast` -- computes an enclosure of the matrix product using the midpoint-radius notation of the matrix [[RUM10]](@ref). +!!! note ":fast option no longer supported" + `:fast` support was removed in `IntervalArithmetic` v0.22. + ### Notes -- By default, `:fast` is used. +- By default, `:slow` is used. - Using `fast` is generally significantly faster, but it may return larger intervals, especially if midpoint and radius have the same order of magnitude (50% overestimate at most) [[RUM99]](@ref). """ function set_multiplication_mode(multype) - multype in (:fast, :slow) || throw(ArgumentError("$multype is not a valid input.")) + multype in (:fast, :slow) || throw(ArgumentError("$multype is not a valid input")) + + @static if vIA >= v"0.22" + multype == :fast && throw(ArgumentError("$multype is not supported anymore")) + end type = MultiplicationType{multype}() @eval *(A::IntervalMatrix, B::IntervalMatrix) = *($type, A, B) @@ -40,81 +47,86 @@ end *(::MultiplicationType{:slow}, A::AbstractMatrix, B::IntervalMatrix) = IntervalMatrix(A * B.mat) *(::MultiplicationType{:slow}, A::IntervalMatrix, B::AbstractMatrix) = IntervalMatrix(A.mat * B) -function *(::MultiplicationType{:fast}, - A::AbstractIntervalMatrix{Interval{T}}, - B::AbstractIntervalMatrix{Interval{T}}) where {T} - Ainf = inf(A) - Asup = sup(A) - Binf = inf(B) - Bsup = sup(B) +# not used anymore since IntervalArithmetic v0.22 +# COV_EXCL_START +@static if vIA < v"0.22" + function *(::MultiplicationType{:fast}, + A::AbstractIntervalMatrix{Interval{T}}, + B::AbstractIntervalMatrix{Interval{T}}) where {T} + Ainf = inf(A) + Asup = sup(A) + Binf = inf(B) + Bsup = sup(B) - mA, mB, R, Csup = setrounding(T, RoundUp) do - mA = Ainf + 0.5 * (Asup - Ainf) - mB = Binf + 0.5 * (Bsup - Binf) + mA, mB, R, Csup = setrounding(T, RoundUp) do + mA = Ainf + 0.5 * (Asup - Ainf) + mB = Binf + 0.5 * (Bsup - Binf) - rA = mA - Ainf - rB = mB - Binf + rA = mA - Ainf + rB = mB - Binf - R = abs.(mA) * rB + rA * (abs.(mB) + rB) - Csup = mA * mB + R + R = abs.(mA) * rB + rA * (abs.(mB) + rB) + Csup = mA * mB + R - return mA, mB, R, Csup - end + return mA, mB, R, Csup + end + + Cinf = setrounding(T, RoundDown) do + return mA * mB - R + end - Cinf = setrounding(T, RoundDown) do - return mA * mB - R + return IntervalMatrix(interval.(Cinf, Csup)) end - return IntervalMatrix(interval.(Cinf, Csup)) -end + function *(::MultiplicationType{:fast}, + A::AbstractMatrix{T}, + B::AbstractIntervalMatrix{Interval{T}}) where {T} + Binf = inf(B) + Bsup = sup(B) -function *(::MultiplicationType{:fast}, - A::AbstractMatrix{T}, - B::AbstractIntervalMatrix{Interval{T}}) where {T} - Binf = inf(B) - Bsup = sup(B) + mB, R, Csup = setrounding(T, RoundUp) do + mB = Binf + 0.5 * (Bsup - Binf) - mB, R, Csup = setrounding(T, RoundUp) do - mB = Binf + 0.5 * (Bsup - Binf) + rB = mB - Binf - rB = mB - Binf + R = abs.(A) * rB + Csup = A * mB + R - R = abs.(A) * rB - Csup = A * mB + R + return mB, R, Csup + end - return mB, R, Csup - end + Cinf = setrounding(T, RoundDown) do + return A * mB - R + end - Cinf = setrounding(T, RoundDown) do - return A * mB - R + return IntervalMatrix(interval.(Cinf, Csup)) end - return IntervalMatrix(interval.(Cinf, Csup)) -end + function *(::MultiplicationType{:fast}, + A::AbstractIntervalMatrix{Interval{T}}, + B::AbstractMatrix{T}) where {T} + Ainf = inf(A) + Asup = sup(A) -function *(::MultiplicationType{:fast}, - A::AbstractIntervalMatrix{Interval{T}}, - B::AbstractMatrix{T}) where {T} - Ainf = inf(A) - Asup = sup(A) + mA, R, Csup = setrounding(T, RoundUp) do + mA = Ainf + 0.5 * (Asup - Ainf) - mA, R, Csup = setrounding(T, RoundUp) do - mA = Ainf + 0.5 * (Asup - Ainf) + rA = mA - Ainf - rA = mA - Ainf + R = rA * abs.(B) + Csup = mA * B + R - R = rA * abs.(B) - Csup = mA * B + R + return mA, R, Csup + end - return mA, R, Csup - end + Cinf = setrounding(T, RoundDown) do + return mA * B - R + end - Cinf = setrounding(T, RoundDown) do - return mA * B - R + return IntervalMatrix((interval.(Cinf, Csup))) end - - return IntervalMatrix((interval.(Cinf, Csup))) end +# COV_EXCL_STOP # function *(::MultiplicationType{:rank1}, # A::AbstractMatrix{Interval{T}}, diff --git a/test/arithmetic.jl b/test/arithmetic.jl index 33ba43a..9cfa1ec 100644 --- a/test/arithmetic.jl +++ b/test/arithmetic.jl @@ -88,5 +88,7 @@ end IntervalMatrix([interval(5, 12.5) interval(-8, 2); interval(-2, 8) interval(5, 12.5)]) @test mid.(A) * A == IntervalMatrix([interval(5, 12.5) interval(-8, 2); interval(-2, 8) interval(5, 12.5)]) + else + @test_throws ArgumentError set_multiplication_mode(:fast) end end