Outsource init code and add tests for IntervalArithmetic; unsupport :fast multiplication

src/IntervalMatrices.jl

@@ -12,54 +12,21 @@ 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
 # ================================
 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

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

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}},

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

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)])

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

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")