Bump IntervalArithmetic to v0.22+

.github/workflows/test-pull-request.yml

@@ -25,7 +25,7 @@ jobs:
         arch:
           - x64
         include:
-          - version: '1.6'  # test on oldest supported version
+          - version: '1.10'  # test on oldest supported version
             arch: x64
             os: ubuntu-latest
     env:

Project.toml

@@ -1,6 +1,6 @@
 name = "LazySets"
 uuid = "b4f0291d-fe17-52bc-9479-3d1a343d9043"
-version = "5.1.0"
+version = "6.0.0"
 
 [deps]
 Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
@@ -22,7 +22,7 @@ StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
 Distributed = "<0.0.1, 1.6"
 ExprTools = "0.1"
 GLPK = "0.11 - 0.15, 1"
-IntervalArithmetic = "0.15 - 0.21, =0.21.2"  # v0.22 removed IntervalBox
+IntervalArithmetic = "0.22 - 1"
 JuMP = "0.21 - 0.23, 1"
 LinearAlgebra = "<0.0.1, 1.6"
 Random = "<0.0.1, 1.6"
@@ -33,4 +33,4 @@ Requires = "0.5, 1"
 SharedArrays = "<0.0.1, 1.6"
 SparseArrays = "<0.0.1, 1.6"
 StaticArraysCore = "1"
-julia = "1.6"
+julia = "1.10"

docs/Project.toml

@@ -8,6 +8,7 @@ DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244"
 ExponentialUtilities = "d4d017d3-3776-5f7e-afef-a10c40355c18"
 GR = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71"
 IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
+IntervalBoxes = "43d83c95-ebbb-40ec-8188-24586a1458ed"
 IntervalMatrices = "5c1f47dc-42dd-5697-8aaa-4d102d140ba9"
 LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
 MiniQhull = "978d7f02-9e05-4691-894f-ae31a51d76ca"
@@ -29,7 +30,8 @@ Documenter = "1"
 DocumenterCitations = "1.3"
 ExponentialUtilities = "1"
 GR = "0"
-IntervalArithmetic = "0.20 - 0.21, =0.21.2"  # v0.22 removed IntervalBox
+IntervalArithmetic = "0.22 - 1"
+IntervalBoxes = "0.2"
 IntervalMatrices = "0.8 - 0.12"
 LaTeXStrings = "1"
 MiniQhull = "0.4"
@@ -40,4 +42,4 @@ RecipesBase = "1"
 StaticArrays = "1"
 SymEngine = "0.7 - 0.13"
 Symbolics = "6.1"
-TaylorModels = "0.6 - 0.8"
+TaylorModels = "0.9"

docs/src/lib/sets/Interval.md

@@ -224,9 +224,3 @@ Inherited from [`AbstractHyperrectangle`](@ref):
 * [`genmat`](@ref genmat(::AbstractHyperrectangle))
 * [`is_interior_point`](@ref is_interior_point(::AbstractVector, ::AbstractHyperrectangle))
 * [`cartesian_product`](@ref cartesian_product(::AbstractHyperrectangle, ::AbstractHyperrectangle))
-
-Some additional functionality is available for `IntervalArithmetic.Interval`s:
-
-```@docs
-fast_interval_pow(::IA.Interval, ::Int)
-```

docs/src/man/tour.md

@@ -96,17 +96,6 @@ vector:
 H = convert(Hyperrectangle, C)
 ```
 
-IntervalArithmetic.jl represents multi-dimensional intervals with the type
-`IntervalBox`, to which we can also `convert`:
-
-```@example tour
-using IntervalArithmetic: IntervalBox
-
-Hbox = convert(IntervalBox, H)
-
-typeof(Hbox)
-```
-
 Hyperrectangles are a special subclass of
 [zonotopes](https://juliareach.github.io/LazySets.jl/dev/lib/sets/Zonotope/).
 

src/Approximations/box_approximation.jl

@@ -421,8 +421,8 @@ function load_taylormodels_box_approximation()
         box_approximation(vTM::Vector{<:TaylorModelN}) = _box_approximation_vTM(vTM)
 
         function _box_approximation_vTM(vTM)
-            B = IA.IntervalBox([evaluate(vTM[i], domain(p)) for (i, p) in enumerate(vTM)]...)
-            return convert(Hyperrectangle, B)
+            bounds = [evaluate(vTM[i], domain(p)) for (i, p) in enumerate(vTM)]
+            return Hyperrectangle(low=IA.inf.(bounds), high=IA.sup.(bounds))
         end
     end
 end  # quote / load_taylormodels_box_approximation

src/Approximations/init.jl

@@ -1,8 +1,8 @@
 function __init__()
     @require Expokit = "a1e7a1ef-7a5d-5822-a38c-be74e1bb89f4" (nothing,)
     @require ExponentialUtilities = "d4d017d3-3776-5f7e-afef-a10c40355c18" (nothing,)
-    @require IntervalMatrices = "5c1f47dc-42dd-5697-8aaa-4d102d140ba9" include("init_IntervalMatrices.jl")
     @require IntervalConstraintProgramming = "138f1668-1576-5ad7-91b9-7425abbf3153" include("init_IntervalConstraintProgramming.jl")
+    @require IntervalMatrices = "5c1f47dc-42dd-5697-8aaa-4d102d140ba9" include("init_IntervalMatrices.jl")
     @require Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9" include("init_Ipopt.jl")
     @require SetProg = "39881422-4b75-5582-a5c7-0feb14562a65" include("init_SetProg.jl")
     # @require StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c" include("init_StaticArraysCore.jl")

src/Approximations/init_IntervalConstraintProgramming.jl

@@ -1,2 +1 @@
-eval(load_paving_overapproximation())
 eval(load_overapproximate_ICP())

src/Approximations/overapproximate.jl

@@ -1,6 +1,5 @@
 using LazySets: block_to_dimension_indices,
                 substitute_blocks,
-                fast_interval_pow,
                 get_constrained_lowdimset
 
 """
@@ -369,7 +368,7 @@ function overapproximate(P::AbstractSparsePolynomialZonotope,
         return UnionSetArray([overapproximate(P, Zonotope)])
     end
 
-    dom = IA.IntervalBox(IA.interval(-1, 1), q)
+    dom = fill(IA.interval(-1, 1), q)
     cells = IA.mince(dom, isnothing(partition) ? nsdiv : partition)
     return UnionSetArray([overapproximate(P, Zonotope, c) for c in cells])
 end
@@ -677,51 +676,6 @@ function overapproximate(P::AbstractSparsePolynomialZonotope{N}, ::Type{<:VPolyt
     return VPolytope(vlist)
 end
 
-# function to be loaded by Requires
-function load_paving_overapproximation()
-    return quote
-        using .IntervalConstraintProgramming: Paving
-
-        """
-            overapproximate(p::Paving{L, N}, dirs::AbstractDirections{N, VN})
-                where {L, N, VN}
-
-        Overapproximate a Paving-type set representation with a polyhedron in constraint
-        representation.
-
-        ### Input
-
-        - `p`    -- paving
-        - `dirs` -- template directions
-
-        ### Output
-
-        An overapproximation of a paving using a polyhedron in constraint representation
-        (`HPolyhedron`) with constraints in direction `dirs`.
-
-        ### Algorithm
-
-        This function takes the union of the elements at the boundary of `p`, first
-        converted into hyperrectangles, and then calculates the support function of the
-        set along each  direction in dirs, to compute the `HPolyhedron` constraints.
-
-        This algorithm requires the IntervalConstraintProgramming package.
-        """
-        function overapproximate(p::Paving{L,N}, dirs::AbstractDirections{N,VN}) where {L,N,VN}
-            # enclose outer approximation
-            Uouter = UnionSetArray(convert.(Hyperrectangle, p.boundary))
-            constraints = [HalfSpace(d, ρ(d, Uouter)) for d in dirs]
-            return HPolyhedron(constraints)
-        end
-
-        # alias with HPolyhedron type as second argument
-        function overapproximate(p::Paving{L,N}, ::Type{<:HPolyhedron},
-                                 dirs::AbstractDirections{N,VN}) where {L,N,VN}
-            return overapproximate(p, dirs)
-        end
-    end
-end  # quote / load_paving_overapproximation
-
 """
 # Extended help
 

src/Approximations/overapproximate_zonotope.jl

@@ -247,7 +247,7 @@ end
 
 """
     overapproximate(P::AbstractSparsePolynomialZonotope, ::Type{<:Zonotope},
-                    dom::IntervalBox)
+                    dom::AbstractVector{<:IA.Interval})
 
 Overapproximate a sparse polynomial zonotope over the parameter domain `dom`
 with a zonotope.
@@ -264,9 +264,10 @@ with a zonotope.
 A zonotope.
 """
 function overapproximate(P::AbstractSparsePolynomialZonotope, ::Type{<:Zonotope},
-                         dom::IA.IntervalBox)
-    @assert dom ⊆ IA.IntervalBox(IA.interval(-1, 1), nparams(P)) "dom should " *
-                                                                 "be a subset of [-1, 1]^q"
+                         dom::AbstractVector{<:IA.Interval})
+    @assert length(dom) == nparams(P) &&
+            all(IA.issubset_interval(vi, IA.interval(-1, 1)) for vi in dom) "dom should be a " *
+                                                                            "subset of [-1, 1]^q"
 
     # handle dependent generators
     G = genmat_dep(P)
@@ -275,12 +276,11 @@ function overapproximate(P::AbstractSparsePolynomialZonotope, ::Type{<:Zonotope}
     Gnew = similar(G)
     @inbounds for (j, g) in enumerate(eachcol(G))
         # monomial value over the domain
-        # α = mapreduce(x -> _fast_interval_pow(x[1],  x[2]), *, zip(dom, E[:, i]))
         α = IA.interval(1, 1)
         for (i, vi) in enumerate(dom)
-            α *= fast_interval_pow(vi, E[i, j])
+            α *= vi^E[i, j]
         end
-        m, r = IA.midpoint_radius(α)
+        m, r = IA.midradius(α)
         cnew .+= m * g
         Gnew[:, j] .= r * g
     end
@@ -396,10 +396,10 @@ function load_taylormodels_overapproximation()
         [-0.5, 0.5]
 
         julia> x₀ = IA.interval(0.0) # expansion point
-        [0, 0]
+        [0.0, 0.0]
 
         julia> D = IA.interval(-3.0, 1.0)
-        [-3, 1]
+        [-3.0, 1.0]
 
         julia> p1 = Taylor1([2.0, 1.0], 2) # define a linear polynomial
          2.0 + 1.0 t + 𝒪(t³)
@@ -440,11 +440,10 @@ function load_taylormodels_overapproximation()
         julia> X = box_approximation(Z)
         Hyperrectangle{Float64, Vector{Float64}, Vector{Float64}}([1.0, -2.1], [2.5, 6.5])
 
-        julia> Y = evaluate(vTM[1], vTM[1].dom) × evaluate(vTM[2], vTM[2].dom)
-        [-1.5, 3.5] × [-8.60001, 4.40001]
-
-        julia> H = convert(Hyperrectangle, Y) # this IntervalBox is the same as X
-        Hyperrectangle{Float64, StaticArraysCore.SVector{2, Float64}, StaticArraysCore.SVector{2, Float64}}([1.0, -2.1000000000000005], [2.5, 6.500000000000001])
+        julia> Y = [evaluate(vTM[1], vTM[1].dom), evaluate(vTM[2], vTM[2].dom)]
+        2-element Vector{IntervalArithmetic.Interval{Float64}}:
+   │     [-1.50001, 3.50001]
+         [-8.60001, 4.40001]
         ```
         However, the zonotope returns better results if we want to approximate the
         Taylor model because it is not axis-aligned:
@@ -553,17 +552,21 @@ function load_taylormodels_overapproximation()
           1.0 x₁ + 𝒪(‖x‖⁹)
           1.0 x₂ + 𝒪(‖x‖⁹)
 
-        julia> x₀ = IA.IntervalBox(0..0, 2) # expansion point
-        [0, 0]²
+        julia> x₀ = fill(0..0, 2) # expansion point
+        2-element Vector{IntervalArithmetic.Interval{Float64}}:
+│        [0.0, 0.0]
+│        [0.0, 0.0]
 
         julia> Dx₁ = IA.interval(0.0, 3.0) # domain for x₁
-        [0, 3]
+        [0.0, 3.0]
 
         julia> Dx₂ = IA.interval(-1.0, 1.0) # domain for x₂
-        [-1, 1]
+        [-1.0, 1.0]
 
-        julia> D = Dx₁ × Dx₂ # take the Cartesian product of the domain on each variable
-        [0, 3] × [-1, 1]
+        julia> D = [Dx₁, Dx₂]
+        2-element Vector{IntervalArithmetic.Interval{Float64}}:
+│         [0.0, 3.0]
+│        [-1.0, 1.0]
 
         julia> r = IA.interval(-0.5, 0.5) # interval remainder
         [-0.5, 0.5]
@@ -629,13 +632,17 @@ function load_taylormodels_overapproximation()
 
                 # build the generators
                 α = mid(rem_nonlin)
-                c[i] = constant_term(pol_lin_norm) + α  # constant terms
-                G[i, 1:n] = get_linear_coeffs(pol_lin_norm)  # linear terms
+                cterm = constant_term(pol_lin_norm)
+                @assert IA.isthin(cterm) "unexpected interval value"
+                c[i] = mid(cterm) + α  # constant terms
+                lin_coeffs = get_linear_coeffs(pol_lin_norm)
+                @assert all(IA.isthin, lin_coeffs) "unexpected interval value"
+                G[i, 1:n] = mid.(lin_coeffs)  # linear terms
                 # interval generator
                 for j in (n + 1):(n + m)
                     G[i, j] = zero(N)
                 end
-                G[i, n + i] = abs(rem_nonlin.hi - α)
+                G[i, n + i] = abs(IA.sup(rem_nonlin) - α)
             end
 
             if remove_redundant_generators
@@ -1160,12 +1167,7 @@ function _overapproximate_zonotope_hyperplane(Z::AbstractZonotope, H::Hyperplane
 
     s = G' * a
     d = b - dot(a, c)
-    @static if VERSION >= v"1.7"
-        sT = s'
-    else
-        sT = s' .+ 0  # convert to Vector (`nullspace` fails for lazy transpose)
-    end
-    V0 = nullspace(sT)
+    V0 = nullspace(s')
 
     cs = s * d / dot(s, s)
     Gs = V0 * V0'
@@ -1194,35 +1196,40 @@ end
 # - the resulting zonotope is G * D + c
 function _overapproximate_zonotope_halfspace_ICP(Z::AbstractZonotope{N},
                                                  H::HalfSpace{N,SingleEntryVector{N}}) where {N}
+    require(@__MODULE__, :IntervalConstraintProgramming; fun_name="overapproximate")
+
     c = center(Z)
     G = genmat(Z)
     p = size(G, 2)
     d = H.a.i
 
     a = G[d, :]
     io = IOBuffer()
-    first = true
     v = H.a.v
     negate = v < zero(N)
     if negate
         v = -v
     end
+    vars = ""
+    first = true
     for (i, ai) in enumerate(a)
         if first
             first = false
         elseif ai >= 0
             write(io, "+")
+            vars *= ", "
         end
         write(io, "$(v * ai)*x$(i)")
+        vars *= "x$(i)"
     end
     if negate
         write(io, ">=$(-H.b + c[d])")
     else
         write(io, "<=$(H.b - c[d])")
     end
     e = Meta.parse(String(take!(io)))
-    X = IA.IntervalBox(IA.interval(-1, 1), p)
-    newD = _contract_zonotope_halfspace_ICP(e, X)
+    X = fill(IA.interval(-1, 1), p)
+    newD = _contract_zonotope_halfspace_ICP(e, X, vars)
     if isempty(newD)
         return EmptySet{N}(length(c))
     end
@@ -1232,15 +1239,22 @@ end
 function load_overapproximate_ICP()
     return quote
         import .IntervalConstraintProgramming as ICP
-
-        function _contract_zonotope_halfspace_ICP(e, X)
+        import .IntervalConstraintProgramming.IntervalBoxes as IB
+
+        function _contract_zonotope_halfspace_ICP(e, X, vars_string)
+            n = length(X)
+            prefix = "IntervalConstraintProgramming.Symbolics.@variables "
+            sym = Meta.parse(prefix * vars_string)
+            vars = eval(quote
+                            $sym
+                        end)
             separator = eval(quote
-                                 ICP.@constraint $e
+                                 ICP.Separator($e, $vars)
                              end)
             # smallest box containing all points in domain X satisfying constraint
-            # (`invokelatest` to avoid world-age issue; `Base.` for VERSION < v"1.9")
-            out, _ = Base.invokelatest(separator, X)
-            return out
+            # (`invokelatest` to avoid world-age issue)
+            boundary, _, _ = invokelatest(separator, IB.IntervalBox(X...))
+            return boundary
         end
     end
 end  # load_overapproximate_ICP()