Update to LazySet v6

.github/CONTRIBUTING.md

@@ -31,7 +31,7 @@ This is something that needs to be done only once, the first time you start cont
 
 **1.** From the Julia REPL in package mode (you can enter package mode by typing `]`) do
 
-```julia
+```jldoctest
 pkg> dev IntervalLinearAlgebra
 ```
 
@@ -149,7 +149,7 @@ If the dependency is quite heavy and used only by some functionalities, you may
 2. In the `[deps]` section of `Project.toml` locate the package you want to make an optional dependency and move the corresponding line to `[extras]`, keep alphabetical ordering.
 3. Add the dependency name to the `test` entry in the `[targets]` section
 4. In the `IntervalLinearAlgebra.jl` file, locate the `__init__` function and add the line
-```julia
+```jldoctest
 @require """Example = "7876af07-990d-54b4-ab0e-23690620f79a" include("file.jl")"""
 ```
 where `file.jl` is the file containing the functions needing `Example.jl`. The line `Example = "7876af07-990d-54b4-ab0e-23690620f79a"` is the same in the Project.toml
@@ -158,10 +158,10 @@ where `file.jl` is the file containing the functions needing `Example.jl`. The l
 ## Documentation guideline
 
 * Documentation is written with [Documenter.jl](https://github.com/JuliaDocs/Documenter.jl). Documentation files are in `docs/src`, generally as markdown file.
-* If you want to include a Julia code example that is **not** executed in the markdown file, use ````` ```julia ````` blocks, e.g.
+* If you want to include a Julia code example that is **not** executed in the markdown file, use ````` ```jldoctest ````` blocks, e.g.
 
 ````
-```julia
+```jldoctest
 a = 1
 b = 2
 ```
@@ -241,7 +241,7 @@ INSERT BIBTEX HERE
 ### Docstrings
 
 * Each exported function should have a docstring. The docstring should roughly follow the following structure
-```julia
+```jldoctest
 """
     funname(param1, param2[, optional_param])
 
@@ -276,7 +276,7 @@ Preferably, as a doctest.
 
 Here is an example
 
-````julia
+````jldoctest
 """
     something(A::Matrix{T}, b::Vector{T}[, tol=1e-10]) where {T<:Interval}
 

.github/workflows/CI.yml

@@ -66,7 +66,7 @@ jobs:
             julia --project=docs -e '
               using Documenter: DocMeta, doctest
               using IntervalLinearAlgebra
-              DocMeta.setdocmeta!(IntervalLinearAlgebra, :DocTestSetup, :(using IntervalLinearAlgebra); recursive=true)
+              DocMeta.setdocmeta!(IntervalLinearAlgebra, :DocTestSetup, :(using IntervalLinearAlgebra, IntervalArithmetic, IntervalArithmetic.Symbols); recursive=true)
               doctest(IntervalLinearAlgebra)'
         - run: julia --project=docs docs/make.jl
           env:

Project.toml

@@ -1,36 +1,37 @@
 name = "IntervalLinearAlgebra"
 uuid = "92cbe1ac-9c24-436b-b0c9-5f7317aedcd5"
-authors = ["Luca Ferranti"]
 version = "0.2.0"
+authors = ["Luca Ferranti"]
 
 [deps]
 CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"
 IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
-IntervalConstraintProgramming = "138f1668-1576-5ad7-91b9-7425abbf3153"
-LazySets = "b4f0291d-fe17-52bc-9479-3d1a343d9043"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 OpenBLASConsistentFPCSR_jll = "6cdc7f73-28fd-5e50-80fb-958a8875b1af"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
+SetRounding = "3cc68bcd-71a2-5612-b932-767ffbe40ab0"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
-[extensions]
-IntervalConstraintProgrammingExt = "IntervalConstraintProgramming"
-LazySetsExt = "LazySets"
-
 [weakdeps]
+IntervalBoxes = "43d83c95-ebbb-40ec-8188-24586a1458ed"
 IntervalConstraintProgramming = "138f1668-1576-5ad7-91b9-7425abbf3153"
 LazySets = "b4f0291d-fe17-52bc-9479-3d1a343d9043"
 
+[extensions]
+IntervalConstraintProgrammingExt = ["IntervalBoxes", "IntervalConstraintProgramming"]
+LazySetsExt = "LazySets"
+
 [compat]
 CommonSolve = "0.2"
-IntervalArithmetic = "0.20.5, 1.0"
-IntervalConstraintProgramming = "0.13, 0.14"
-LazySets = "2 - 5"
+IntervalArithmetic = "1"
+IntervalBoxes = "0.3"
+IntervalConstraintProgramming = "0.13, 0.14, 0.15"
+LazySets = "2 - 6"
 OpenBLASConsistentFPCSR_jll = "0.3.29 - 0.3.30"
 Reexport = "1"
+SetRounding = "0.2"
 StaticArrays = "0.12, 1"
 julia = "1.10"
 
 [extras]
-IntervalConstraintProgramming = "138f1668-1576-5ad7-91b9-7425abbf3153"
 LazySets = "b4f0291d-fe17-52bc-9479-3d1a343d9043"

docs/Project.toml

@@ -1,5 +1,7 @@
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
+IntervalLinearAlgebra = "92cbe1ac-9c24-436b-b0c9-5f7317aedcd5"
 LazySets = "b4f0291d-fe17-52bc-9479-3d1a343d9043"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"

docs/make.jl

@@ -15,7 +15,7 @@ for (root, _, files) in walkdir(litdir)
         end
     end
 end
-DocMeta.setdocmeta!(IntervalLinearAlgebra, :DocTestSetup, :(using IntervalLinearAlgebra); recursive=true)
+DocMeta.setdocmeta!(IntervalLinearAlgebra, :DocTestSetup, :(using IntervalLinearAlgebra, IntervalArithmetic, IntervalArithmetic.Symbols); recursive=true)
 
 makedocs(;
     modules=[IntervalLinearAlgebra],

docs/src/CONTRIBUTING.md

@@ -304,7 +304,7 @@ The function uses the *something sometimes somewhere* algorithm proposed by Some
 
 ### Example
 
-```jldoctest
+```julia
 julia> A = [1..2 3..4;5..6 7..8]
 2×2 Matrix{Interval{Float64}}:
 [1, 2]  [3, 4]

ext/IntervalConstraintProgrammingExt.jl

@@ -7,61 +7,38 @@ else
     using IntervalLinearAlgebra
     using IntervalConstraintProgramming
 end
+using IntervalBoxes: IntervalBox
 
-"""
-returns the unrolled expression for \$|a ⋅x - b|\$
-\$a\$ and \$x\$ must be vectors of the same length and \$b\$ is scalar.
-The absolue value in the equation is taken elementwise.
-"""
-function oettli_lhs(a, b, x)
-    ex = :( $(a[1])*$(x[1]))
-    for i = 2:length(x)
-        ex = :( $ex + $(a[i])*$(x[i]))
-    end
-    return :(abs($ex - $b))
-end
+# Access Symbolics through IntervalConstraintProgramming's dependency
+const Symbolics = IntervalConstraintProgramming.Symbolics
+const Num = Symbolics.Num
 
 """
-returns the unrolled expression for \$a ⋅|x| + b\$
-\$a\$ and \$x\$ must be vectors of the same length and \$b\$ is scalar.
-The absolue value in the equation is taken elementwise.
-"""
-function oettli_rhs(a, b, x)
-    ex = :( $(a[1])*abs($(x[1])))
-    for i = 2:length(x)
-        ex = :( $ex + $(a[i])*abs($(x[i])))
-    end
-    return :($ex + $b)
-end
-
-"""
-Returns the separator for the constraint `|a_c ⋅x - b_c| - a_r ⋅|x| - b_r <= 0`.
-
-`a` and `x` must be vectors of the same length and `b` is scalar.
-
-The absolue values in the equation are taken elementwise.
+Returns the separator for the Oettli-Prager constraint for the `i`-th row:
+`|a_c ⋅x - b_c| - a_r ⋅|x| - b_r <= 0`.
 
-`a_c` and `a_r` are vectors containing midpoints and radii of the intervals in `a`. Similar `b_c` and `b_r`.
+`a_c` and `a_r` are vectors containing midpoints and radii of the intervals in `a`.
+Similarly `b_c` and `b_r`.
 """
 function oettli_eq(a, b, x)
     ac = mid.(a)
     ar = IntervalArithmetic.radius.(a)
     bc = mid(b)
     br = IntervalArithmetic.radius(b)
-    lhs = oettli_lhs(ac, bc, x)
-    rhs = oettli_rhs(ar, br, x)
-    ex = :(@constraint $lhs - $rhs <= 0)
-    @eval $ex
+    lhs = abs(sum(ac[i] * x[i] for i in eachindex(x)) - bc)
+    rhs = sum(ar[i] * abs(x[i]) for i in eachindex(x)) + br
+    return constraint(lhs - rhs <= 0)
 end
 
 function (op::IntervalLinearAlgebra.NonLinearOettliPrager)(A, b, X::IntervalBox)
-    vars = ntuple(i -> Symbol(:x, i), length(b))
-    separators = [oettli_eq(A[i,:], b[i], vars) for i in 1:length(b)]
+    n = length(b)
+    vars = [Num(Symbolics.Sym{Real}(Symbol(:x, i))) for i in 1:n]
+    separators = [oettli_eq(A[i,:], b[i], vars) for i in 1:n]
     S = reduce(∩, separators)
     return Base.invokelatest(pave, S, X, op.tol)
 end
 
-(op::IntervalLinearAlgebra.NonLinearOettliPrager)(A, b, X=enclose(A, b)) = op(A, b, IntervalBox(X))
+(op::IntervalLinearAlgebra.NonLinearOettliPrager)(A, b, X=enclose(A, b)) = op(A, b, IntervalBox(X...))
 
 IntervalLinearAlgebra._default_precondition(_, ::NonLinearOettliPrager) = NoPrecondition()
 

src/IntervalLinearAlgebra.jl

@@ -1,6 +1,6 @@
 module IntervalLinearAlgebra
 
-using StaticArrays, Reexport
+using StaticArrays, Reexport, SetRounding
 using LinearAlgebra: checksquare
 
 import Base: +, -, *, /, \, ==,
@@ -9,6 +9,7 @@ import Base: +, -, *, /, \, ==,
 import CommonSolve: solve
 
 @reexport using LinearAlgebra, IntervalArithmetic
+using IntervalArithmetic.Symbols: var"..", ±
 
 const IA = IntervalArithmetic
 
@@ -79,6 +80,4 @@ function  __init__()
     end
 end
 
-set_multiplication_mode(config[:multiplication])
-
 end

src/classify.jl

@@ -12,16 +12,16 @@ For more details see section 4.6 of [[HOR19]](@ref).
 ```jldoctest
 julia> A = [2..4 -2..1; -1..2 2..4]
 2×2 Matrix{Interval{Float64}}:
-  [2, 4]  [-2, 1]
- [-1, 2]   [2, 4]
+  [2.0, 4.0]_com  [-2.0, 1.0]_com
+ [-1.0, 2.0]_com   [2.0, 4.0]_com
 
 julia> is_strongly_regular(A)
 true
 
 julia> A = [0..2 1..1;-1.. -1 0..2]
 2×2 Matrix{Interval{Float64}}:
-   [0, 2]  [1, 1]
- [-1, -1]  [0, 2]
+  [0.0, 2.0]_com   [1.0, 1.0]_com
+ [-1.0, -1.0]_com  [0.0, 2.0]_com
 
 julia> is_strongly_regular(A)
 false
@@ -49,16 +49,16 @@ For more details see section 4.4 of [[HOR19]](@ref).
 ```jldoctest
 julia> A = [2..4 -1..1; -1..1 2..4]
 2×2 Matrix{Interval{Float64}}:
-  [2, 4]  [-1, 1]
- [-1, 1]   [2, 4]
+  [2.0, 4.0]_com  [-1.0, 1.0]_com
+ [-1.0, 1.0]_com   [2.0, 4.0]_com
 
 julia> is_H_matrix(A)
 true
 
 julia> A = [2..4 -2..1; -1..2 2..4]
 2×2 Matrix{Interval{Float64}}:
-  [2, 4]  [-2, 1]
- [-1, 2]   [2, 4]
+  [2.0, 4.0]_com  [-2.0, 1.0]_com
+ [-1.0, 2.0]_com   [2.0, 4.0]_com
 
 julia> is_H_matrix(A)
 false
@@ -87,16 +87,16 @@ For more details see section 4.5 of [[HOR19]](@ref).
 ```jldoctest
 julia> A = [2..4 -1..1; -1..1 2..4]
 2×2 Matrix{Interval{Float64}}:
-  [2, 4]  [-1, 1]
- [-1, 1]   [2, 4]
+  [2.0, 4.0]_com  [-1.0, 1.0]_com
+ [-1.0, 1.0]_com   [2.0, 4.0]_com
 
 julia> is_strictly_diagonally_dominant(A)
 true
 
 julia> A = [2..4 -2..1; -1..2 2..4]
 2×2 Matrix{Interval{Float64}}:
-  [2, 4]  [-2, 1]
- [-1, 2]   [2, 4]
+  [2.0, 4.0]_com  [-2.0, 1.0]_com
+ [-1.0, 2.0]_com   [2.0, 4.0]_com
 
 julia> is_strictly_diagonally_dominant(A)
 false
@@ -107,7 +107,7 @@ function is_strictly_diagonally_dominant(A::AbstractMatrix{T}) where {T<:Interva
     m == n || return false
 
     @inbounds for i=1:m
-        sum_mag = sum(Interval(mag(A[i, k])) for k=1:n if k ≠ i)
+        sum_mag = sum(IA.interval(mag(A[i, k])) for k=1:n if k ≠ i)
         mig(A[i, i]) ≤ inf(sum_mag) && return false
 
     end
@@ -125,16 +125,16 @@ for all ``i≠j``. For more details see section 4.2 of [[HOR19]](@ref).
 ```jldoctest
 julia> A = [2..4 -2.. -1; -2.. -1 2..4]
 2×2 Matrix{Interval{Float64}}:
-   [2, 4]  [-2, -1]
- [-2, -1]    [2, 4]
+  [2.0, 4.0]_com   [-2.0, -1.0]_com
+ [-2.0, -1.0]_com   [2.0, 4.0]_com
 
 julia> is_Z_matrix(A)
 true
 
 julia> A = [2..4 -2..1; -1..2 2..4]
 2×2 Matrix{Interval{Float64}}:
-  [2, 4]  [-2, 1]
- [-1, 2]   [2, 4]
+  [2.0, 4.0]_com  [-2.0, 1.0]_com
+ [-1.0, 2.0]_com   [2.0, 4.0]_com
 
 julia> is_Z_matrix(A)
 false
@@ -166,16 +166,16 @@ non-negative inverse. For more details see section 4.2 of [[HOR19]](@ref).
 ```jldoctest
 julia> A = [2..2 -1..0; -1..0 2..2]
 2×2 Matrix{Interval{Float64}}:
-  [2, 2]  [-1, 0]
- [-1, 0]   [2, 2]
+  [2.0, 2.0]_com  [-1.0, 0.0]_com
+ [-1.0, 0.0]_com   [2.0, 2.0]_com
 
 julia> is_M_matrix(A)
 true
 
 julia> A = [2..4 -2..1; -1..2 2..4]
 2×2 Matrix{Interval{Float64}}:
-  [2, 4]  [-2, 1]
- [-1, 2]   [2, 4]
+  [2.0, 4.0]_com  [-2.0, 1.0]_com
+ [-1.0, 2.0]_com   [2.0, 4.0]_com
 
 julia> is_M_matrix(A)
 false

src/eigenvalues/interval_eigenvalues.jl

@@ -35,15 +35,15 @@ utilize normal floating point computations.
 ```jldoctest
 julia> A = [0 -1 -1; 2 -1.399.. -0.001 0; 1 0.5 -1]
 3×3 Matrix{Interval{Float64}}:
- [0, 0]  [-1, -1]                       [-1, -1]
- [2, 2]       [-1.39901, -0.000999999]    [0, 0]
- [1, 1]        [0.5, 0.5]               [-1, -1]
+ [0.0, 0.0]_com_NG  [-1.0, -1.0]_com_NG   [-1.0, -1.0]_com_NG
+ [2.0, 2.0]_com_NG  [-1.399, -0.001]_com   [0.0, 0.0]_com_NG
+ [1.0, 1.0]_com_NG   [0.5, 0.5]_com_NG    [-1.0, -1.0]_com_NG
 
 julia> eigenbox(A)
-[-1.90679, 0.970154] + [-2.51903, 2.51903]im
+[-1.90678, 0.970154]_com_NG + im*[-2.51903, 2.51903]_com_NG
 
 julia> eigenbox(A, Hertz())
-[-1.64732, 0.520456] + [-2.1112, 2.1112]im
+[-1.64731, 0.520455]_com_NG + im*[-2.11119, 2.11119]_com_NG
 ```
 """
 function eigenbox(A::Symmetric{Interval{T}, Matrix{Interval{T}}}, ::Rohn) where {T}
@@ -54,7 +54,7 @@ function eigenbox(A::Symmetric{Interval{T}, Matrix{Interval{T}}}, ::Rohn) where
     ρ = eigmax(AΔ)
     λmax = eigmax(Ac)
     λmin = eigmin(Ac)
-    return Interval(λmin - ρ, λmax + ρ)
+    return IA.interval(λmin - ρ, λmax + ρ)
 
 end
 
@@ -86,7 +86,7 @@ function eigenbox(A::Symmetric{Interval{T}, Matrix{Interval{T}}}, ::Hertz) where
         λmin = min(λmin, candmin)
         λmax = max(λmax, candmax)
     end
-    return IA.Interval(λmin, λmax)
+    return IA.interval(λmin, λmax)
 end
 
 function eigenbox(A::AbstractMatrix{Interval{T}},