Bump IntervalArithmetic to v0.22+

Project.toml

@@ -7,6 +7,7 @@ CarlemanLinearization = "4803f6b2-022a-4c1b-a771-522a3413ec86"
 CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"
 HybridSystems = "2207ec0c-686c-5054-b4d2-543502888820"
 IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
+IntervalBoxes = "43d83c95-ebbb-40ec-8188-24586a1458ed"
 IntervalMatrices = "5c1f47dc-42dd-5697-8aaa-4d102d140ba9"
 LazySets = "b4f0291d-fe17-52bc-9479-3d1a343d9043"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -23,13 +24,20 @@ TaylorIntegration = "92b13dbe-c966-51a2-8445-caca9f8a7d42"
 TaylorModels = "314ce334-5f6e-57ae-acf6-00b6e903104a"
 TaylorSeries = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea"
 
+[sources]
+CarlemanLinearization = {path = "/home/christian/.julia/dev/CarlemanLinearization"}
+IntervalBoxes = {path = "/home/christian/.julia/dev/IntervalBoxes"}
+LazySets = {path = "/home/christian/.julia/dev/LazySets"}
+MathematicalSystems = {path = "/home/christian/.julia/dev/MathematicalSystems"}
+
 [compat]
 CarlemanLinearization = "0.3 - 0.4"
 CommonSolve = "0.2"
 HybridSystems = "0.4"
-IntervalArithmetic = "0.16 - 0.20, =0.20.9"  # new versions require updates and are incompatible with dependencies
+IntervalArithmetic = "0.22 - 1"
+IntervalBoxes = "0.2"
 IntervalMatrices = "0.11 - 0.12"
-LazySets = "2.14, 3, 4, 5"
+LazySets = "6"
 LinearAlgebra = "<0.0.1, 1.6"
 MathematicalSystems = "0.14.1"
 Parameters = "0.10 - 0.12"
@@ -40,7 +48,7 @@ Reexport = "0.2, 1"
 Requires = "0.5, 1"
 SparseArrays = "<0.0.1, 1.6"
 StaticArrays = "0.12, 1"
-TaylorIntegration = "0.15 - 0.16"
-TaylorModels = "0.7 - 0.8"
-TaylorSeries = "0.17 - 0.18"
-julia = "1.6"
+TaylorIntegration = "0.17"
+TaylorModels = "0.9"
+TaylorSeries = "0.20.6"
+julia = "1.10"

docs/Project.toml

@@ -15,15 +15,18 @@ Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 ReachabilityBase = "379f33d0-9447-4353-bd03-d664070e549f"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 
+[sources]
+LazySets = {path = "/home/christian/.julia/dev/LazySets"}
+
 [compat]
 DisplayAs = "0.1"
 Documenter = "1"
 DocumenterCitations = "1.3"
 ExponentialUtilities = "1"
-IntervalArithmetic = "=0.20.9"  # new versions require updates and are incompatible with dependencies
+IntervalArithmetic = "0.22 - 1"
 JLD2 = "0.4 - 0.6"
 LaTeXStrings = "1"
-LazySets = "2.14, 3, 4, 5"
+LazySets = "6"
 Literate = "2"
 MathematicalSystems = "0.14.1"
 Optim = "0.15 - 0.22, 1"

docs/make.jl

@@ -224,13 +224,13 @@ makedocs(; sitename="ReachabilityAnalysis.jl",
                 "Further examples" => [
                                        #
                                        "Overview"               => "man/examples_overview.md",
-                                       "Duffing oscillator"     => "generated_examples/DuffingOscillator.md",
+                                    #    "Duffing oscillator"     => "generated_examples/DuffingOscillator.md",  # TODO temporarily deactivated
                                        "Laub-Loomis"            => "generated_examples/LaubLoomis.md",
                                        "Production-Destruction" => "generated_examples/ProductionDestruction.md",
                                        "Brusselator"            => "generated_examples/Brusselator.md",
                                        "Quadrotor"              => "generated_examples/Quadrotor.md",
                                        "Epidemic (SEIR) model"  => "generated_examples/SEIR.md",
-                                       "Spacecraft"             => "generated_examples/Spacecraft.md"
+                                    #    "Spacecraft"             => "generated_examples/Spacecraft.md"  # TODO temporarily deactivated
                                        #
                                        ],
                 "API Reference" => [

docs/src/man/faq.md

@@ -139,7 +139,7 @@ with the `trajectories` keyword argument.
 using ReachabilityAnalysis, OrdinaryDiffEq
 
 # formulate initial-value problem
-prob = @ivp(x' = 1.01x, x(0) ∈ 0 .. 0.5)
+prob = @ivp(x' = 1.01x, x(0) ∈ Interval(0, 0.5))
 
 # solve the flowpipe using a default algorithm, and also compute trajectories
 sol = solve(prob, tspan=(0.0, 1.0), ensemble=true, trajectories=250)
@@ -223,7 +223,7 @@ boxes intersecting the time point `3.0` decrease by a factor 2.5x.
 ```@example cosine
 fig = plot(f(0.1)(3.0), vars=(0, 1), xlab="time", ylab="x(t)", lab="ΔT=0.1", color=:lightblue)
 
-I(Δt, t) = -ρ([-1.0, 0.0], f(Δt)(t)) .. ρ([1.0, 0.0], f(Δt)(t)) |> Interval
+I(Δt, t) = Interval(-ρ([-1.0, 0.0], f(Δt)(t)), ρ([1.0, 0.0], f(Δt)(t)))
 
 I01 = I(0.1, 3.0)
 plot!(y -> max(I01), xlims=(2.9, 3.1), lw=3.0, style=:dash, color=:lightblue, lab="Δx = $(I01.dat)")
@@ -291,20 +291,6 @@ Although it is in principle possible to  ODE solvers for
 The section [Some common gotchas](@ref) of the user manual details do's and dont's
 for the `@taylorize` macro to speedup reachability computations using Taylor models.
 
-### A note on interval types
-
-When using intervals as set representation, `ReachabilityAnalysis.jl` relies on
-rigorous floating-point arithmetic implemented in pure Julia in the library [IntervalArithmetic.jl](https://github.com/JuliaIntervals/IntervalArithmetic.jl) (we often use `const IA = IntervalArithmetic` as an abbreviation).
-The main struct defined in the library is `IA.Interval` (and the corresponding
-multi-dimensional interval is `IA.IntervalBox`). Internally, the set `LazySets.Interval`
-is **wrapper-type** of `IA.Interval` and these two types should not be confused,
-although our user APIs extensively use [duck typing](https://en.wikipedia.org/wiki/Duck_typing),
-in the sense that `x(0) ∈ 0 .. 1` (`IA.Interval` type) and `x(0) ∈ Interval(0, 1)` are valid.
-
-On a technical level, the reason to have `LazySets.Interval` as a wrapper type
-of `IA.Interval` is that Julia doesn't allow multiple inheritance, but it was a design
-choice that intervals should belong to the `LazySets` type hierarchy.
-
 ## References
 
 [^v1]: Version 1.0 of the language was released in August 2018; see [this blog post](https://julialang.org/blog/2018/08/one-point-zero/).

docs/src/man/introduction.md

@@ -15,7 +15,7 @@ with initial condition ``x(0) ∈ [0.45, 0.55]``. We can solve the problem as fo
 using ReachabilityAnalysis, Plots
 
 # define the initial-value problem
-prob = @ivp(x' = -x, x(0) ∈ 0.45 .. 0.55)
+prob = @ivp(x' = -x, x(0) ∈ Interval(0.45, 0.55))
 
 # solve it
 sol = solve(prob, T=4.0)

docs/src/tutorials/linear_methods/discrete_time.md

@@ -229,14 +229,14 @@ We can also pick, or filter, those reach-sets whose intersection is non-empty wi
 a given time interval.
 
 ```@example discrete_propagation
-F(2.5 .. 3.0) # returns a view
+F((2.5, 3.0)) # returns a view
 ```
 
 We can also filter by a condition on the time span using `tstart` and `tend`:
 
 ```@example discrete_propagation
 # get all reach-sets from time t = 4 onwards
-aux = F(4 .. tend(F));
+aux = F((4, tend(F)));
 
 length(aux)
 ```

docs/src/tutorials/taylor_methods/introduction.md

@@ -1,5 +1,4 @@
 ```@meta
-DocTestSetup = :(using ReachabilityAnalysis)
 CurrentModule = ReachabilityAnalysis
 ```
 
@@ -20,29 +19,29 @@ x'(t) = -x(t) ~ \sin(t),\qquad t ≥ 0.
 Standard integration schemes fail to produce helpful solutions if the initial state is an interval. We illustrate this point
 by solving the given differential equation with the `Tsit5` algorithm from [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl) suite.
 
-```@example nonlinear_univariate
+```@example nonlinear_univariate_bland
 using OrdinaryDiffEq, IntervalArithmetic
 
 # initial condition
-x₀ = [-1 .. 1]
+X0 = [interval(-1, 1)]
 
 # define the problem
 function f(dx, x, p, t)
     dx[1] = -x[1] * sin(t)
 end
 
 # pass to solvers
-prob = ODEProblem(f, x₀, (0.0, 2.0))
+prob = ODEProblem(f, X0, (0.0, 2.0))
 sol = solve(prob, Tsit5(), adaptive=false, dt=0.05, reltol=1e-6)
 nothing # hide
 ```
-There is no plot recipe readily available so we create it by hand using [`LazySets.jl`](https://github.com/JuliaReach/LazySets.jl).
 
-```@example nonlinear_univariate
-using LazySets, Plots
-using LazySets: Interval
+Now we plot the result.
+
+```@example nonlinear_univariate_bland
+using Plots
 
-out = [Interval(sol.t[i]) × Interval(sol.u[i][1]) for i in 1:20]
+out = [[interval(sol.t[i]), interval(sol.u[i][1])] for i in 1:20]
 
 fig = plot(out, xlab="t", ylab="x(t)", lw=3.0, alpha=1., c=:black, marker=:none, lab="", title="Standard integrator with an interval initial condition")
 
@@ -63,7 +62,7 @@ function f(dx, x, p, t)
 end
 
 # define the set of initial states
-X0 = -1 .. 1
+X0 = Interval(-1, 1)
 
 # define the initial-value problem
 prob = @ivp(x' = f(x), x(0) ∈ X0, dim=1)
@@ -139,7 +138,7 @@ Filtering by the time span also works with time intervals; in the following exam
 `1.0 .. 3.0` is non-empty:
 
 ```@example nonlinear_univariate
-length(sol(1.0 .. 3.0))
+length(sol((1.0, 3.0)))
 ```
 
 On the other hand, evaluating over a given time point or time interval can be achieved using the `evaluate` function:

examples/Brusselator/Brusselator.jl

@@ -49,7 +49,7 @@ end
 # The initial set is ``x ∈ [0.8, 1]``, ``y ∈ [0, 0.2]``. The time horizon is
 # ``18``. These settings are taken from [ChenAS13](@citet).
 
-U₀ = (0.8 .. 1.0) × (0.0 .. 0.2)
+U₀ = Hyperrectangle(low=[0.8, 0], high=[1, 0.2])
 prob = @ivp(u' = brusselator!(u), u(0) ∈ U₀, dim:2)
 
 T = 18.0;

examples/DuffingOscillator/DuffingOscillator.jl

@@ -21,7 +21,7 @@ const ω = 1.2
 
     f = γ * cos(ω * t)
 
-    du[1] = u[2]
+    du[1] = v
     du[2] = -α * x - δ * v - β * x^3 + f
     return du
 end

examples/EMBrake/EMBrake.jl

@@ -34,6 +34,7 @@
 #    constant over time: ``u(0) ∈ \mathcal{U}``, ``\dot{u}(t) = 0.``
 
 using ReachabilityAnalysis, SparseArrays
+using IntervalArithmetic: interval
 
 # ## Common functionality between model variants
 
@@ -93,7 +94,7 @@ function embrake_pv_1(; Tsample=1.E-4, ζ=1e-6, Δ=3.0, x0=0.05)
     K = 0.02
     drot = 0.1
     i = 113.1167
-    p = 504.0 + (-Δ .. Δ)
+    p = 504.0 + interval(-Δ, Δ)
 
     ## state variables: [I, x, xe, xc]
     A = IntervalMatrix([-p 0 KP/L KI/L;
@@ -111,7 +112,7 @@ end;
 
 function embrake_pv_2(; Tsample=1.E-4, ζ=1e-6, x0=0.05, χ=5.0)
     ## model's constants
-    Δ = -χ / 100 .. χ / 100
+    Δ = interval(-χ / 100, χ / 100)
     L = 1.e-3 * (1 + Δ)
     KP = 10000.0 * (1 + Δ)
     KI = 1000.0 * (1 + Δ)

examples/Lorenz/Lorenz.jl

@@ -69,19 +69,19 @@ fig = plot(solz; vars=(0, 1), xlab="t", ylab="x")  #!jl
 # It is apparent by inspection that variable ``x(t)`` does not exceed ``20`` in
 # the computed time span:
 
-fig = plot(solz(0.0 .. 1.5); vars=(0, 1), xlab="t", ylab="x", lw=0)  #!jl
+fig = plot(solz((0.0, 1.5)); vars=(0, 1), xlab="t", ylab="x", lw=0)  #!jl
 plot!(fig, x -> 20.0; c=:red, xlims=(0.0, 1.5), lab="")  #!jl
 
 #!jl DisplayAs.Text(DisplayAs.PNG(fig))  #hide
 
 # We can prove that this is the case by evaluating the support function of the
 # flowpipe along direction ``[1, 0, 0]``:
 
-@assert ρ([1.0, 0, 0], solz(0 .. 1.5)) < 20 "the property should be proven"
+@assert ρ([1.0, 0, 0], solz((0, 1.5))) < 20 "the property should be proven"
 
 #-
 
-ρ([1.0, 0, 0], solz(0 .. 1.5))  #!jl
+ρ([1.0, 0, 0], solz((0, 1.5)))  #!jl
 
 # In a similar fashion, we can compute extremal values of variable ``y(t)``:
 

examples/LotkaVolterra/LotkaVolterra.jl

@@ -84,8 +84,8 @@ end
 
 # ## Specification
 
-p_int = (0.99 .. 1.01) × (0.99 .. 1.01) × (2.99 .. 3.01) × (0.99 .. 1.01) × (0.099 .. 0.101)
-U0 = cartesian_product(Singleton([1.0, 1.0]), convert(Hyperrectangle, p_int))
+p_int = Hyperrectangle(low=[0.99, 0.99, 2.99, 0.99, 0.099], high=[1.01, 1.01, 3.01, 1.01, 0.101])
+U0 = cartesian_product(Singleton([1.0, 1.0]), p_int)
 prob = @ivp(u' = lotkavolterra_parametric!(u), dim:7, u(0) ∈ U0);
 
 # ## Analysis

examples/OpAmp/OpAmp.jl

@@ -110,7 +110,7 @@ fig = plot(sol_nondet; vars=(0, 1), xlab=L"t", ylab=L"e_{out}",  #!jl
 # we solve for three different initial conditions of increasing width, for a
 # shorter time horizon.
 
-Δt = 0 .. 0.04  #!jl
+Δt = (0, 0.04)  #!jl
 
 fig = plot(; xlab=L"t", ylab=L"e_{out}(t)", title="Solution for nondeterministic input")  #!jl
 

examples/ProductionDestruction/ProductionDestruction.jl

@@ -102,8 +102,8 @@ function prod_dest_verif(sol; T=100.0, target=10.0)
 
     ## check that the target belongs to the Minkowski sum of the projections
     ## onto each coordinate
-    B = convert(IntervalBox, X) # get a product-of-intervals representation
-    contains_target = target ∈ sum(B)
+    S = concretize(MinkowskiSumArray([project(set(X), [i]) for i in 1:3]))
+    contains_target = [target] ∈ S
 
     return nonnegative && contains_target, vol
 end;

examples/Quadrotor/Quadrotor.jl

@@ -166,12 +166,12 @@ const v3 = SingleEntryVector(3, 12, 1.0)
 
     ## Condition: b1 = (x[3] < 1.4) for all time
     unsafe1 = HalfSpace(-v3, -1.4) # unsafe: -x3 <= -1.4
-    b1 = all([isdisjoint(unsafe1, set(R)) for R in sol(0.0 .. tf)])
+    b1 = all([isdisjoint(unsafe1, set(R)) for R in sol((0.0, tf))])
     #b1 = ρ(v3, sol) < 1.4
 
     ## Condition: x[3] > 0.9 for t ≥ 1.0
     unsafe2 = HalfSpace(v3, 0.9) # unsafe: x3 <= 0.9
-    b2 = all([isdisjoint(unsafe2, set(R)) for R in sol(1.0 .. tf)])
+    b2 = all([isdisjoint(unsafe2, set(R)) for R in sol((1.0, tf))])
 
     ## Condition: x[3] ⊆ Interval(0.98, 1.02) for t = 5.0
     b3 = set(project(sol[end]; vars=(3))) ⊆ Interval(0.98, 1.02)

examples/SEIR/SEIR.jl

@@ -54,11 +54,7 @@ end
 
 E₀ = 1e-4
 x₀ = [1 - E₀, E₀, 0, 0]
-α = 0.2 ± 0.01
-β = 1.0 ± 0.0
-γ = 0.5 ± 0.01
-p = [α, β, γ]
-X0 = IntervalBox(vcat(x₀, p));
+X0 = Hyperrectangle(vcat(x₀, [0.2, 1, 0.5]), vcat(zeros(4), [0.01, 0, 0.01]))
 prob = @ivp(x' = seir!(x), dim:7, x(0) ∈ X0);
 
 # ## Analysis

examples/Spacecraft/Spacecraft.jl

@@ -320,7 +320,7 @@ function solve_spacecraft(prob; k=25, s=missing)
                   intersection_method=TemplateHullIntersection(BoxDirections(5)),
                   clustering_method=LazyClustering(1),
                   disjointness_method=BoxEnclosure())
-    sol12jump = overapproximate(sol12[2](120 .. 150), Zonotope)
+    sol12jump = overapproximate(sol12[2]((120, 150)), Zonotope)
     t0 = tstart(sol12jump[1])
     sol12jump_c = cluster(sol12jump, 1:length(sol12jump), BoxClustering(k, s))
 

src/Algorithms/CARLIN/kronecker.jl

@@ -130,7 +130,7 @@ Return a hyperrectangle with the interval powers `[x, x^2, …, x^pow]`.
 A hyperrectangle such that the `i`-th dimension is the interval `x^i`.
 """
 function kron_pow_stack(x::IA.Interval, pow::Int)
-    return convert(Hyperrectangle, IntervalBox([kron_pow(x, i) for i in 1:pow]))
+    return convert(Hyperrectangle, IntervalBox([kron_pow(x, i) for i in 1:pow]...))
 end
 
 """

src/Algorithms/FLOWSTAR/post.jl

@@ -50,7 +50,16 @@ function post(alg::FLOWSTAR{ST,OT,PT,IT}, ivp::IVP{<:AbstractContinuousSystem},
 
         # initial set
         X0 = initial_state(ivp)
-        dom = isa(X0, IntervalBox) ? X0 : convert(IntervalBox, overapproximate(X0, Hyperrectangle))
+        if X0 isa AbstractVector{<:IA.Interval}
+            dom = X0
+        elseif X0 isa IntervalBox
+            dom = [X0[i] for i in 1:length(X0)]
+        elseif X0 isa LazySet
+            IB = convert(IntervalBox, overapproximate(X0, Hyperrectangle))
+            dom = [IB[i] for i in 1:length(IB)]
+        else
+            throw(ArgumentError("invalid `X0` of type $(typeof(X0))"))
+        end
 
         name = "flowstar" # TODO parse name of f
 
@@ -85,7 +94,7 @@ function post(alg::FLOWSTAR{ST,OT,PT,IT}, ivp::IVP{<:AbstractContinuousSystem},
     counter = t0
     for Fi in flow
         dt = domain(first(Fi))
-        δt = TimeInterval(counter + dt)
+        δt = TimeIntervalC(counter + dt)
         counter += sup(dt)
         Ri = TaylorModelReachSet(Fi, δt + Δt0)
         push!(F, Ri)

src/Algorithms/QINT/reach_homog.jl

@@ -45,7 +45,7 @@ function reach_homog_QINT(; a, b, c, # right-hand side: f(x) = ax^2 + bx + c
 
         # solve linear reachability
         prob = @ivp(x' = α * x + β, x(0) ∈ X0)
-        sol = post(INT(; δ=δ), prob, IA.interval(0.0, Δ); Δt0=Δti)
+        sol = post(INT(; δ=δ), prob, TimeIntervalC(0.0, Δ); Δt0=Δti)
 
         # compute admissible linearization error
         kθ = 1 / (exp(α * Δ) - 1) * α * θ * Δ

src/Algorithms/TMJets/TMJets21a/TMJets21a.jl

@@ -35,7 +35,7 @@ and `TaylorSeries.jl` is used to work with truncated Taylor series.
     abstol::N = DEFAULT_ABS_TOL_TMJETS
     maxsteps::Int = DEFAULT_MAX_STEPS_TMJETS
     adaptive::Bool = true
-    minabstol::N = Float64(TM._DEF_MINABSTOL)
+    minabstol::N = Float64(_DEF_MINABSTOL)
     absorb::Bool = false
     disjointness::DM = ZonotopeEnclosure()
 end