Merge pull request #905 from JuliaReach/schillic/903

#903 - Fix/add `setops` option for `Forward`/`FirstOrder`

Author: schillic

src/Discretization/FirstOrderModule.jl

@@ -4,7 +4,8 @@
 module FirstOrderModule
 
 using ..DiscretizationModule
-using ..Exponentiation: _exp, BaseExp
+using ..Exponentiation: _exp, _alias, BaseExp
+using ..ApplySetops: _apply_setops
 using LinearAlgebra
 using MathematicalSystems
 using LazySets
@@ -36,18 +37,20 @@ The transformations are [LeGuernicG10](@cite):
 
 where ``B_ε`` is the input ball of radius ``ε`` centered in the origin.
 """
-struct FirstOrder{EM} <: AbstractApproximationModel
+struct FirstOrder{EM,SO} <: AbstractApproximationModel
     exp::EM
+    setops::SO
 end
 
 # convenience constructor
-function FirstOrder(; exp=BaseExp)
-    return FirstOrder(exp)
+function FirstOrder(; exp=BaseExp, setops=:lazy)
+    return FirstOrder(_alias(exp), _alias(setops))
 end
 
 function Base.show(io::IO, alg::FirstOrder)
     print(io, "`FirstOrder` approximation model with:\n")
     print(io, "    - exponentiation method: $(alg.exp)\n")
+    print(io, "    - set operations method: $(alg.setops)\n")
     return nothing
 end
 
@@ -70,6 +73,7 @@ function discretize(ivp::IVP{<:CLCS,<:AbstractZonotope}, δ, alg::FirstOrder)
     Xδ = Φ * X0
     α = (exp(norm_A * δ) - 1 - norm_A * δ) * norm(X0, Inf)
     Ω0 = ConvexHull(X0, Xδ + BallInf(zeros(n), α))
+    Ω0 = _apply_setops(Ω0, alg.setops)
 
     # create result
     X = stateset(ivp)
@@ -97,10 +101,12 @@ function discretize(ivp::IVP{<:CLCCS,<:AbstractZonotope}, δ, alg::FirstOrder)
     norm_U_over_A = isapproxzero(norm_A) ? 0 : norm(U, Inf) / norm_A
     α = factor * (norm(X0, Inf) + norm_U_over_A)
     Ω0 = ConvexHull(X0, Xδ + δ * U + BallInf(zeros(n), α))
+    Ω0 = _apply_setops(Ω0, alg.setops)
 
     # discretize inputs
     β = factor * norm_U_over_A
     V = δ * U + BallInf(zeros(n), β)
+    V = _apply_setops(V, alg.setops)
 
     # create result
     B = IdentityMultiple(one(eltype(A)), size(A, 1))

src/Discretization/ForwardModule.jl

@@ -87,8 +87,10 @@ function discretize(ivp::IVP{<:CLCS,<:LazySet}, δ, alg::Forward)
     Φcache = sum(A) == abs(sum(A)) ? Φ : nothing
     P2A_abs = Φ₂(A_abs, δ, alg.exp, alg.inv, Φcache)
     E₊ = sih(P2A_abs * sih((A * A) * X0, alg.sih), alg.sih)
+
     Ω0 = ConvexHull(X0, Φ * X0 ⊕ E₊)
     Ω0 = _apply_setops(Ω0, alg)
+
     X = stateset(ivp)
     Sdis = ConstrainedLinearDiscreteSystem(Φ, X)
     return InitialValueProblem(Sdis, Ω0)
@@ -115,6 +117,7 @@ function discretize(ivp::IVP{<:CLCS,Interval{N}}, δ, ::Forward) where {N}
     E⁺ = (Φ - one(N) - aδ) * convert(Interval, symmetric_interval_hull(X0)).dat
 
     Ω0 = Interval(hull(X0.dat, Φ * X0.dat + E⁺))
+
     X = stateset(ivp)
     # the system constructor creates a matrix
     Sdis = ConstrainedLinearDiscreteSystem(Φ, X)
@@ -143,6 +146,7 @@ function discretize(ivp::IVP{<:CLCCS,<:LazySet}, δ, alg::Forward)
 
     # discretize inputs
     V = δ * U ⊕ Eψ0
+    V = _apply_setops(V, alg.setops)
 
     Ω0 = ConvexHull(X0, Φ * X0 ⊕ V ⊕ E⁺)
     Ω0 = _apply_setops(Ω0, alg.setops)