Merge pull request #887 from JuliaReach/schillic/discrete_post

Add discrete post for linear algorithms
JuliaReach · Jan 9, 2025 · 762beb4 · 762beb4
2 parents 0e4056f + bf3ccd2
commit 762beb4
Show file tree

Hide file tree

Showing 14 changed files with 287 additions and 144 deletions.
diff --git a/src/Algorithms/A20/post.jl b/src/Algorithms/A20/post.jl
@@ -1,12 +1,11 @@
-function post(alg::A20{N}, ivp::IVP{<:AbstractContinuousSystem}, tspan;
-              Δt0::TimeInterval=zeroI, kwargs...) where {N}
-    @unpack δ, max_order = alg
+# continuous post
+function post(alg::A20, ivp::IVP{<:AbstractContinuousSystem}, tspan;
+              Δt0::TimeInterval=zeroI, kwargs...)
+    δ = alg.δ
 
     NSTEPS = get(kwargs, :NSTEPS, compute_nsteps(δ, tspan))
 
     # normalize system to canonical form
-    # x' = Ax, x in X, or
-    # x' = Ax + u, x in X, u in U
     ivp_norm = _normalize(ivp)
 
     # homogenize the initial-value problem
@@ -15,13 +14,30 @@ function post(alg::A20{N}, ivp::IVP{<:AbstractContinuousSystem}, tspan;
     end
 
     # discretize system
-    ivp_discr = discretize(ivp_norm, δ, approx_model)
-    Φ = state_matrix(ivp_discr)
-    Ω0 = initial_state(ivp_discr)
-    X = stateset(ivp_discr)
+    ivp_discr = discretize(ivp_norm, δ, alg.approx_model)
+
+    return post(alg, ivp_discr, NSTEPS; Δt0=Δt0, kwargs...)
+end
+
+# discrete post
+function post(alg::A20{N}, ivp::IVP{<:AbstractDiscreteSystem}, NSTEPS=nothing;
+              Δt0::TimeInterval=zeroI, kwargs...) where {N}
+    @unpack δ, max_order = alg
+
+    if isnothing(NSTEPS)
+        if haskey(kwargs, :NSTEPS)
+            NSTEPS = kwargs[:NSTEPS]
+        else
+            throw(ArgumentError("`NSTEPS` not specified"))
+        end
+    end
+
+    Φ = state_matrix(ivp)
+    Ω0 = initial_state(ivp)
+    X = stateset(ivp)
 
     # true <=> there is no input, i.e. the system is of the form x' = Ax, x ∈ X
-    got_homogeneous = !hasinput(ivp_discr)
+    got_homogeneous = !hasinput(ivp)
 
     # this algorithm requires Ω0 to be a zonotope
     Ω0 = _convert_or_overapproximate(Zonotope, Ω0)
@@ -53,7 +69,7 @@ function post(alg::A20{N}, ivp::IVP{<:AbstractContinuousSystem}, tspan;
                             disjointness_method)
     else
         # TODO: implement preallocate option for this scenario
-        U = inputset(ivp_discr)
+        U = inputset(ivp)
         @assert isa(U, LazySet) "expected input of type `<:LazySet`, but got $(typeof(U))"
         U = _convert_or_overapproximate(Zonotope, U)
         reach_inhomog_GLGM06!(F, Ω0, Φ, NSTEPS, δ, max_order, X, U, reduction_method, Δt0,

diff --git a/src/Algorithms/ASB07/post.jl b/src/Algorithms/ASB07/post.jl
@@ -1,22 +1,38 @@
-function post(alg::ASB07{N}, ivp::IVP{<:AbstractContinuousSystem}, tspan;
-              Δt0::TimeInterval=zeroI, kwargs...) where {N}
-    @unpack δ, approx_model, max_order, reduction_method, static, recursive, dim, ngens = alg
+# continuous post
+function post(alg::ASB07, ivp::IVP{<:AbstractContinuousSystem}, tspan;
+              Δt0::TimeInterval=zeroI, kwargs...)
+    δ = alg.δ
 
     NSTEPS = get(kwargs, :NSTEPS, compute_nsteps(δ, tspan))
 
     # normalize system to canonical form
-    # x' = Ax, x in X, or
-    # x' = Ax + u, x in X, u in U
     ivp_norm = _normalize(ivp)
 
     # discretize system
-    ivp_discr = discretize(ivp_norm, δ, approx_model)
-    Φ = state_matrix(ivp_discr)
-    Ω0 = initial_state(ivp_discr)
-    X = stateset(ivp_discr)
+    ivp_discr = discretize(ivp_norm, δ, alg.approx_model)
+
+    return post(alg, ivp_discr, NSTEPS; Δt0=Δt0, kwargs...)
+end
+
+# discrete post
+function post(alg::ASB07{N}, ivp::IVP{<:AbstractDiscreteSystem}, NSTEPS=nothing;
+              Δt0::TimeInterval=zeroI, kwargs...) where {N}
+    @unpack δ, approx_model, max_order, reduction_method, static, recursive, dim, ngens = alg
+
+    if isnothing(NSTEPS)
+        if haskey(kwargs, :NSTEPS)
+            NSTEPS = kwargs[:NSTEPS]
+        else
+            throw(ArgumentError("`NSTEPS` not specified"))
+        end
+    end
+
+    Φ = state_matrix(ivp)
+    Ω0 = initial_state(ivp)
+    X = stateset(ivp)
 
     # true <=> there is no input, i.e. the system is of the form x' = Ax, x ∈ X
-    got_homogeneous = !hasinput(ivp_discr)
+    got_homogeneous = !hasinput(ivp)
 
     # this algorithm requires Ω0 to be a zonotope
     Ω0 = _convert_or_overapproximate(Zonotope, Ω0)
@@ -34,7 +50,7 @@ function post(alg::ASB07{N}, ivp::IVP{<:AbstractContinuousSystem}, tspan;
     if got_homogeneous
         reach_homog_ASB07!(F, Ω0, Φ, NSTEPS, δ, max_order, X, recursive, reduction_method, Δt0)
     else
-        U = inputset(ivp_discr)
+        U = inputset(ivp)
         @assert isa(U, LazySet) "expected input of type `<:LazySet`, but got $(typeof(U))"
         U = _convert_or_overapproximate(Zonotope, U)
         reach_inhomog_ASB07!(F, Ω0, Φ, NSTEPS, δ, max_order, X, U, recursive, reduction_method, Δt0)

diff --git a/src/Algorithms/BFFPSV18/post.jl b/src/Algorithms/BFFPSV18/post.jl
@@ -1,13 +1,11 @@
-function post(alg::BFFPSV18{N,ST}, ivp::IVP{<:AbstractContinuousSystem}, tspan;
-              Δt0::TimeInterval=zeroI, kwargs...) where {N,ST}
-    @unpack δ, approx_model, vars, block_indices,
-    row_blocks, column_blocks = alg
+# continuous post
+function post(alg::BFFPSV18, ivp::IVP{<:AbstractContinuousSystem}, tspan;
+              Δt0::TimeInterval=zeroI, kwargs...)
+    δ = alg.δ
 
     NSTEPS = get(kwargs, :NSTEPS, compute_nsteps(δ, tspan))
 
     # normalize system to canonical form
-    # x' = Ax, x in X, or
-    # x' = Ax + u, x in X, u in U
     ivp_norm = _normalize(ivp)
 
     # homogenize the initial-value problem
@@ -16,19 +14,35 @@ function post(alg::BFFPSV18{N,ST}, ivp::IVP{<:AbstractContinuousSystem}, tspan;
     end
 
     # discretize system
-    ivp_discr = discretize(ivp_norm, δ, approx_model)
-    Φ = state_matrix(ivp_discr)
-    Ω0 = initial_state(ivp_discr)
-    X = stateset(ivp_discr)
+    ivp_discr = discretize(ivp_norm, δ, alg.approx_model)
+
+    return post(alg, ivp_discr, NSTEPS; Δt0=Δt0, kwargs...)
+end
+
+# discrete post
+function post(alg::BFFPSV18{N,ST}, ivp::IVP{<:AbstractDiscreteSystem}, NSTEPS=nothing;
+              Δt0::TimeInterval=zeroI, kwargs...) where {N,ST}
+    @unpack δ, approx_model, vars, block_indices,
+    row_blocks, column_blocks = alg
+
+    if isnothing(NSTEPS)
+        if haskey(kwargs, :NSTEPS)
+            NSTEPS = kwargs[:NSTEPS]
+        else
+            throw(ArgumentError("`NSTEPS` not specified"))
+        end
+    end
+
+    Φ = state_matrix(ivp)
+    Ω0 = initial_state(ivp)
+    X = stateset(ivp)
 
     # true <=> there is no input, i.e. the system is of the form x' = Ax, x ∈ X
-    got_homogeneous = !hasinput(ivp_discr)
+    got_homogeneous = !hasinput(ivp)
 
     # decompose the initial states into a cartesian product
     # TODO add option to do the lazy decomposition
     Xhat0 = _decompose(Ω0, column_blocks, ST)
-    Φ = state_matrix(ivp_discr)
-    X = stateset(ivp_discr) # invariant
 
     # force using sparse type for the matrix exponential
     if alg.sparse
@@ -51,7 +65,7 @@ function post(alg::BFFPSV18{N,ST}, ivp::IVP{<:AbstractContinuousSystem}, tspan;
                               row_blocks, column_blocks, Δt0, viewval)
 
     else
-        U = inputset(ivp_discr)
+        U = inputset(ivp)
         @assert isa(U, LazySet) "expected input of type `<:LazySet`, but got $(typeof(U))"
         reach_inhomog_BFFPSV18!(F, Xhat0, Φ, NSTEPS, δ, X, U, ST,
                                 vars, block_indices,

diff --git a/src/Algorithms/BOX/post.jl b/src/Algorithms/BOX/post.jl
@@ -1,6 +1,7 @@
-function post(alg::BOX{N}, ivp::IVP{<:AbstractContinuousSystem}, tspan;
-              Δt0::TimeInterval=zeroI, kwargs...) where {N}
-    @unpack δ, approx_model, static, dim, recursive = alg
+# continuous post
+function post(alg::BOX, ivp::IVP{<:AbstractContinuousSystem}, tspan;
+              Δt0::TimeInterval=zeroI, kwargs...)
+    δ = alg.δ
 
     NSTEPS = get(kwargs, :NSTEPS, compute_nsteps(δ, tspan))
 
@@ -13,13 +14,30 @@ function post(alg::BOX{N}, ivp::IVP{<:AbstractContinuousSystem}, tspan;
     end
 
     # discretize system
-    ivp_discr = discretize(ivp_norm, δ, approx_model)
-    Φ = state_matrix(ivp_discr)
-    Ω0 = initial_state(ivp_discr)
-    X = stateset(ivp_discr)
+    ivp_discr = discretize(ivp_norm, δ, alg.approx_model)
+
+    return post(alg, ivp_discr, NSTEPS; Δt0=Δt0, kwargs...)
+end
+
+# discrete post
+function post(alg::BOX{N}, ivp::IVP{<:AbstractDiscreteSystem}, NSTEPS=nothing;
+              Δt0::TimeInterval=zeroI, kwargs...) where {N}
+    @unpack δ, approx_model, static, dim, recursive = alg
+
+    if isnothing(NSTEPS)
+        if haskey(kwargs, :NSTEPS)
+            NSTEPS = kwargs[:NSTEPS]
+        else
+            throw(ArgumentError("`NSTEPS` not specified"))
+        end
+    end
+
+    Φ = state_matrix(ivp)
+    Ω0 = initial_state(ivp)
+    X = stateset(ivp)
 
     # true <=> there is no input, i.e. the system is of the form x' = Ax, x ∈ X
-    got_homogeneous = !hasinput(ivp_discr)
+    got_homogeneous = !hasinput(ivp)
 
     # this algorithm requires Ω0 to be hyperrectangle
     Ω0 = _overapproximate(Ω0, Hyperrectangle)
@@ -36,7 +54,7 @@ function post(alg::BOX{N}, ivp::IVP{<:AbstractContinuousSystem}, tspan;
     if got_homogeneous
         reach_homog_BOX!(F, Ω0, Φ, NSTEPS, δ, X, recursive, Δt0)
     else
-        U = inputset(ivp_discr)
+        U = inputset(ivp)
         @assert isa(U, LazySet) "expected input of type `<:LazySet`, but got $(typeof(U))"
         # TODO: can we use support function evaluations for the input set?
         U = overapproximate(U, Hyperrectangle)

diff --git a/src/Algorithms/GLGM06/post.jl b/src/Algorithms/GLGM06/post.jl
@@ -1,13 +1,11 @@
-# continuous post for GLGM06 using Zonotope set representation
-function post(alg::GLGM06{N}, ivp::IVP{<:AbstractContinuousSystem}, tspan;
-              Δt0::TimeInterval=zeroI, kwargs...) where {N}
+# continuous post
+function post(alg::GLGM06, ivp::IVP{<:AbstractContinuousSystem}, tspan;
+              Δt0::TimeInterval=zeroI, kwargs...)
     δ = alg.δ
 
     NSTEPS = get(kwargs, :NSTEPS, compute_nsteps(δ, tspan))
 
     # normalize system to canonical form
-    # x' = Ax, x in X, or
-    # x' = Ax + u, x in X, u in U
     ivp_norm = _normalize(ivp)
 
     # homogenize the initial-value problem
@@ -21,6 +19,7 @@ function post(alg::GLGM06{N}, ivp::IVP{<:AbstractContinuousSystem}, tspan;
     return post(alg, ivp_discr, NSTEPS; Δt0=Δt0, kwargs...)
 end
 
+# discrete post
 function post(alg::GLGM06{N}, ivp::IVP{<:AbstractDiscreteSystem}, NSTEPS=nothing;
               Δt0::TimeInterval=zeroI, kwargs...) where {N}
     @unpack δ, approx_model, max_order, static, dim, ngens,

diff --git a/src/Algorithms/INT/post.jl b/src/Algorithms/INT/post.jl
@@ -1,11 +1,11 @@
-# this algorithms assumes that the initial-value problem is one-dimensional
-function post(alg::INT{N}, ivp::IVP{<:AbstractContinuousSystem}, tspan;
-              Δt0::TimeInterval=zeroI, kwargs...) where {N}
+# continuous post
+function post(alg::INT, ivp::IVP{<:AbstractContinuousSystem}, tspan;
+              Δt0::TimeInterval=zeroI, kwargs...)
     n = statedim(ivp)
     n == 1 || throw(ArgumentError("this algorithm applies to one-dimensional " *
                                   "systems, but this initial-value problem is $n-dimensional"))
 
-    @unpack δ, approx_model = alg
+    δ = alg.δ
 
     NSTEPS = get(kwargs, :NSTEPS, compute_nsteps(δ, tspan))
 
@@ -18,16 +18,33 @@ function post(alg::INT{N}, ivp::IVP{<:AbstractContinuousSystem}, tspan;
     end
 
     # discretize system
-    ivp_discr = discretize(ivp_norm, δ, approx_model)
-    Φ = state_matrix(ivp_discr)[1, 1]
-    Ω0 = initial_state(ivp_discr)
-    X = stateset(ivp_discr)
+    ivp_discr = discretize(ivp_norm, δ, alg.approx_model)
+
+    return post(alg, ivp_discr, NSTEPS; Δt0=Δt0, kwargs...)
+end
+
+# discrete post
+function post(alg::INT{N}, ivp::IVP{<:AbstractDiscreteSystem}, NSTEPS=nothing;
+              Δt0::TimeInterval=zeroI, kwargs...) where {N}
+    @unpack δ, approx_model = alg
+
+    if isnothing(NSTEPS)
+        if haskey(kwargs, :NSTEPS)
+            NSTEPS = kwargs[:NSTEPS]
+        else
+            throw(ArgumentError("`NSTEPS` not specified"))
+        end
+    end
+
+    Φ = state_matrix(ivp)[1, 1]
+    Ω0 = initial_state(ivp)
+    X = stateset(ivp)
 
     # this algorithm requires Ω0 to be an interval
     Ω0 = overapproximate(Ω0, Interval)
 
     # true <=> there is no input, i.e. the system is of the form x' = Ax, x ∈ X
-    got_homogeneous = !hasinput(ivp_discr)
+    got_homogeneous = !hasinput(ivp)
 
     # preallocate output flowpipe
     IT = typeof(Ω0)
@@ -37,7 +54,7 @@ function post(alg::INT{N}, ivp::IVP{<:AbstractContinuousSystem}, tspan;
     if got_homogeneous
         reach_homog_INT!(F, Ω0, Φ, NSTEPS, δ, X, Δt0)
     else
-        U = inputset(ivp_discr)
+        U = inputset(ivp)
         @assert isa(U, LazySet)
         U = overapproximate(U, Interval) # TODO guarantee this on the discretization?
         reach_inhomog_INT!(F, Ω0, Φ, NSTEPS, δ, X, U, Δt0)

diff --git a/src/Algorithms/LGG09/post.jl b/src/Algorithms/LGG09/post.jl
@@ -1,8 +1,10 @@
-function post(alg::LGG09{N,AM,VN,TN}, ivp::IVP{<:AbstractContinuousSystem}, tspan;
-              Δt0::TimeInterval=zeroI, kwargs...) where {N,AM,VN,TN}
-    @unpack δ, approx_model, template, static, threaded, vars = alg
+# continuous post
+function post(alg::LGG09, ivp::IVP{<:AbstractContinuousSystem}, tspan;
+              Δt0::TimeInterval=zeroI, kwargs...)
+    δ = alg.δ
 
     # dimension check
+    template = alg.template
     @assert statedim(ivp) == dim(template) "the problems' dimension $(statedim(ivp)) " *
                                            "doesn't match the dimension of the template directions, $(dim(template))"
 
@@ -17,10 +19,27 @@ function post(alg::LGG09{N,AM,VN,TN}, ivp::IVP{<:AbstractContinuousSystem}, tspa
     end
 
     # discretize system
-    ivp_discr = discretize(ivp_norm, δ, approx_model)
-    Φ = state_matrix(ivp_discr)
-    Ω₀ = initial_state(ivp_discr)
-    X = stateset(ivp_discr)
+    ivp_discr = discretize(ivp_norm, δ, alg.approx_model)
+
+    return post(alg, ivp_discr, NSTEPS; Δt0=Δt0, kwargs...)
+end
+
+# discrete post
+function post(alg::LGG09{N,AM,VN,TN}, ivp::IVP{<:AbstractDiscreteSystem}, NSTEPS=nothing;
+              Δt0::TimeInterval=zeroI, kwargs...) where {N,AM,VN,TN}
+    @unpack δ, approx_model, template, static, threaded, vars = alg
+
+    if isnothing(NSTEPS)
+        if haskey(kwargs, :NSTEPS)
+            NSTEPS = kwargs[:NSTEPS]
+        else
+            throw(ArgumentError("`NSTEPS` not specified"))
+        end
+    end
+
+    Φ = state_matrix(ivp)
+    Ω₀ = initial_state(ivp)
+    X = stateset(ivp)
 
     if alg.sparse # ad-hoc conversion of Φ to sparse representation
         Φ = sparse(Φ)
@@ -29,7 +48,7 @@ function post(alg::LGG09{N,AM,VN,TN}, ivp::IVP{<:AbstractContinuousSystem}, tspa
     cacheval = Val(alg.cache)
 
     # true <=> there is no input, i.e. the system is of the form x' = Ax, x ∈ X
-    got_homogeneous = !hasinput(ivp_discr)
+    got_homogeneous = !hasinput(ivp)
 
     # NOTE: option :static currently ignored!
 
@@ -41,7 +60,7 @@ function post(alg::LGG09{N,AM,VN,TN}, ivp::IVP{<:AbstractContinuousSystem}, tspa
     if got_homogeneous
         ρℓ = reach_homog_LGG09!(F, template, Ω₀, Φ, NSTEPS, δ, X, Δt0, cacheval, Val(alg.threaded))
     else
-        U = inputset(ivp_discr)
+        U = inputset(ivp)
         @assert isa(U, LazySet)
         ρℓ = reach_inhomog_LGG09!(F, template, Ω₀, Φ, NSTEPS, δ, X, U, Δt0, cacheval,
                                   Val(alg.threaded))