Support IntervalArithmetic v0.22

Project.toml

@@ -10,7 +10,7 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 
 [compat]
-IntervalArithmetic = "0.15, 0.16, 0.17, 0.18, 0.19, 0.20, 0.21"
+IntervalArithmetic = "0.15, 0.16, 0.17, 0.18, 0.19, 0.20, 0.21, 0.22"
 PkgVersion = "0.3.3"
 Reexport = "0.2, 1"
 julia = "1.1"

docs/Project.toml

@@ -1,5 +1,7 @@
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
 
 [compat]
 Documenter = "1"
+IntervalArithmetic = "0.21"

docs/src/index.md

@@ -51,7 +51,7 @@ An *interval matrix* is a matrix whose coefficients are intervals. For instance,
 ```jldoctest quickstart
 julia> using IntervalMatrices
 
-julia> A = IntervalMatrix([0..1 1..2; 2..3 -4.. -2])
+julia> A = IntervalMatrix([interval(0, 1) interval(1, 2); interval(2, 3) interval(-4, -2)])
 2×2 IntervalMatrix{Float64, Interval{Float64}, Matrix{Interval{Float64}}}:
  [0.0, 1.0]   [1.0, 2.0]
  [2.0, 3.0]  [-4.0, -2.0]
@@ -81,7 +81,7 @@ julia> 2A
 Or an interval multiple of $A$,
 
 ```jldoctest quickstart
-julia> (-1.0..1.0) * A
+julia> interval(-1, 1) * A
 2×2 IntervalMatrix{Float64, Interval{Float64}, Matrix{Interval{Float64}}}:
  [-1.0, 1.0]  [-2.0, 2.0]
  [-3.0, 3.0]  [-4.0, 4.0]
@@ -107,8 +107,8 @@ $e^{At} - I$. Then, at $t = 1.0$,
 ```jldoctest quickstart
 julia> A + 1/2 * A^2
 2×2 IntervalMatrix{Float64, Interval{Float64}, Matrix{Interval{Float64}}}:
-  [0.5, 4.50001]  [-3.25001, 2.50001]
- [-4.25001, 3.0]  [-2.25001, 9.0]
+  [1.0, 4.50001]  [-3.0, 2.0]
+ [-4.0, 2.50001]  [-1.0, 9.0]
 ```
 However, that result is not tight. The computation can be performed exactly via
 single-use expressions implemented in this library:

src/IntervalMatrices.jl

@@ -24,6 +24,34 @@ else
     # IntervalArithmetic v0.21 requires interval, but prior versions did not offer this method
     IntervalArithmetic.interval(a::Complex) = complex(interval(real(a)), interval(imag(a)))
 end
+if vIA >= v"0.22"
+    import Base: intersect
+    export ±, midpoint_radius
+
+    function ±(x::Number, y::Number)
+        return x + interval(-y, y)
+    end
+
+    function Base.:(==)(A::AbstractMatrix{<:Interval}, B::AbstractMatrix{<:Interval})
+        return size(A) == size(B) && all(map((a, b) -> isequal_interval(a, b), A, B))
+    end
+
+    function intersect(a::Interval{T}, b::Interval{T}) where T
+        lo = max(inf(a), inf(b))
+        hi = min(sup(a), sup(b))
+        if lo > hi
+            return emptyinterval(T)
+        end
+        return interval(lo, hi)
+    end
+    intersect(a::Interval{T}, b::Interval{S}) where {T,S} = intersect(promote(a, b)...)
+
+    Base.in(x::Number, y::Interval) = inf(y) <= x <= sup(y)
+else
+    import IntervalArithmetic: ±, midpoint_radius
+
+    issubset_interval(x::Interval, y::Interval) = issubset(x, y)
+end
 
 # ================================
 # Type defining an interval matrix

src/affine.jl

@@ -16,12 +16,12 @@ where ``A₀`` and ``A₁`` are real (or complex) matrices, and ``λ`` is an int
 
 ### Examples
 
-The matrix ``I + [1 1; -1 1] * (0 .. 1)`` is:
+The matrix ``I + [1 1; -1 1] * interval(0, 1)`` is:
 
 ```jldoctest
 julia> using LinearAlgebra
 
-julia> P = AffineIntervalMatrix1(Matrix(1.0I, 2, 2), [1 1; -1 1.], 0 .. 1);
+julia> P = AffineIntervalMatrix1(Matrix(1.0I, 2, 2), [1 1; -1 1.], interval(0, 1));
 
 julia> P
 2×2 AffineIntervalMatrix1{Float64, Interval{Float64}, Matrix{Float64}, Matrix{Float64}}:
@@ -76,7 +76,7 @@ contains one matrix proportional to an interval.
 
 ### Examples
 
-The affine matrix ``I + [1 1; -1 1] * (0 .. 1) + [0 1; 1 0] * (2 .. 3)`` is:
+The affine matrix ``I + [1 1; -1 1] * interval(0, 1) + [0 1; 1 0] * interval(2, 3)`` is:
 
 ```jldoctest
 julia> using LinearAlgebra
@@ -85,7 +85,7 @@ julia> A0 = Matrix(1.0I, 2, 2);
 
 julia> A1 = [1 1; -1 1.]; A2 = [0 1; 1 0];
 
-julia> λ1 = 0 .. 1; λ2 = 2 .. 3;
+julia> λ1 = interval(0, 1); λ2 = interval(2, 3);
 
 julia> P = AffineIntervalMatrix(A0, [A1, A2], [λ1, λ2])
 2×2 AffineIntervalMatrix{Float64, Interval{Float64}, Matrix{Float64}, Matrix{Float64}, Vector{Matrix{Float64}}, Vector{Interval{Float64}}}:

src/matrix.jl

@@ -1,6 +1,6 @@
 import Base: similar, ∈, ⊆, ∩, ∪, real, imag
 import Random: rand
-import IntervalArithmetic: ±, inf, sup, mid, diam, radius, midpoint_radius, hull
+import IntervalArithmetic: inf, sup, mid, diam, radius, hull
 
 """
     AbstractIntervalMatrix{IT} <: AbstractMatrix{IT}

src/operations/mult.jl

@@ -1,4 +1,4 @@
-const config = Dict(:multiplication => :fast)
+const config = Dict(:multiplication => :slow)
 
 struct MultiplicationType{T} end
 

src/operations/power.jl

@@ -19,7 +19,7 @@ The internal representation (such as the fields) are subject to future changes.
 ### Examples
 
 ```jldoctest
-julia> A = IntervalMatrix([2.0..2.0 2.0..3.0; 0.0..0.0 -1.0..1.0])
+julia> A = IntervalMatrix([interval(2, 2) interval(2, 3); interval(0, 0) interval(-1, 1)])
 2×2 IntervalMatrix{Float64, Interval{Float64}, Matrix{Interval{Float64}}}:
  [2.0, 2.0]   [2.0, 3.0]
  [0.0, 0.0]  [-1.0, 1.0]

src/operations/setops.jl

@@ -52,7 +52,7 @@ function ⊆(A::AbstractIntervalMatrix, B::AbstractIntervalMatrix)
 
     m, n = size(A)
     @inbounds for j in 1:n, i in 1:m
-        if !(A[i, j] ⊆ B[i, j])
+        if !issubset_interval(A[i, j], B[i, j])
             return false
         end
     end

test/Project.toml

@@ -1,4 +1,5 @@
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+PkgVersion = "eebad327-c553-4316-9ea0-9fa01ccd7688"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

test/affine.jl

@@ -1,7 +1,7 @@
 @testset "Affine interval matrix in one interval" begin
     A0 = [1 0; 0 1.0]
     A1 = [0 1; 1 0.0]
-    λ = 0 .. 1
+    λ = interval(0, 1)
     P = AffineIntervalMatrix1(A0, A1, λ)
 
     @test size(P) == (2, 2)
@@ -23,7 +23,7 @@ end
     A0 = [1 0; 0 1.0]
     A1 = [0 1; 1 0.0]
     A2 = copy(A1)
-    λ1 = 0 .. 1
+    λ1 = interval(0, 1)
     λ2 = copy(λ1)
     P = AffineIntervalMatrix(A0, [A1, A2], [λ1, λ2])
 

test/arithmetic.jl

@@ -1,42 +1,42 @@
 @testset "Interval arithmetic" begin
-    a = -1.5 ± 0.5
-    b = -1 .. 1
-    @test a + b == -3.0 .. 0.0
-    @test a * b == -2.0 .. 2.0
-    @test a * b + a == -4.0 .. 1.0
-    @test a * (b + 1) == -4.0 .. 0.0
+    a = interval(-2, -1)
+    b = interval(-1, 1)
+    @test a + b == interval(-3, 0)
+    @test a * b == interval(-2, 2)
+    @test a * b + a == interval(-4, 1)
+    @test a * (b + 1) == interval(-4, 0)
 end
 
 @testset "Interval matrix arithmetic" begin
-    a = 1.0 .. 1.3
-    b = 2.0 .. 3.5
-    c = -0.5 ± 0.5
-    d = 0.0 ± 0.1
-    a₊ = 2.0 .. 2.6
-    b₊ = 4.0 .. 7.0
-    c₊ = -2.0 .. 0.0
-    d₊ = -0.2 .. 0.2
-    a₋ = -0.3 .. 0.3
-    b₋ = -1.5 .. 1.5
-    c₋ = -1.0 .. 1.0
-    d₋ = -0.2 .. 0.2
+    a = interval(1, 1.3)
+    b = interval(2, 3.5)
+    c = interval(-1, 0)
+    d = interval(-0.1, 0.1)
+    a₊ = interval(2, 2.6)
+    b₊ = interval(4, 7)
+    c₊ = interval(-2, 0)
+    d₊ = interval(-0.2, 0.2)
+    a₋ = interval(-0.3, 0.3)
+    b₋ = interval(-1.5, 1.5)
+    c₋ = interval(-1, 1)
+    d₋ = interval(-0.2, 0.2)
     A = IntervalMatrix([a b; c d])
 
     B = A + A
     @test B isa IntervalMatrix && B == IntervalMatrix([a₊ b₊; c₊ d₊])
 
     B = A - A
-    @test B isa IntervalMatrix && B[1, 1].lo ≈ a₋.lo && B[1, 1].hi ≈ a₋.hi &&
+    @test B isa IntervalMatrix && inf(B[1, 1]) ≈ inf(a₋) && sup(B[1, 1]) ≈ sup(a₋) &&
           B[1, 2] == b₋ && B[2, 1] == c₋ && B[2, 2] == d₋
 
     # arithmetic with an interval or number
-    for x in (0.0 .. 2.0, 2.0)
+    for x in (interval(0, 2), 2.0)
         @test A + x == x + A == IntervalMatrix([a+x b+x; c+x d+x])
     end
 
     # multiply interval and interval matrix
-    x = 0.0 .. 2.0
-    Ax = IntervalMatrix([0.0..2.6 0.0..7.0; -2.0..0.0 -0.2..0.2])
+    x = interval(0, 2)
+    Ax = IntervalMatrix([interval(0, 2.6) interval(0, 7); interval(-2, 0) interval(-0.2, 0.2)])
     @test x * A == A * x == Ax
     # multiply scalar and interval matrix
     x = 1.0
@@ -57,26 +57,29 @@ end
     @test Ainf - A isa IntervalMatrix
     @test A * Ainf isa IntervalMatrix
     @test Ainf * A isa IntervalMatrix
-    @test A \ Ainf isa IntervalMatrix
-    @test Ainf \ A isa IntervalMatrix
+    @static if PkgVersion.Version(IntervalMatrices.IntervalArithmetic) < v"0.22"
+        @test A \ Ainf isa IntervalMatrix
+        @test Ainf \ A isa IntervalMatrix
+    end
 end
 
 @testset "Matrix multiplication" begin
-
     # test default settings
-    @test get_multiplication_mode() == Dict(:multiplication => :fast)
+    @test get_multiplication_mode() == Dict(:multiplication => :slow)
 
-    A = IntervalMatrix([2..4 -2..1; -1..2 2..4])
+    A = IntervalMatrix([interval(2, 4) interval(-2, 1); interval(-1, 2) interval(2, 4)])
     set_multiplication_mode(:slow)
-    @test A * A == IntervalMatrix([0..18 -16..8; -8..16 0..18])
-    @test A * mid.(A) == IntervalMatrix([5..12.5 -8..2; -2..8 5..12.5])
-    @test mid.(A) * A == IntervalMatrix([5..12.5 -8..2; -2..8 5..12.5])
+    @test A * A == IntervalMatrix([interval(0, 18) interval(-16, 8); interval(-8, 16) interval(0, 18)])
+    @test A * mid.(A) == 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)])
 
     # set_multiplication_mode(:rank1)
-    # @test A * A == [0..18 -16..8; -8..16 0..18]
+    # @test A * A == [interval(0, 18) interval(-16, 8); interval(-8, 16) interval(0, 18)]
 
-    set_multiplication_mode(:fast)
-    @test A * A == IntervalMatrix([-2..19.5 -16..10; -10..16 -2..19.5])
-    @test A * mid.(A) == IntervalMatrix([5..12.5 -8..2; -2..8 5..12.5])
-    @test mid.(A) * A == IntervalMatrix([5..12.5 -8..2; -2..8 5..12.5])
+    @static if PkgVersion.Version(IntervalMatrices.IntervalArithmetic) < v"0.22"
+        set_multiplication_mode(:fast)
+        @test A * A == IntervalMatrix([interval(-2, 19.5) interval(-16, 10); interval(-10, 16) interval(-2, 19.5)])
+        @test A * mid.(A) == 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)])
+    end
 end

test/constructors.jl

@@ -1,5 +1,5 @@
 @testset "Interval matrix construction" begin
-    m1 = IntervalMatrix([-1.1..0.9 -4.1 .. -3.9; 3.8..4.2 0.0..0.9])
+    m1 = IntervalMatrix([interval(-1.1, 0.9) interval(-4.1, -3.9); interval(3.8, 4.2) interval(0, 0.9)])
     m2 = IntervalMatrix{Float64}(undef, 2, 2)
     @test m2 isa IntervalMatrix{Float64} && size(m2) == (2, 2)
     m3 = similar(m1)
@@ -11,20 +11,20 @@
     A = [1 2; 3 4]
     B = [1 2; 4 5]
 
-    @test IntervalMatrix(A, B) == IntervalMatrix([1..1 2..2; 3..4 4..5])
+    @test IntervalMatrix(A, B) == IntervalMatrix([interval(1, 1) interval(2, 2); interval(3, 4) interval(4, 5)])
 
-    @test A ± B == IntervalMatrix([0..2 0..4; -1..7 -1..9])
+    @test A ± B == IntervalMatrix([interval(0, 2) interval(0, 4); interval(-1, 7) interval(-1, 9)])
 end
 
 @testset "Interval matrix midpoint_radius" begin
-    m = IntervalMatrix([-1.1..0.9 -4.1 .. -3.9; 3.9..4.1 -1.1..0.9])
+    m = IntervalMatrix([interval(-1.1, 0.9) interval(-4.1, -3.9); interval(3.9, 4.1) interval(-1.1, 0.9)])
     c, s = midpoint_radius(m)
     @test c ≈ [-0.1 -4; 4 -0.1]
     @test s ≈ [1 0.1; 0.1 1]
 end
 
 @testset "Complex interval matrices" begin
-    m1 = IntervalMatrix([(1 .. 2)+im * (3 .. 4) 1; 2 3])
+    m1 = IntervalMatrix([interval(1, 2) + im * interval(3, 4) 1; 2 3])
     @test m1 isa IntervalMatrix{Float64}
     @test eltype(m1) == Complex{Interval{Float64}}
 

test/exponential.jl

@@ -8,9 +8,9 @@ using IntervalMatrices: TaylorOverapproximation,
                         _exp_remainder_series
 
 @testset "Interval matrix exponential" begin
-    @test quadratic_expansion(-3 .. 3, 1.0, 2.0) == interval(-0.125, 21)
+    @test quadratic_expansion(interval(-3, 3), 1.0, 2.0) == interval(-0.125, 21)
 
-    M = IntervalMatrix([-1.1..0.9 -4.1 .. -3.9; 3.9..4.1 -1.1..0.9])
+    M = IntervalMatrix([interval(-1.1, 0.9) interval(-4.1, -3.9); interval(3.9, 4.1) interval(-1.1, 0.9)])
 
     for i in 0:4
         _truncated_exponential_series(M, 1.0, i)
@@ -36,22 +36,22 @@ using IntervalMatrices: TaylorOverapproximation,
 end
 
 @testset "Interval matrix correction terms" begin
-    m = IntervalMatrix([-1.1..0.9 -4.1 .. -3.9; 3.9..4.1 -1.1..0.9])
+    m = IntervalMatrix([interval(-1.1, 0.9) interval(-4.1, -3.9); interval(3.9, 4.1) interval(-1.1, 0.9)])
     f = correction_hull(m, 1e-3, 5)
     f2 = input_correction(m, 1e-3, 5)
     f = correction_hull(mid(m), 1e-3, 5)
 end
 
 @testset "Interval matrix square" begin
-    m = IntervalMatrix([-1.1..0.9 -4.1 .. -3.9; 3.9..4.1 -1.1..0.9])
+    m = IntervalMatrix([interval(-1.1, 0.9) interval(-4.1, -3.9); interval(3.9, 4.1) interval(-1.1, 0.9)])
 
     a = m * m
     b = square(m)
     @test b ⊆ a
 end
 
 @testset "Interval matrix power" begin
-    m = IntervalMatrix([2.0..2.0 2.0..3.0; 0.0..0.0 -1.0..1.0])
+    m = IntervalMatrix([interval(2, 2) interval(2, 3); interval(0, 0) interval(-1, 1)])
     pow = IntervalMatrixPower(m)
 
     @test base(pow) === m
@@ -85,7 +85,11 @@ end
     A = Matrix(transmission_line())
     expA = exp(A)
     @test opnorm(expA, Inf) < 1e-6
-    Aint = IntervalMatrix(interval.(A))
-    expA_int = exp_overapproximation(Aint, 1.0, 10)
-    @test opnorm(expA_int, Inf) > 1e122
+
+    @static if PkgVersion.Version(IntervalMatrices.IntervalArithmetic) < v"0.22"
+        # in newer versions, too large numbers lead to empty interval
+        Aint = IntervalMatrix(interval.(A))
+        expA_int = exp_overapproximation(Aint, 1.0, 10)
+        @test opnorm(expA_int, Inf) > 1e122
+    end
 end

test/runtests.jl

@@ -1,10 +1,17 @@
 using IntervalMatrices, Test, LinearAlgebra, SparseArrays
+import PkgVersion
 
 using IntervalMatrices: _truncated_exponential_series,
                         horner,
                         scale_and_square,
                         correction_hull,
                         input_correction
+
+@static if PkgVersion.Version(IntervalMatrices.IntervalArithmetic) >= v"0.22"
+    # equality test for convenience
+    Base.:(==)(x::Interval, y::Interval) = isequal_interval(x, y)
+end
+
 include("models.jl")
 
 include("constructors.jl")