Outsource discretization methods to modules

docs/src/man/algorithms/ORBIT.md

@@ -44,13 +44,17 @@ Before considering the computation of $\Phi$ and $\Phi_1$, note that the method
 x(k\delta) = \Phi^k x_0 + \sum_{i=0}^{k-1} \Phi^i v_{k-i},\qquad k = 1,\ldots, N
 ```
 
-The matrix $\Phi = e^{A\delta}$ can be evaluated in different ways, using the function [`_exp`](@ref):
+The matrix $\Phi = e^{A\delta}$ can be evaluated in different ways, using the function
+[`ReachabilityAnalysis.Exponentiation._exp`](@ref):
 
-(1) `method=:base` uses Julia's built-in implementation (if `method = :base`),
+(1) `method = :base` uses Julia's built-in implementation (if `method = :base`),
 
 (2) `method = :lazy` uses a lazy wrapper of the matrix exponential which is then evaluated using Krylov subspace methods.
 
-Method (1) is the default method. Method (2) is particularly useful to work with very large and sparse matrices (e.g. typically of order `n > 2000`). Evaluation of $\Phi_1(u, \delta)$ is available through the function [`Φ₁`](@ref). Two implementations are available:
+Method (1) is the default method. Method (2) is particularly useful to work with
+very large and sparse matrices (e.g. typically of order `n > 2000`). Evaluation
+of $\Phi_1(u, \delta)$ is available through the function
+[`ReachabilityAnalysis.Exponentiation.Φ₁`](@ref). Two implementations are available:
 
 (1) If the coefficients matrix $A$ is invertible, then the integral is equivalent to computing $A^{-1}(e^{A\delta} - I)$.
 

docs/src/tutorials/linear_methods/dense_time.md

@@ -123,14 +123,15 @@ $\delta \to 0$ (in Hausdorff distance).
 ```
 
 Algorithms implementing conservative time discretization can be used from the
-[`discretize(ivp::IVP, δ, alg::AbstractApproximationModel)`](@ref) function.
+[`discretize(ivp::IVP, δ, alg::ReachabilityAnalysis.DiscretizationModule.AbstractApproximationModel)`](@ref)
+function.
 Set-based conservative discretization of a continuous-time initial value problem
 into a discrete-time problem.
 This function receives three inputs: the initial value problem (``ivp`) for a
 linear ODE in canonical form, (e.g. the system returned by `normalize`);
 the step-size (`δ`), and the algorithm (`alg`) used to compute the approximation model.
-Do `subtypes(ReachabilityAnalysis.AbstractApproximationModel)` to see the
-available approximation models.
+Do `subtypes(ReachabilityAnalysis.DiscretizationModule.AbstractApproximationModel)`
+to see the available approximation models.
 
 The output of a discretization is a new initial value problem of a discrete system.
 Different approximation algorithms and their respective options are described

src/Algorithms/A20/A20.jl

@@ -1,3 +1,5 @@
+using ..Overapproximate: _convert_or_overapproximate
+
 """
     A20{N} <: AbstractContinuousPost
 

src/Algorithms/ASB07/ASB07.jl

@@ -1,3 +1,5 @@
+using ..Overapproximate: _convert_or_overapproximate, _split, _overapproximate_interval_linear_map
+
 """
     ASB07{N, AM, RM, S, R} <: AbstractContinuousPost
 

src/Algorithms/GLGM06/GLGM06.jl

@@ -1,3 +1,5 @@
+using ..Overapproximate: _convert_or_overapproximate
+
 """
     GLGM06{N, AM, S, D, NG, P, RM} <: AbstractContinuousPost
 

src/Algorithms/ORBIT/ORBIT.jl

@@ -11,8 +11,8 @@ for singleton initial conditions and proper input sets.
 
 The type fields are:
 
-- `N`   -- number type of the step-size
-- `AM`  -- type of the approximation model
+- `N`  -- number type of the step-size
+- `AM` -- type of the approximation model
 
 ## Notes
 

src/Continuous/normalization.jl

@@ -39,10 +39,6 @@ end
 
 *(M::AbstractMatrix, input::ConstantInput) = ConstantInput(M * input.U)
 
-# convenience functions
-next_set(inputs::ConstantInput) = collect(nextinput(inputs, 1))[1]
-next_set(inputs::AbstractInput, state::Int64) = collect(nextinput(inputs, state))[1]
-
 """
     add_dimension(A::AbstractMatrix, m=1)
 
@@ -56,12 +52,14 @@ Append one or more zero rows and columns to a matrix.
 ### Examples
 
 ```jldoctest add_dimension_test
+julia> using ReachabilityAnalysis: add_dimension
+
 julia> A = [0.4 0.25; 0.46 -0.67]
 2×2 Matrix{Float64}:
  0.4    0.25
  0.46  -0.67
 
-julia> ReachabilityAnalysis.add_dimension(A)
+julia> add_dimension(A)
 3×3 Matrix{Float64}:
  0.4    0.25  0.0
  0.46  -0.67  0.0
@@ -70,7 +68,7 @@ julia> ReachabilityAnalysis.add_dimension(A)
 To append more than one zero row-column, use the second argument `m`:
 
 ```jldoctest add_dimension_test
-julia> ReachabilityAnalysis.add_dimension(A, 2)
+julia> add_dimension(A, 2)
 4×4 Matrix{Float64}:
  0.4    0.25  0.0  0.0
  0.46  -0.67  0.0  0.0
@@ -96,19 +94,21 @@ Adds an extra dimension to a LazySet through a Cartesian product.
 ### Examples
 
 ```jldoctest add_dimension_set
+julia> using ReachabilityAnalysis: add_dimension
+
 julia> X = BallInf(ones(9), 0.5);
 
 julia> dim(X)
 9
 
-julia> Xext = ReachabilityAnalysis.add_dimension(X);
+julia> Xext = add_dimension(X);
 
 julia> dim(Xext)
 10
 
 julia> X = ZeroSet(4);
 
-julia> dim(ReachabilityAnalysis.add_dimension(X))
+julia> dim(add_dimension(X))
 5
 
 julia> typeof(X)
@@ -118,7 +118,7 @@ ZeroSet{Float64}
 More than one dimension can be added passing the second argument:
 
 ```jldoctest add_dimension_set
-julia> Xext = ReachabilityAnalysis.add_dimension(BallInf(zeros(10), 0.1), 4);
+julia> Xext = add_dimension(BallInf(zeros(10), 0.1), 4);
 
 julia> dim(Xext)
 14
@@ -152,13 +152,15 @@ Adds an extra dimension to a initial-value system.
 ```jldoctest add_dimension_cont_sys
 julia> using MathematicalSystems, SparseArrays
 
+julia> using ReachabilityAnalysis: add_dimension
+
 julia> A = sprandn(3, 3, 0.5);
 
 julia> X0 = BallInf(zeros(3), 1.0);
 
 julia> s = InitialValueProblem(LinearContinuousSystem(A), X0);
 
-julia> sext = ReachabilityAnalysis.add_dimension(s);
+julia> sext = add_dimension(s);
 
 julia> statedim(sext)
 4
@@ -169,16 +171,18 @@ If there is an input set, it is also extended:
 ```jldoctest add_dimension_cont_sys
 julia> using LinearAlgebra
 
+julia> using ReachabilityAnalysis.DiscretizationModule: next_set
+
 julia> U = ConstantInput(Ball2(ones(3), 0.1));
 
 julia> s = InitialValueProblem(ConstrainedLinearControlContinuousSystem(A, Matrix(1.0I, size(A)), nothing, U), X0);
 
-julia> sext = ReachabilityAnalysis.add_dimension(s);
+julia> sext = add_dimension(s);
 
 julia> statedim(sext)
 4
 
-julia> dim(ReachabilityAnalysis.next_set(inputset(sext)))
+julia> dim(next_set(inputset(sext)))
 4
 ```
 
@@ -191,24 +195,24 @@ julia> U = VaryingInput([Ball2(ones(3), 0.1 * i) for i in 1:3]);
 
 julia> s = InitialValueProblem(ConstrainedLinearControlContinuousSystem(A, Matrix(1.0I, size(A)), nothing, U), X0);
 
-julia> sext = ReachabilityAnalysis.add_dimension(s);
+julia> sext = add_dimension(s);
 
 julia> statedim(sext)
 4
 
-julia> dim(ReachabilityAnalysis.next_set(inputset(sext), 1))
+julia> dim(next_set(inputset(sext), 1))
 4
 ```
 
 Extending a varying input set with more than one extra dimension:
 
 ```jldoctest add_dimension_cont_sys
-julia> sext = ReachabilityAnalysis.add_dimension(s, 7);
+julia> sext = add_dimension(s, 7);
 
 julia> statedim(sext)
 10
 
-julia> dim(ReachabilityAnalysis.next_set(inputset(sext), 1))
+julia> dim(next_set(inputset(sext), 1))
 10
 ```
 """

src/Discretization/ApplySetops.jl

@@ -0,0 +1,70 @@
+# =========================================================================
+# Alternatives to apply the set operation depending on the desired output
+# =========================================================================
+module ApplySetops
+
+using ..DiscretizationModule: AbstractApproximationModel, hasbackend
+using LazySets
+using LazySets.Approximations: AbstractDirections
+
+export _apply_setops
+
+function _apply_setops(X, alg::AbstractApproximationModel)
+    if hasbackend(alg)
+        _apply_setops(X, alg.setops, alg.backend)
+    else
+        _apply_setops(X, alg.setops)
+    end
+end
+
+_apply_setops(X::LazySet, ::Val{:lazy}) = X  # no-op
+_apply_setops(X::LazySet, ::Val{:concrete}) = concretize(X) # concrete set
+_apply_setops(X, template::AbstractDirections) = overapproximate(X, template) # template oa
+_apply_setops(X::LazySet, ::Val{:box}) = box_approximation(X) # box oa
+
+_apply_setops(M::AbstractMatrix, X::LazySet, ::Val{:lazy}) = M * X
+_apply_setops(M::AbstractMatrix, X::LazySet, ::Val{:concrete}) = linear_map(M, X)
+
+# evantually we should use concretize, but requires fast fallback operations in 2D
+# such as Minkowski sum not yet available
+function _apply_setops(X::ConvexHull{N,AT,MS}, ::Val{:vrep},
+                       backend=nothing) where {N,
+                                               AT<:AbstractPolytope{N},
+                                               LM<:LinearMap{N,AT,N},
+                                               MS<:MinkowskiSum{N,LM}}
+    n = dim(X)
+    VT = n == 2 ? VPolygon : VPolytope
+
+    # CH(A, B) := CH(X₀, ΦX₀ ⊕ E₊)
+    A = X.X
+    B = X.Y
+    X₀ = convert(VT, A)
+
+    if n == 2
+        ΦX₀ = convert(VT, B.X)
+        E₊ = convert(VT, B.Y)
+        out = convex_hull(X₀, minkowski_sum(ΦX₀, E₊))
+    else
+        # generic conversion to VPolytope is missing, see LazySets#2467
+        ΦX₀ = VPolytope(vertices_list(B.X; prune=false))
+        E₊ = convert(VT, B.Y)
+        aux = minkowski_sum(ΦX₀, E₊; apply_convex_hull=false)
+        out = convex_hull(X₀, aux; backend=backend)
+    end
+
+    return out
+end
+
+# give X = CH(X₀, ΦX₀ ⊕ E₊), return a zonotope overapproximation
+function _apply_setops(X::ConvexHull{N,AT,MS}, ::Val{:zono},
+                       backend=nothing) where {N,
+                                               AT<:AbstractZonotope{N},
+                                               LM<:LinearMap{N,AT,N},
+                                               MS<:MinkowskiSum{N,LM}}
+    # CH(A, B) := CH(X₀, ΦX₀ ⊕ E₊)
+    A = X.X
+    B = X.Y
+    return overapproximate(CH(A, concretize(B)), Zonotope)
+end
+
+end  # module

src/Discretization/BackwardModule.jl

@@ -4,7 +4,7 @@
 module BackwardModule
 
 using ..DiscretizationModule
-using ..Exponentiation: _alias
+using ..Exponentiation: _alias, BaseExp
 
 export Backward
 
@@ -55,11 +55,12 @@ function Base.show(io::IO, alg::Backward)
     print(io, "    - set operations method: $(alg.setops)\n")
     print(io, "    - symmetric interval hull method: $(alg.sih)\n")
     print(io, "    - invertibility assumption: $(alg.inv)\n")
-    return print(io, "    - polyhedral computations backend: $(alg.backend)\n")
+    print(io, "    - polyhedral computations backend: $(alg.backend)\n")
+    return nothing
 end
 
-Base.show(io::IO, m::MIME"text/plain", alg::Backward) = print(io, alg)
+Base.show(io::IO, ::MIME"text/plain", alg::Backward) = print(io, alg)
 
 # TODO: add corresponding `discrete` methods <<<<<
 
-end
+end  # module

src/Discretization/CorrectionHull.jl → src/Discretization/CorrectionHullModule.jl

@@ -1,14 +1,26 @@
 # ===============================================================
 # Discretize using the correction hull of the matrix exponential
 # ===============================================================
-
+module CorrectionHullModule
+
+using ..DiscretizationModule
+using ..Exponentiation: _exp, _alias, IntervalExpAlg, interval_matrix
+using ..Overapproximate: _convert_or_overapproximate, _overapproximate
+using IntervalMatrices: IntervalMatrix, correction_hull, input_correction, _exp_remainder
+using MathematicalSystems
+using LazySets
+using LazySets: LinearMap
+import IntervalArithmetic as IA
+using LinearAlgebra
 using Reexport
 
-using IntervalMatrices: correction_hull, input_correction
+export CorrectionHull
 
-using ..Exponentiation: _exp, _alias
 @reexport import ..DiscretizationModule: discretize
 
+const CLCS = ConstrainedLinearContinuousSystem
+const CLCCS = ConstrainedLinearControlContinuousSystem
+
 """
     CorrectionHull{EM} <: AbstractApproximationModel
 
@@ -47,10 +59,11 @@ end
 function Base.show(io::IO, alg::CorrectionHull)
     print(io, "`CorrectionHull` approximation model with:\n")
     print(io, "    - exponentiation method: $(alg.exp)\n")
-    return print(io, "    - order: $(alg.order)\n")
+    print(io, "    - order: $(alg.order)\n")
+    return nothing
 end
 
-Base.show(io::IO, m::MIME"text/plain", alg::CorrectionHull) = print(io, alg)
+Base.show(io::IO, ::MIME"text/plain", alg::CorrectionHull) = print(io, alg)
 
 # -----------------------------------------------------------------
 # Correction hull: homogeneous case x' = Ax, x in X
@@ -87,7 +100,7 @@ end
 
 # F(δ) without the E(δ) correction term
 function _correction_hull_without_E(A, δ, p)
-    timeint(δ, i) = interval((i^(-i / (i - 1)) - i^(-1 / (i - 1))) * δ^i, 0)
+    timeint(δ, i) = IA.interval((i^(-i / (i - 1)) - i^(-1 / (i - 1))) * δ^i, 0)
     F = sum(map(x -> timeint(δ, i) * x, A^i / factorial(i)) for i in 2:p)
     return IntervalMatrix(F)
 end
@@ -180,7 +193,7 @@ function _Cδ(A, δ, order)
     n = size(A, 1)
     A² = A * A
     if isa(A, IntervalMatrix)
-        Iδ = IntervalMatrix(Diagonal(fill(IntervalArithmetic.interval(δ), n)))
+        Iδ = IntervalMatrix(Diagonal(fill(IA.interval(δ), n)))
     else
         Iδ = Matrix(δ * I, n, n)
     end
@@ -200,6 +213,8 @@ function _Cδ(A, δ, order)
             M += (δⁱ⁺¹ / αᵢ₊₁) * Aⁱ
         end
     end
-    E = IntervalMatrices._exp_remainder(A, δ, order; n=n)
+    E = _exp_remainder(A, δ, order; n=n)
     return M + E * δ
 end
+
+end  # module